EPJ manuscript No.
(will be inserted by the editor)
Aspects of stochastic resonance in reaction–diffusion
systems: The nonequilibrium-potential approach
Horacio S. Wio1,a and Roberto R. Deza2
1 Instituto de F́ısica de Cantabria, Universidad de Cantabria and CSIC
E-39005 Santander, Spain
2 Departamento de F́ısica, FCEyN, Universidad Nacional de Mar del Plata
Deán Funes 3350, 7600 Mar del Plata, Argentina
Abstract. We analyze several aspects of the phenomenon of stochastic resonance
in reaction–diffusion systems, exploiting the nonequilibrium potential’s frame-
work. The generalization of this formalism (sketched in the appendix) to extended
systems is first carried out in the context of a simplified scalar model, for which sta-
tionary patterns can be found analytically. We first show how system-size stochas-
tic resonance arises naturally in this framework, and then how the phenomenon
of array-enhanced stochastic resonance can be further enhanced by letting the
diffusion coefficient depend on the field. A yet less trivial generalization is ex-
emplified by a stylized version of the FitzHugh–Nagumo system, a paradigm of
the activator–inhibitor class. After discussing for this system the second aspect
enumerated above, we derive from it—through an adiabatic-like elimination of
the inhibitor field—an effective scalar model that includes a nonlocal contribu-
tion. Studying the role played by the range of the nonlocal kernel and its effect
on stochastic resonance, we find an optimal range that maximizes the system’s
response.
1 Introduction
Stochastic resonance (SR) is nowadays a paradigm of the constructive effects of fluctuations on
nonlinear systems [1,2]. Sketchily, the phenomenon occurs whenever the Kramers’ rate for the
transition between attractors matches the typical frequency of a signal which is incapable by
itself to trigger that transition (i.e. it is subthreshold). Whereas several measures of SR can be
defined [the signal-to-noise ratio (SNR) and the spectral amplification factor (SAF) being the
main ones], theoretical analysis is usually carried on in terms of the two-state approximation
[1]. Since its discovery a quarter of century ago—and besides exploring related phenomena
like e.g., coherence resonance [3]—interest has gradually shifted towards increasingly complex
systems, networks and nonlinear media being the main directions. Instances of this trend are
the experiments carried out to explore the role of SR in sensory and other biological functions
[4], and experiments in chemical systems [5].
Our concern throughout this review will be with nonlinear media that can be described
as reaction–diffusion (RD) systems, namely those that can be thought of as a collection of
diffusively coupled nonlinear units. The possibility of enhancing the system’s response through
the coupling of those units [6,7,8,9,10,11,12] has been among the issues explored during the
last decade, together with the “naturalness” problem (how does nature manage to make the
a e-mail: wio@ifca.unican.es
2 Will be inserted by the editor
system’s response less dependent on a fine tuning of the noise intensity) or that of searching
for different ways to control the phenomenon [13,14].
In dissipative dynamical systems, the very notion of Lyapunov’s function is as useful as
that of attractor itself. Even when often it cannot be explicitly computed (because integrability
conditions are not readily met), it allows for picturesque reasoning in terms of “energy land-
scapes” or “attraction basins”. When those dynamical systems are submitted to forces that can
be modeled (à la Langevin) as stochastic, a new meaning—statistical in nature, akin to the
notion of “free energy”—is added to the picture (in fact, it is worth mentioning that the first
function known to have the Lyapunov property was Boltzmann’s H-function). Moreover, even
for vanishing noise intensity, the very existence of stochastic terms (a “transport matrix”) can
render the system well conditioned regarding integrability conditions. That was the rationale
behind the definition of nonequilibrium potential (NEP) [16,17,18,19], two approaches to which
are described in the Appendix. Such NEP is a special Lyapunov’s function of the associated
deterministic system, which for nonequilibrium systems plays a role similar to that played by
a thermodynamic potential in equilibrium thermodynamics [16]. It is closely related to the sta-
tionary solution of the system’s Fokker–Planck equation, and characterizes the global properties
of the dynamics: attractors, linear and relative stability of these attractors, height of the barri-
ers separating attraction basins. In addition, it allows to evaluate the transition rates among
the different attractors [16,17,18,19,20]. Regarding the problem of SR in extended systems, it
was shown that the knowledge of the NEP allows to obtain a rather complete picture of the
behavior of the output signal-to-noise ratio (SNR). The novelty with nonequilibrium extended
systems is that even pointlike attractors in the medium’s infinite-dimensional phase-space can
be nontrivial field configurations (real-space patterns).
In a series of recent papers we have studied the SR phenomenon for the transitions between
two different patterns [8,9,10,11,12,15], exploiting the concept of nonequilibrium potential. In
this review we discuss some recent results concerning different aspects of SR in RD systems. In
Sec. 2 we discuss the phenomenon of system-size stochastic resonance (SSSR), and show how
can it be analyzed and understood within a NEP framework [21]. In Sec. 3, after reviewing a
recent study on the enhancement of the SNR found for a scalar system with density-dependent
diffusivity [12], we discuss its extension [22] to an array of FitzHugh–Nagumo units [23]. In
Sec. 4—through an adiabatic-like elimination of the inhibitor field in an activator–inhibitor
system—an effective scalar system with a nonlocal term is derived, and the role of the local and
nonlocal interactions on the SR response studied. The main conclusions are finally summarized
in Sec. 5.
2 System-size stochastic resonance
Recent studies on biological models of the Hodgkin–Huxley type [24,25] have shown that ion
concentrations along cell membranes display intrinsic SR-like phenomena as the number of ion
channels is varied. A related result [26] shows that even in the absence of external forcing, the
regularity of the collective firing of a set of coupled excitable FitzHugh–Nagumo units is optimal
for a given number of elements. From a physics point of view, the same phenomenon—called
system-size stochastic resonance (SSSR)—has also been found in an Ising model as well as in
a set of globally coupled units described by a φ4 theory [27]. It has been even shown to arise in
opinion formation models [28].
Since the SSSR phenomenon is peculiar to extended systems, there is an obvious interest
in describing it within a NEP framework, that offers a very general framework for the study
of the dependence of SR and related phenomena on any of the system’s parameters. Here we
discuss in some detail a one-component (“scalar”) RD system [21], and briefly refer to other
cases analyzed in [27] and [29].
2.1 Review of the scalar model
For one-variable dynamical systems, the Lyapunov’s function can always be found by quadra-
ture. This property can be readily translated to scalar RD systems: the Lyapunov’s functional
Will be inserted by the editor 3
F [φ] fulfills the “potential” condition ∂tφ(y, t) = −δF [φ]/δφ(y, t), where δ/δφ(y, t) indicates a
functional derivative. This is also the NEP for a scalar transport matrix (i.e. a multiple of the
unit matrix in the medium’s infinite-dimensional phase space).
The specific model we shall focus on here has a piecewise linear reaction term, that mimics
general bistable RD models [23], e.g. those with a cubic-like reaction term. In the following,
we shall exploit some of the results on the influence of general (partially reflective or albedo)
boundary conditions found in [30], as well as previous studies of the NEP [17] and of SR
[8,9,11,15]. The particular dimensionless form of the deterministic model we start with is [30,8,9]
− φ+ φhΘ(φ− φc), (1)
where Θ(z) is the Heaviside step function, φc is the value at which the piecewise-linear “reaction
term” has the jump, and D is a phenomenological “diffusion coefficient”, not necessarily related
to γ in Sec. 2.2 [in Graham’s approach, Eq. (15), D/(∆y)2 would be the matrix elements Qνν±1
in a discretization of the Laplacian]. All the effects of the parameters that keep the system away
from equilibrium (such as the electric current in electrothermal devices like the ballast resistor
[23,30], or some external reactant concentration in chemical models) are now included in φc.
For the system to display a bistable behavior, it must be 0 < φc < φh.
We consider here the class of static structures φ(y) studied in [30]. They are even solutions
to the stationary (∂tφ = 0) version of Eq. (1) in the bounded domain [−yL, yL], with equal
albedo boundary conditions (b.c.) at both ends
∂φ(y, t)
y=±yL
= ∓k φ(±yL, t).
k > 0 is called the albedo parameter: the limit k → 0 yieds Neumann’s b.c. ∂yφ(y, t)|y=±yL = 0
and the k →∞ one, Dirichlet’s b.c. φ(±yL, t) = 0.
The explicit form of these static patterns is
φ(y) = φh

sinh(yc)ρ′
k, yL+y√
k, yL√
−yL ≤ y ≤ −yc,
1− cosh(y)ρ
yL−y±c√
k, yL√
−yc ≤ y ≤ yc,
sinh(yc)ρ′
k, yL−y√
k, yL√
yc ≤ y ≤ yL,
where ρ(k, ζ) = sinh(ζ) + k cosh(ζ) and ρ′(k, ζ) = ∂ρ/∂ζ. The coordinate values y±c at which
φ(yc) = φc are
y±c =
yL − ln
Z2 + 1− k2
1 + k
with Z =
. (3)
Each real solution y±c < yL to Eq. (3) represents a structure with a central “activated” zone
(φ > φc) and two lateral “resting” regions (φ < φc). Figure 5 in [30] displays the relation yc/yL
vs k, for several values of φc/φh.
Typical shapes of the arising patterns are shown in Fig. 1. Through a linear stability analysis
it has been shown [30] that the structure with the smallest “excited” region [that is with
yc = y+c , denoted by φu(y)] is unstable, whereas the other one [with yc = y
c , denoted by φ1(y)]
is linearly stable. The trivial homogeneous solution φ0(y) = 0 exists and is linearly stable for
any parameter set. These two linearly stable solutions (φ0 and φ1) are the only stable static
structures under albedo b.c. We will concentrate on the region of values of φc/φh, yL and k
where φ1 exists.
For the finite system with albedo b.c., the NEP is a functional of φ and a function of k, yL
and φc/φh. It has the expression [17]
F([φ], φc/φh, k, yL) =∫ yL
∫ φ(y,t)
[−φ′ + φhΘ(φ′ − φc)] dφ′ +
∂φ(y, t)
φ2(y, t)
4 Will be inserted by the editor
Fig. 1. Inhomogeneous static solutions to Eq. (1) for yL = D = φh = 1.0 and φc = 0.193 (highlighted
by the dashed horizontal line). Dash-dotted upper curve: φ1(y) for k = 1.0 (in this case, yc ≡ y−c > yL).
Lower curves: φ1(y) (solid line, for which +y
c is highlighted) and φu(y) (dotted line, for which −y+c is
highlighted) for k = 7.0.
Fig. 2. NEP F([φ], k, yL) evaluated at the inhomogeneous stationary solutions φ1(y) (lower branch)
and φu(y) (upper branch), as a function of: (a) system’s size yL, with k = 3.0; (b) albedo parameter k,
with yL = 1.2. The remaining parameters are D = 1.0, φc/φh = 0.193. The NEP for the homogeneous
stationary solution φ0(y) coincides with the horizontal axis.
When the NEP is evaluated at the inhomogeneous static solutions of Eq. (1) [Eqs. (2) and
(3)] it takes the explicit form [8,17]
Fu,1(φc/φh, k, yL) = F([φu,1], φc/φh, k, yL) (4)
= φ2h
−y±c (1− 2φcφh
+ sinh
y±c /
) ρ(k, (yL − y±c )/√D)
k, yL/
while at the trivial solution φ0 ≡ 0 it is F([φ0], yL) = F0 = 0.
Figure 2a depicts the nonequilibrium potential F([φ], yL) as a function of the system’s size
yL, keeping the albedo parameter k and the ratio φc/φh fixed. The curves correspond to the NEP
Fu,1(yL), whereas F0 coincides with the x–axis. Our focus is the “bistable zone” yL ' 0.72,
where φ1(y) exists. The unstable structure φu(y) is a saddle point for F [φ] (in the medium’s
infinite-dimensional phase space), so its NEP Fu(yL) > 0 (upper branch in Fig. 2a). On the
other hand, both φ1(y) (lower branch) and φ0 (x–axis) are local minima of the NEP.
One immediately notices that ∆F1(yL) ≡ Fu(yL)−F1(yL) is an (almost linearly) increasing
function of yL (this has a profound implication for SSSR, as we shall see below). Equation (3) has
real solutions only for yL ' 0.72. This corresponds to a supercritical saddle–node bifurcation,
at which both inhomogeneous structures pop up. Now, the most important feature in Fig. 2a is
that F1(yL) vanishes at a certain system’s size y∗L (≈ 1.0 for the given values of k and φc/φh).
At that point, the stable inhomogeneous structure φ1(y) and the trivial solution φ0(y) exchange
their relative stabilities.
Will be inserted by the editor 5
For completeness, in Fig. 2b we plot F([φ], k) for the same value of φc/φh as in Fig. 2a.
For the chosen value of yL it is always F1(k) < 0, and correspondingly there is no “stability
exchange” as a function of k. Also, the initially large ∆F1(k) decreases with k as Fu(k) and
F1(k) tend (for k →∞) to the values corresponding to Dirichlet’s b.c. [17].
2.2 Results for SSSR
By including an additive spatiotemporal noise source ξ(y, t) [15,33], Eq. (1) becomes a stochastic
partial differential equation for the random field φ(y, t). The simplest assumptions about ξ(y, t)
are that it is Gaussian, with zero mean and a correlation function given by 〈ξ(y, t)ξ(y′, t′)〉 =
2γ δ(t− t′)δ(y − y′), where γ denotes the noise strength.
As discussed in [8,9,11,15], known results for activation processes in multidimensional sys-
tems [31] allow us to estimate the activation rate using the following Kramers’-like expression
for the mean first-passage time for the transitions between attractors
〈τi〉 = τ0 exp
∆F i(yL)
where ∆F i(yL) = Fu(yL) − F i(yL), i = 0, 1. The prefactor τ0 is usually determined by the
curvature of F [φ] at its extrema. On one hand, it is typically several orders of magnitude smaller
than the average time 〈τ〉, while on the other it does not change significatively when varying the
system’s parameters around the “bistable point” y∗L, where F([φ0], y
L) = F([φ1], y
L). Hence we
may simplify the analysis by assuming here that τ0 is constant, and scale it out of our results.
The behavior of 〈τ〉 as a function of k and φc/φh has been shown in [8,9,17].
As done in [8], we now assume that the system is (adiabatically) subject to an external
harmonic variation of the parameter φc: φc(t) = φ0c + δφc cos(ωt) [9,15], and exploit the “two-
state approximation” [1] as in [9,11,15]. Such approximation reduces the whole dynamics on
the bistable potential landscape to one where the transitions occur only between the states
associated to the bottom of each well, hence the only dynamical contents resides in the transition
rates. Up to first order in the amplitude δφc (assumed to be small in order that the periodic
input be sub-threshold) the transition rates Wi adopt the form
µi ∓ αi
cos(ωt)
, (5)
where (at constant φh) µi ≈ exp[−∆F i(φ0c , yL)/γ] and αi ≈ µi (∂∆F i/∂φc|φ0c ), i = 0, 1. The
quantity inside parentheses can be obtained analytically using Eq. (4). These results allow to
calculate the autocorrelation function, the power spectrum density and finally the SNR, that
we indicate by R. The detailed calculation can be found in the appendix of [11]. Up to the
relevant order (the second) in the signal amplitude δφc, we obtain
4µ0µ1
(α0µ1 + α1µ0)2
µ0 + µ1
µ0 + µ1
Φ, (6)
where we have used the form of the αi to reduce the expression, and defined Φ = [2φh yc(yL)]2.
Figure 3 (left) is a plot of R as a function of the noise intensity γ for a fixed system’s length
yL, displaying the typical maximum that has become the fingerprint of the SR phenomenon.
In Fig. 3 (right), the roles of γ and yL are exchanged (R is plotted as a function of yL for fixed
γ). Such a response is the expected one for a system exhibiting SSSR. In both cases, the values
of k and φ0c/φh are kept fixed.
Within the NEP context and in this kind of systems, the phenomenon arises due to the
breakdown of the NEP’s symmetry. This means that (as shown in Fig. 2) when varying yL,
both attractors can exchange their relative stability. For yL = y∗L ≈ 1 both stable structures—
the inhomogeneous one φ1(y) and the trivial one φ0—have the same value for the NEP. For
yL < y
L, φ1(y) becomes less stable than φ0 so transitions from φ1(y) to φ0 are more frequent
6 Will be inserted by the editor
Fig. 3. Left: SNR vs the noise intensity γ, with system’s size yL = 1.1. Right: SNR vs yL, for γ = 0.1.
The remaining parameters are k = 3.0, D = 1.0, and φ0c/φh = 0.193.
(the barrier is lower) than in the reverse direction, thus reducing the system’s response. When
yL ∼ 0.72, φ1(y) and φu(y) coalesce and disappear, and the response is strictly zero (within
the linear response scheme implicit in the two-state approximation). When yL > y∗L, φ1(y)
becomes more stable than φ0, making now transitions from φ0 to φ1(y) more frequent than in
the reverse direction, and reducing again the system’s response. Clearly, the system’s response
has a maximum when both attractors have the same stability (yL = y∗L), and decays when
departing from that situation. Hence, for this system and within this framework, SSSR arises
as a particular case of the more general discussion done in [11]. It should not come as a surprise
to find an analogy with the mechanism of double stochastic coherence described in [32], where
the NEP’s symmetry is induced by (an additional, multiplicative) noise.
By comparing figures 2a and 3 it becomes apparent that the value of yL at which the SNR
has its maximum differs slightly from y∗L (where the crossing between F
1 and F0 takes place).
The origin of this discrepancy is the following: whereas on qualitative grounds we have argued
that the maximum of the SNR should be related to the potential being symmetric (both wells
having the same “energy”) [11], the exact condition is that the transition rates between both
wells be equal. In general, due to small differences between the curvatures at the bottom of each
well, those rates become equal for values of yL slightly different from the one at the symmetric
case. Although by adopting here a constant value of τ0 we have assumed equal curvatures, there
is still a difference in the values of the αi, since the ∂∆F i/∂φc|φ0c , i = 0, 1 are slightly different
(a fact reflected in the dependence of Φ on yL).
Additional light can be shed on the phenomenon when viewed from a different angle. Figure
4 is a plot of R as a function of k at fixed values of γ, yL and φ0c/φh. It exhibits a broad resonance
since for k not too large (indicating a high reflectiveness at the boundary or a reduced exchange
with the environment) R increases with k, whereas it slightly decreases for larger k values (where
the system’s boundaries become absorbent). An explanation of this behavior in terms of the
NEP has been given in [29]: as already observed, the NEP’s symmetry is broken for this value of
yL; moreover, whereas the lower branch in Fig. 2b goes rapidly towards the value corresponding
to Dirichlet’s b.c. the upper branch keeps increasing, thus degrading the SNR. In any case, the
fact that the resonance is broad indicates the robustness of the system’s response with regard
to k, a parameter that (together with γ) encodes the coupling with the environment.
We stress the fact that the NEP framework put forward in this review allows to study
SSSR between whole patterns. The explanation offered in [27] to the phenomenon resorted to a
collective variable X ≈ (1/N)
j=1 xj , and to the fact that the noise in the effective stochastic
differential equation for X scaled with size. In [29] it was shown that all the cases discussed in
[27] can be put within the same NEP framework than the above studied scalar model. In fact,
the aforementioned almost linear increasing dependence of ∆F1(yL) on yL can be interpreted
as a noise scaling with size. There are however situations where the NEP’s symmetry is retained
as the system’s size is varied. We may then speak of a genuinely noise-scaling SSSR, in contrast
to the cases that could be called NEP symmetry breaking SSSR [29].
Will be inserted by the editor 7
Fig. 4. SNR vs k for yL = 1.2, γ = 0.1 and φ
c/φh = 0.193.
3 Case of selective coupling
In this section we analyze SR in two extended systems with density-dependent diffusive-like cou-
pling: an extension of the scalar RD model considered in Sec. 2 [12], and an array of FitzHugh–
Nagumo [23] units.
3.1 Scalar model
Here we extend the one-component RD model discussed in Sec. 2 by letting the diffusive pa-
rameter D in Eq. (1) depend on the field φ(x, t). As a matter of fact, since in the ballast resistor
[23,30] the thermal conductivity is a function of the energy density, the resulting equation for
the temperature field includes a temperature-dependent diffusion coefficient in a natural way.
The form of the governing equation is now
∂tφ(x, t) = ∂x [D(φ)∂xφ] + f(φ) + ξ(x, t), (7)
with ξ(x, t) and f(φ) as in Sec. 2.
As it was done for the reaction term, a simple choice (that retains however the qualitative
features of the system) is to consider the following dependence of the diffusion term on the field
variable
D(φ) = D0[1 + hΘ(φ− φc)].
For simplicity, here we choose the same threshold φc for the reaction term and the diffusion
coefficient.
We assume the system to lie in a bounded domain [−L,L], with Dirichlet b.c. at both ends:
φ(±L, t) = 0. The form of the patterns is analogous to what has been obtained in Sec. 2, the
only difference being that in the present case dφ/dx|xc is discontinuous and the area of the
“activated” central zone depends on h.
As before, the indicated patterns are extrema of the NEP: the unstable pattern φu(x) is
a saddle-point of this functional, separating the attractors φ0(x) and φs(x). For the case of a
field-dependent diffusion coefficient D(φ(x, t)) as described by Eq. (7), the NEP reads [12]
F [φ] =
D(φ′)f(φ′) dφ′ +
Given that ∂tφ = −[1/D(φ)](δF/δφ), one finds dF/dt = −
(δF/δφ)2 dx ≤ 0, thus warranting
the Lyapunov’s functional property.
Whereas in Sec. 2 we kept φ0c constant and varied yL, we now vary instead φ
c at constant
L. Similarly as before, both linearly stable states have the same value of the NEP (i.e., they
are equally stable) at some value φ∗c of the threshold. The way F [φ] depends on φ0c resembles
the dependence on yL shown in Sec. 2, but now ∆Fs is an (almost linearly) decreasing function
8 Will be inserted by the editor
0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14
γ (arb. units)
-0.2 -0.1 0.0 0.1 0.2
h (arb.units)
 (arb. units)
0.90 0.95 1.00 1.05 1.10
Fig. 5. Left: SNR vs the noise intensity γ, for D0 = 1.0 and h = 0.0 (full line), −0.25 (dashed line)
and 0.25 (dotted line). Right: the maximum Rmax of the SNR curve as a function of the selectiveness
h of the coupling, for D0 = 0.9 (dashed line), 1.0 (full line) and 1.1 (dotted line). The arrows a and b
indicate the response gain due respectively to a homogeneous increase of the coupling and to a selective
one. The larger gain in the second case is apparent. The inset shows the dependence of Rmax on D0 for
h = −0.25 (lower line), 0.0 and 0.25 (upper line). The remaining parameters are L = 1.0, δφc = 0.01
and Ω = 0.01.
of φ0c and both inhomogeneous structures coalesce and disappear through a subcritical saddle–
node bifurcation. As in the previous case, we analyze only the neighborhood of φ∗c . We shall
moreover consider only the neighborhood of h = 0, where the main trends of the effect can be
captured.
Figure 5 (left) depicts the dependence of R on the noise intensity γ for three values of h,
each curve displaying the typical SR maximum. Figure 5 (right) is a plot of the value Rmax
of these maxima as a function of h. The dramatic increase of Rmax (several dB for a small
positive variation of h) is apparent, and shows the strong effect that the selective coupling (or
field-dependent diffusivity) has on the system’s response.
It must be noted that the only two approximations made in order to derive our results—
namely the Kramers-like expression and the two-level approximation used for the evaluation of
the correlation function [12]—break down for large positive values of h because for increasing
selectivity the curves of F [φ] vs φ0c shift towards the left [12], which in turn means that the
barrier separating the attractors at φ∗c tends to zero. This effect is basically the same as the one
discussed in Refs. [9,11] in connection with global diffusivity D0. It is also worth noting that
except for the two aforementioned approximations, all the previous results (e.g. the profiles
of the stationary patterns and the corresponding values of the nonequilibrium potential) are
analytically exact.
3.2 FitzHugh–Nagumo model
Here we study an array of FitzHugh–Nagumo [23] units, with a density-dependent (diffusive-
like) coupling. The NEP for this system was found within the excitable regime and for particular
values of the coupling strength [19]. In the general case, however, the form of the NEP has not
been found yet. Hence, we have resorted to a study based on numerical simulations, analyzing
the influence of different parameters on the system’s response. Nevertheless, the idea of the
existence of such a NEP has always underlied this study. The results show that the enhancement
of the SNR found for the scalar system [12] is robust, and that the indicated non-homogeneous
coupling could clearly contribute to enhance the SR phenomenon in more general situations.
We consider a simplified version of the FitzHugh–Nagumo model [12,17,23], which has been
useful for gaining qualitative insight into the excitable and oscillatory dynamics in neural and
chemical systems. It consist of two variables:
– a (fast) activator field u, that in the case of neural systems represents the voltage variable,
while in chemical systems represents the concentration of an autocatalytic species.
Will be inserted by the editor 9
– an inhibitor field v, associated (within a neural context) to the concentration of potassium
ions in the medium, and that in a general chemical reaction inhibits the generation of the
u species.
Instead of considering the usual cubic-like nonlinear form, we use a piecewise linear version
∂u(x, t)
Du(u)
+ f(u)− v + ξ(x, t) (8)
∂v(x, t)
Dv(v)
+ β u− α v, (9)
where f(u) = −u + Θ(u − φc), and ξ(x, t) is a δ–correlated white Gaussian noise, as before.
γ indicates the noise intensity and φc is the “discontinuity” point, at which the piecewise
linearized function f(u) presents a jump. The parameter � = τu/τv � 1 indicates the timescale
ratio between the (fast) activator and the inhibitor. Together with α and β, it is chosen to
correspond to the excitable regime. We consider Dirichlet b.c. at x = ±L. Although the results
are qualitatively the same as those that could appear considering the usual FitzHugh–Nagumo
equations, this simplified version allows us to compare directly with the previous analytical
results for this system [12,15].
As in [12], we assume that the diffusion coefficient Du(u) is not constant, but depends on
the field u according to Du(u) = D0u [1 + hΘ(u− φc)]. This form implies that the value of
Du(u) depends “selectively” on whether the field u fulfills u > φc or u < φc. D0u is the value of
the diffusion constant without such “selective” term, and h indicates the size of the difference
between the diffusion constants in both regions [if h = 0 then Du(u) = D0u constant]. Dv(v) is
the diffusion for the inhibitor v, that we assume to be homogeneously constant.
This system is known to exhibit two stable stationary patterns. One of them is u(x) = 0,
v(x) = 0, while the other is one with nonzero values. Further, we consider, as before, that
an external, periodic, signal enters into the system through the value of the threshold φc,
φc(t) = φ0c + δφ cos(ωt), where ω is the signal frequency, and δφ its intensity.
All the results were obtained through numerical simulations of the system. The continuous
version of the system indicated by Eqs. (8), (9) was transformed into a second-order spatially
discrete one
u̇i = Dui∆ui + (Dui+1 −Dui−1)(ui+1 + ui−1) + f(ui)− vi + ξi(t),
v̇i = Dv∆vi + β ui − α vi,
with ∆φi ≡ (φi−1 + φi+1 − 2φi). The extensive numerical simulations performed for a set of
equations were done exploiting Heun’s algorithm [33].
In this spatially-extended system there are different ways of measuring the overall system’s
response to the external signal. In particular, we show the evaluated output SNR in two different
ways (the units being given in dB):
– SNR for the middle element of the chain evaluated over the dynamical evolution of uN/2,
that we call SNR2 (however, having Dirichlet b.c., the local response depends on the distance
to the boundaries).
– In order to measure the overall response of the system to the external signal, we computed
the SNR as follows: We digitized the system’s dynamics to a dichotomic process s(t): At
time t the system has an associated value of s(t) = 1 (0) if the Hilbert distance to pattern
1 (0) is lower than to the other pattern. Stated in mathematical terms, we computed the
distance D2[·, ·] defined by
D2[f, g] =
dx [f(x)− g(x)]2
in L2([−L,L]), the Hilbert space of the real-valued functions in that interval. At time t, a
digitized process is computed by means of
s(t) =
1 if D2 [Pu1 (x), u(x, t)] < D2 [Pu0 (x), u(x, t)]
0 if D2 [Pu1 (x), u(x, t)] ≥ D2 [Pu0 (x), u(x, t)]
10 Will be inserted by the editor
Model Numerical
α β � φ0c δφ ω Dv ∆t N
0.3 0.4 0.03 0.52 0.4 5π/8 1.0 10−3 51
Table 1. Fixed parameters for the FitzHugh–Nagumo model with Du(u)
Fig. 6. SNR2 and SNRp vs the noise intensity γ, for D
u = 0.3 and h = −2.0 (+), −1.0 (4), 0.0 (�),
1.0 (�), and 2.0 (©). Both measures reveal a systematic enhancement of the SNR as h increases.
-2 0 2 4 6 8 10
-2 0 2 4 6 8 10
Fig. 7. SNR2 and SNRp vs the selectiveness of coupling h, for D
u = 0.3 and γ = 0.032 (©), 0.32 (�),
0.6 (�), 1.2 (4), and 3.2 (×).
We call this “global-like” measure SNRp.
The parameters kept fixed have been summarized in Table 1. The simulation was repeated 250
times for each parameter set, and the SNR was computed by recourse of the average power
spectral density.
Figure 6 depicts the results for the different SNR measures we have previously defined, as
functions of the noise intensity γ. For both measures it is apparent that there is an enhancement
of the response for h > 0, when compared with the h = 0 case, while for h < 0 the response is
smaller.
In Fig. 7 we show again the response’s measures, but now as functions of h. We have plotted
the maximum of each SNR curve for D0u = 0.3, and γ = 0.01, 0.1, and 0.3. It is clear that there
exists an optimal value of γ for which the response is largest. The rapid fall in the response for
h < 0 is also apparent.
In Fig. 8 we show the dependance of SNR on h, for different values of the diffusion which
depends on the activator density D0u. It is apparent that the response becomes larger when
the value of D0u is larger. However, as was discussed in [8,15], it is clear that for still larger
values of D0u, the symmetry of the underlying potential (that is the relative stability between
the attractors) is broken and the response finally falls down.
The previous figures clearly show that the response to the external signal grows with the
“selectiveness” of the coupling, showing the robustness of the phenomenon presented in [12,15].
Will be inserted by the editor 11
-2 0 2 4 6 8 10
-2 0 2 4 6 8 10
Fig. 8. SNR2 and SNRp vs h for γ = 0.032 and D
u = 0.0 (©), 0.1 (�), and 0.3 (�).
4 Nonlocal Interaction
Let us consider again a system like the one described by Eqs. (8), (9), but now assume that
Du and Dv are constant. In Ref. [18] it was assumed that the inhibitor-like field has a diffusive
transport behavior, and is fast enough that can be adiabatically eliminated, thus yielding an
effective scalar RD equation with a nonlocal term, characterized by a diffusive kernel G(x, x′).
After briefly reviewing the derivation of the NEP for that situation, we shall assume in this
section that the transport mechanism of the adiabatically eliminated inhibitor-like field is of
nondiffusive character, thus yielding a kernel H(x, x′) that is more localized in space, and with
a controllable interaction range.
Following Ref. [18], let the system be defined by
∂tu(x, t) = Du∂
xu(x, t)− u(x, t) +Θ[u(x, t)− a]− v(x, t),
�−1 ∂tv(x, t) = Dv∂
xv(x, t) + βu(x, t)− αv(x, t), (10)
where � was defined after Eqs. (8), (9), and let it be confined to the domain [−L,L], with
Dirichlet b.c. at both ends: u(±L, t) = v(±L, t) = 0. Contrarily to the standard hypothesis,
we now assume that the inhibitor is much faster than the activator (i.e. τv � τu). In the limit
�→∞, we can rewrite Eq. (10) as
∂tu(x, t) = Du∂
xu(x, t)− u(x, t) +Θ[u(x, t)− a]− v(x, t),
0 = Dv∂
xv(x, t) + βu(x, t)− αv(x, t).
In the last pair of equations we can eliminate the inhibitor (which is now slaved to the activator)
by solving the second equation using the Green’s function method
[−Dv∂2x + α]G(x, x
′) = δ(x− x′),
v(x) = β
dx′G(x, x′)u(x′),
where the Green’s function G(x, x′) is given by
G(x, x′) =
[sinh k(L− x′)/ sinh 2kL] sinh k(L+ x) x < x′,
[sinh k(L+ x′)/ sinh 2kL] sinh k(L− x) x > x′,
with k = (α/Dv)1/2. This slaving procedure reduces our system to a nonlocal equation for the
activator only, that has the form
∂u(x, t)
∂2u(x, t)
+ f(u)− β
G(x, x′)u(x′) dx′. (11)
12 Will be inserted by the editor
From this equation, and taking into account the symmetry of the Green’s function G(x, x′), we
can obtain the Lyapunov functional for this system, which has the form
F [u] =
f(w) dw +
dx′G(x, x′)u(x′)u(x)
. (12)
This spatial nonlocal term in the NEP takes into account the repulsion between activated
zones. When two activated zones come near each other, the exponential tails of the inhibitor
concentration overlap, increasing its concentration between both activated zones and creating
an effective repulsion between them. Hence the Green’s function plays the role of an exponential
screening between the activated zones. In Ref. [9] the knowledge of such NEP was exploited to
study SR on the system indicated by Eq. (11).
The starting point of our present analysis will be the effective, nonlocal and stochastic RD
equation for the real (activator-like) field φ(x, t), analogous to Eq. (11), defined in the one
dimensional domain x ∈ [−L,L] by
+ f(φ)− β
H(x, x′)φ(x′) dx′ + ξ(x, t),
where the diffusivityD is constant, and we assume a cubic nonlinear term f(φ) = φ (φ−b) (2−φ).
Here ξ(x, t) is an additive Gaussian white noise, as in the previous cases. As before, the system
is subject to Dirichlet b.c. φ(±L, t) = 0.
Similarly to Eq. (11), this system could be written in a variational form, with the func-
tional F [φ] given by Eq. (12). As anticipated, here we consider a nondiffusive kernel, with a
controllable interaction range. In order to keep our analysis simple we propose the following
H(x, x′) =
1/2; |x− x′| ≤ l
0; |x− x′| > l, (13)
which allows the analysis by just varying the interaction range 2l.
The new effective RD equation contains local and nonlocal couplings (corresponding to the
diffusive and the nonlocal contribution, respectively). The last one contains the nonlocal kernel
given by Eq. (13), with a variable range 2l, that corresponds to the interaction of the field at
points x′ ∈ [x − l, x + l]. However, such points will contribute if and only if they are inside
the domain [−L,L]. We are now in position to study the role played by the nonlocal kernel
(particularly by its range 2l) on the SR phenomenon.
As before, the SR between stationary solutions was investigated in terms of the two-state
approach (all the details about the procedure and the evaluation of the SNR can be found
in Refs. [11,15]). As usual, we subject our system to a weak external signal b = b0 + A(t) =
b0 +∆b cos(ωst), rocking the NEP [15]. In order to have a subthreshold signal, the amplitude
∆b should satisfy ∆b� b. We have chosen b0 as the value of b at which F [φs] = F [φ = 0] = 0
when β = 0.
Up to first order in the small amplitude ∆b, the transition rates Wi and the functions αi have
the form indicated in Eqs. (5), (6). But now Φ, that depends on the inhomogeneous attractor
φs, has the form
After fixing the length of the system L and the kernel range 2l we can use the above indicated
expressions to find the SNR, that shows the usual bell-shaped form of stochastic resonance as
a function of the noise intensity.
In Fig. 9 we show the dependence of the SNR on the kernel range 2l for a fixed value of 2L.
There is a nonmonotonic behavior in the system’s response against variation of 2l, that can be
explained by the following facts:
– On one hand, the transition rates are decreasing functions of the range 2l for fixed 2L.
Therefore, the ratio µ1µ2/(µ1 +µ2) in the expression for the SNR also reflects this behavior.
Will be inserted by the editor 13
– The other factor in this expression has a maximum for a kernel range that corresponds to
“first neighbor” sites, namely x = x′ ± l.
Fig. 9. SNR as a function of 2l for 2L fixed. The parameters are D = 0.6, 2L = 6.25, β = 0.02 and
b0 = 0.719123.
The maximum in the system’s response as a function of the kernel range is due to the interplay
between these two factors. From this analysis of the comparative weights of the local (diffusive)
and nonlocal terms contributing to the SR response, it is apparent that the range of such
nonlocal kernel has an optimum value yielding a maximum for the SNR [36].
5 Conclusions
We have discussed three different aspects of the phenomenon of stochastic resonance in reaction–
diffusion systems, within the nonequilibrium potential’s framework. In first place we have dis-
cussed system-size SR in a scalar model. Even though we have not shown the details here, it
has been also possible to also study other cases [29]. In particular, a model of globally coupled
nonlinear oscillators discussed in [27], showing that it can also be described within the NEP
framework, with SSSR arising through an “effective” scaling of the noise intensity with the
system’s size.
In second place we presented a study of SR in systems with a density-dependent (diffusive-
like) coupling. We initially discuss the case of a scalar system [12], and afterwards extent the
analysis to an array of FitzHugh–Nagumo units, with a field-dependent activator diffusion [22].
For the second system, when both diffusions are constant (that is: Du > 0 and Dv = 0), has
a known form of the NEP [15]. However, in the general case we have not been able to find
the form of the NEP (but the idea of such a NEP is always underlying our analysis) and have
to resort to an analysis based on numerical simulations. The result shows that the system’s
response is enhanced due to the particular form of the non-homogeneous coupling. From such
results, we can conclude that the phenomenon of enhancement of the SNR, due to a selectivity
in the coupling, initially found for a scalar system [12] is robust, and that the indicated non-
homogeneous coupling could clearly contribute to enhance the SR phenomenon in very general
systems.
Finally, we analyzed an activator-like field including a nonlocal contribution that arise
through an effective adiabatic elimination of an auxiliary (inhibitor-like) field. By exploiting the
knowledge of the nonequilibrium potential in such a case, we have analyzed the dependence of
the SNR on the nonlocal interaction kernel range, founding that there is an optimal value of the
kernel range, yielding a maximum in the system’s response, corresponding to a very localized
interaction.
The indicated results clearly show that the “nonequilibrium potential” (even if not known
in detail [10]) offers a very useful framework to analyze a wide spectrum of characteristics
associated to SR in spatially extended or coupled systems. For instance, within this framework,
14 Will be inserted by the editor
the phenomenon of SSSR looks—as other aspects of SR in extended systems [11]—as a natural
consequence of a breaking of the symmetry of the NEP [27].
Appendix: Brief review of the nonequilibrium potential scheme
Loosely speaking, the notion of NEP is an extension to nonequilibrium situations of that of
equilibrium thermodynamic potential. In order to introduce it, we consider a general system of
nonlinear stochastic equations (admitting the possibility of multiplicative noises)
q̇ν = Kν(q) + gνi (q) ξi(t), ν = 1, . . . , n; (14)
where repeated indices are summed over. Equation (14) is stated in the sense of Itô. The {ξi(t)},
i = 1, . . . ,m ≤ n are mutually independent sources of Gaussian white noise with typical strength
Graham’s approach
The Fokker–Planck equation corresponding to Eq. (14) takes the form
Kν(q)P +
∂qν ∂qµ
Qνµ(q)P (15)
where P (q, t; γ) is the probability density of observing q = (q1, . . . , qn) at time t for noise
intensity γ, and Qνµ(q) = gνi (q) g
i (q) is the matrix of transport coefficients of the system,
which is symmetric and non-negative. In the long time limit (t→∞), the solution of Eq. (15)
tends to the stationary distribution Pst(q). According to [16], the NEP Φ(q) associated to Eq.
(15) is defined by
Φ(q) = − lim
γ lnPst(q, γ). (16)
In other words
Pst(q) d
nq = Z(q) exp
+O(γ)
where Φ(q) is the NEP of the system and Z(q) is defined as the limit
lnZ(q) = lim
lnPst(q, γ) +
Here dΩq = dnq/
G(q) is the invariant volume element in the q-space and G(q) is the deter-
minant of the contravariant metric tensor (for the Euclidean metric it is G = 1). It was shown
[16] that Φ(q) is the solution of a Hamilton–Jacobi-like equation (HJE)
Kν(q)
Qνµ(q)
and Z(q) is the solution of a linear first-order partial differential equation depending on Φ(q)
(not shown here).
Equation (16) and the normalization condition ensure that Φ is bounded from below. Fur-
thermore, from the separation of the streaming velocity of the probability flow in the steady
state into conservative and dissipative parts, it follows that
dΦ(q)
= Kν(q)
∂Φ(q)
Qνµ(q)
i.e. Φ is a LF for the dynamics of the system when fluctuations are neglected. Under the
deterministic dynamics, q̇ν = Kν(q), Φ decreases monotonically and takes a minimum value on
attractors. In particular, Φ must be constant on all extended attractors (such as limit cycles or
strange attractors) [16].
Will be inserted by the editor 15
Ao’s approach
An alternative way to look into this problem is due to Ao [37]. Let us refer again to the system
in Eq. (14). Following [37], we introduce now the auxiliary matrix
H−1(q) = S(q) + A(q),
with S(q) a symmetric matrix while A(a) is an antisymmetric one. H−1(q) is now used to
rewrite the initial system as
H−1(q) q̇ = H−1(q) K(q, t) + H−1(q) ξ(q, t) = −∇Φ(q, t) + η(q, t),
where −∇Φ(q, t) = H−1(q) K(q, t), η(q, t) = H−1(q) ξ(q, t) and H−1(q) q̇ = −∇Φ(q, t) +
η(q, t). The new stochastic variables, η(q, t), fulfill
〈η(q, t)ηT (q, t′)〉 = 2 S(q) δ(t− t′) = 2 G−1(q)Q(q)[G−1(q)]T δ(t− t′),
that imposes a condition on the arbitrary definition of S(q) as we have
[S(q) + A(q)] Q(q) [S(q)−A(q)] = S(q).
As shown by Ao, the last equation also implies H(q) + HT (q) = 2 Q(q). We also have
H−1(q) K(q, t)
= ∇× [−∇Φ(q, t)] = 0. From the previous equations, we have in principle
all the needed conditions to determine H−1(q), and to obtain from it the potential Φ(q, t). The
interesting feature of this approach is that it resorts neither to Pst(q) nor to the small-noise
limit, thus being applicable in principle to more general situations.
Acknowledgements
The authors acknowledge the collaboration of B. von Haeften, S. Bouzat, C. J. Tessone, G. G.
Izús, M. Kuperman, S. Mangioni, A. Sánchez, F. Castelpoggi, in different aspects and/or stages
of this work. HSW thanks the European Commission for the award of a Marie Curie Chair at
the Universidad de Cantabria, Spain.
References
1. L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni, Rev. Mod. Phys. 70, 223 (1998).
2. T. Wellens, V. Shatokin, and A. Buchleitner, Rep. Prog. Phys. 67, 45 (2004).
3. A. S. Pikovsky and J. Kurths, Phys. Rev. Lett. 78, 775 (1997).
4. J. K. Douglas et al., Nature 365, 337 (1993); J. J. Collins et al., Nature 376, 236 (1995); S. M.
Bezrukov and I. Vodyanoy, Nature 378, 362 (1995).
5. A. Guderian, G. Dechert, K. Zeyer, and F. Schneider, J. Phys. Chem. 100, 4437 (1996); A. Förster,
M. Merget, and F. Schneider, J. Phys. Chem. 100, 4442 (1996); W. Hohmann, J. Müller, and F.
W. Schneider; J. Phys. Chem. 100, 5388 (1996).
6. J. F. Lindner et al., Phys. Rev. E 53, 2081 (1996).
7. A. R. Bulsara and G. Schmera, Phys. Rev. E 47, 3734 (1993); P. Jung, U. Behn, E. Pantazelou,
and F. Moss, Phys. Rev. A 46, R1709 (1992); P. Jung and G. Mayer-Kress, Phys. Rev. Lett. 74,
2130 (1995); J. F. Lindner et al., Phys. Rev. Lett. 75, 3 (1995); F. Marchesoni, L. Gammaitoni,
and A. R. Bulsara, Phys. Rev. Lett. 76, 2609 (1996).
8. H. S. Wio, Phys. Rev. E 54, R3075 (1996); H. S. Wio and F. Castelpoggi, Unsolved Problems of
Noise, Proc. Conf. UPoN’96, C. R. Doering, L. B. Kiss, and M. Schlesinger, Eds. (World Scientific,
Singapore, 1997), pg. 229; F. Castelpoggi and H. S. Wio, Europhys. Lett. 38, 91 (1997).
9. F. Castelpoggi and H. S. Wio, Phys. Rev. E 57, 5112 (1998).
10. H. S. Wio et al., Physica A 257, 275 (1998); M. Kuperman, H. S. Wio, G. Izús, and R. Deza, Phys.
Rev. E 57, 5122 (1998).
11. S. Bouzat and H. S. Wio, Phys. Rev. E 59, 5142 (1999).
16 Will be inserted by the editor
12. B. von Haeften, R. Deza, and H. S. Wio, Phys. Rev. Lett. 84, 404 (2000).
13. C. J. Tessone, H. S. Wio, and P. Hänggi, Phys. Rev. E 62, 4623 (2000).
14. M. A. Fuentes, R. Toral, and H. S. Wio, Physica A 295, 114 (2001).
15. H. S. Wio, S. Bouzat and B. von Haeften, in Proc. 21st IUPAP Int. Conf. on Statistical Physics,
STATPHYS21, A. Robledo and M. Barbosa (Eds.), Physica A 306C, 140-156 (2002).
16. R. Graham, in Instabilities and Nonequilibrium Structures, Eds. E. Tirapegui and D. Villaroel (D.
Reidel, Dordrecht, 1987); R. Graham and T. Tél, Phys. Rev. A 42, 4661 (1990); R. Graham and T.
Tél, in Instabilities and Nonequilibrium Structures III, E. Tirapegui and W. Zeller, eds. (Kluwert,
1991); O. Descalzi and R. Graham, Phys. Lett. A 170, 84 (1992); O. Descalzi and R. Graham, Z.
Phys. B 93, 509 (1994); H. S. Wio, in 4th Granada Seminar in Computational Physics, Eds. P.
Garrido and J. Marro (Springer-Verlag, Berlin, 1997), pg.135.
17. G. Izús et al., Phys. Rev. E 52, 129 (1995); G. Izús et al., Int. J. Mod. Phys. B 10, 1273 (1996).
18. G. Drazer and H. S. Wio, Physica A 240, 571 (1997).
19. D. H. Zanette, H. S. Wio and R. Deza, Phys. Rev. E 53, 353 (1996); F. Castelpoggi, H. S. Wio
and D. H. Zanette, Int. J. Mod. Phys. B 11, 1717 (1997); G. Izús, R. Deza and H. S. Wio, Phys.
Rev. E 58, 93 (1998); S. Bouzat and H. S. Wio, Phys. Lett. A 247, 297 (1998); G. G. Izús, R. R.
Deza and H. S. Wio, Comp. Phys. Comm. 121-122, 406 (1999).
20. B. von Haeften et al., Phys. Rev. E 69, 021107 (2004); B. von Haeften et al., in Noise in Complex
Systems and Stochastic Dynamics, Z. Gingl, J. M. Sancho, L. Schimansky-Geier, and J. Kerstez
(Eds.), Proc. SPIE 5471, 258-265 (2004).
21. H. S. Wio, in Noise and Fluctuations, Proc. 18th Int. Conf. on Noise and Fluctuations-ICNF2005,
Eds. T. Gonzalez, J. Mateos, D. Pardo, AIP 780, 55 (2005).
22. C. J. Tessone and H. S. Wio, Physica A 374, 46 (2006).
23. H. S. Wio, An Introduction to Stochastic Processes and Nonequilibrium Statistical Physics (World
Scientific, Singapore, 1994); A. S. Mikhailov, Foundations of Synergetics I (Springer-Verlag, Berlin,
1990).
24. G. Schmid, I. Goychuk, and P. Hänggi, Europhys. Lett. 56, 22 (2001); ibid., Phys. Biol. 1, 61
(2004).
25. P. Jung and J. W. Shuai, Europhys. Lett. 56, 29 (2001); J. W. Shuai and P. Jung, Phys. Rev. Lett.
88, 068102 (2003).
26. R. Toral, C. Mirasso, and J. Gunton, Europhys. Lett. 61, 162 (2003).
27. A. Pikovsky, A. Zaikin, and M. A. de la Casa, Phys. Rev. Lett. 88, 050601 (2002).
28. C. J. Tessone and R. Toral, Physica A 351, 106-116 (2005).
29. B. von Haeften, G. G. Izús, and H. S. Wio, Phys. Rev. E 72, 021101 (2005).
30. C. L. Schat and H. S. Wio, Physica A 180, 295 (1992).
31. P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
32. A. Zaikin et al., Phys. Rev. Lett. 90, 030601 (2003).
33. J. Garćıa-Ojalvo and J. M. Sancho, Noise in Spatially Extended Systems (Springer-Verlag, New
York, 1999).
34. It is worthwhile noting that when the parameter h is large enough, under some circumstances the
coupling term might become negative, in what is known as “inhibitory coupling” [35]. This is a
very interesting kind of coupling that has attracted much attention, both in neural and chemical
context, that we will not discuss here.
35. P. Dayan and L. F. Abbott, Theoretical neuroscience: Computational and mathematical modeling
of neural systems (MIT Press, Cambridge, 2001).
36. B. von Haeften and H. S. Wio, Physica A 376, 199 (2007).
37. P. Ao, cond-mat/0302081 (2003); P. Ao, J. Phys. A 37, L25-L30 (2004).
http://arxiv.org/abs/cond-mat/0302081
	Introduction
	System-size stochastic resonance
	Case of selective coupling
	Nonlocal Interaction
	Conclusions

  We analyze several aspects of the phenomenon of stochastic resonance in
reaction-diffusion systems, exploiting the nonequilibrium potential's
framework. The generalization of this formalism (sketched in the appendix) to
extended systems is first carried out in the context of a simplified scalar
model, for which stationary patterns can be found analytically. We first show
how system-size stochastic resonance arises naturally in this framework, and
then how the phenomenon of array-enhanced stochastic resonance can be further
enhanced by letting the diffusion coefficient depend on the field. A yet less
trivial generalization is exemplified by a stylized version of the
FitzHugh-Nagumo system, a paradigm of the activator-inhibitor class. After
discussing for this system the second aspect enumerated above, we derive from
it -through an adiabatic-like elimination of the inhibitor field- an effective
scalar model that includes a nonlocal contribution. Studying the role played by
the range of the nonlocal kernel and its effect on stochastic resonance, we
find an optimal range that maximizes the system's response.

Introduction
Stochastic resonance (SR) is nowadays a paradigm of the constructive effects of fluctuations on
nonlinear systems [1,2]. Sketchily, the phenomenon occurs whenever the Kramers’ rate for the
transition between attractors matches the typical frequency of a signal which is incapable by
itself to trigger that transition (i.e. it is subthreshold). Whereas several measures of SR can be
defined [the signal-to-noise ratio (SNR) and the spectral amplification factor (SAF) being the
main ones], theoretical analysis is usually carried on in terms of the two-state approximation
[1]. Since its discovery a quarter of century ago—and besides exploring related phenomena
like e.g., coherence resonance [3]—interest has gradually shifted towards increasingly complex
systems, networks and nonlinear media being the main directions. Instances of this trend are
the experiments carried out to explore the role of SR in sensory and other biological functions
[4], and experiments in chemical systems [5].
Our concern throughout this review will be with nonlinear media that can be described
as reaction–diffusion (RD) systems, namely those that can be thought of as a collection of
diffusively coupled nonlinear units. The possibility of enhancing the system’s response through
the coupling of those units [6,7,8,9,10,11,12] has been among the issues explored during the
last decade, together with the “naturalness” problem (how does nature manage to make the
a e-mail: wio@ifca.unican.es
2 Will be inserted by the editor
system’s response less dependent on a fine tuning of the noise intensity) or that of searching
for different ways to control the phenomenon [13,14].
In dissipative dynamical systems, the very notion of Lyapunov’s function is as useful as
that of attractor itself. Even when often it cannot be explicitly computed (because integrability
conditions are not readily met), it allows for picturesque reasoning in terms of “energy land-
scapes” or “attraction basins”. When those dynamical systems are submitted to forces that can
be modeled (à la Langevin) as stochastic, a new meaning—statistical in nature, akin to the
notion of “free energy”—is added to the picture (in fact, it is worth mentioning that the first
function known to have the Lyapunov property was Boltzmann’s H-function). Moreover, even
for vanishing noise intensity, the very existence of stochastic terms (a “transport matrix”) can
render the system well conditioned regarding integrability conditions. That was the rationale
behind the definition of nonequilibrium potential (NEP) [16,17,18,19], two approaches to which
are described in the Appendix. Such NEP is a special Lyapunov’s function of the associated
deterministic system, which for nonequilibrium systems plays a role similar to that played by
a thermodynamic potential in equilibrium thermodynamics [16]. It is closely related to the sta-
tionary solution of the system’s Fokker–Planck equation, and characterizes the global properties
of the dynamics: attractors, linear and relative stability of these attractors, height of the barri-
ers separating attraction basins. In addition, it allows to evaluate the transition rates among
the different attractors [16,17,18,19,20]. Regarding the problem of SR in extended systems, it
was shown that the knowledge of the NEP allows to obtain a rather complete picture of the
behavior of the output signal-to-noise ratio (SNR). The novelty with nonequilibrium extended
systems is that even pointlike attractors in the medium’s infinite-dimensional phase-space can
be nontrivial field configurations (real-space patterns).
In a series of recent papers we have studied the SR phenomenon for the transitions between
two different patterns [8,9,10,11,12,15], exploiting the concept of nonequilibrium potential. In
this review we discuss some recent results concerning different aspects of SR in RD systems. In
Sec. 2 we discuss the phenomenon of system-size stochastic resonance (SSSR), and show how
can it be analyzed and understood within a NEP framework [21]. In Sec. 3, after reviewing a
recent study on the enhancement of the SNR found for a scalar system with density-dependent
diffusivity [12], we discuss its extension [22] to an array of FitzHugh–Nagumo units [23]. In
Sec. 4—through an adiabatic-like elimination of the inhibitor field in an activator–inhibitor
system—an effective scalar system with a nonlocal term is derived, and the role of the local and
nonlocal interactions on the SR response studied. The main conclusions are finally summarized
in Sec. 5.
2 System-size stochastic resonance
Recent studies on biological models of the Hodgkin–Huxley type [24,25] have shown that ion
concentrations along cell membranes display intrinsic SR-like phenomena as the number of ion
channels is varied. A related result [26] shows that even in the absence of external forcing, the
regularity of the collective firing of a set of coupled excitable FitzHugh–Nagumo units is optimal
for a given number of elements. From a physics point of view, the same phenomenon—called
system-size stochastic resonance (SSSR)—has also been found in an Ising model as well as in
a set of globally coupled units described by a φ4 theory [27]. It has been even shown to arise in
opinion formation models [28].
Since the SSSR phenomenon is peculiar to extended systems, there is an obvious interest
in describing it within a NEP framework, that offers a very general framework for the study
of the dependence of SR and related phenomena on any of the system’s parameters. Here we
discuss in some detail a one-component (“scalar”) RD system [21], and briefly refer to other
cases analyzed in [27] and [29].
2.1 Review of the scalar model
For one-variable dynamical systems, the Lyapunov’s function can always be found by quadra-
ture. This property can be readily translated to scalar RD systems: the Lyapunov’s functional
Will be inserted by the editor 3
F [φ] fulfills the “potential” condition ∂tφ(y, t) = −δF [φ]/δφ(y, t), where δ/δφ(y, t) indicates a
functional derivative. This is also the NEP for a scalar transport matrix (i.e. a multiple of the
unit matrix in the medium’s infinite-dimensional phase space).
The specific model we shall focus on here has a piecewise linear reaction term, that mimics
general bistable RD models [23], e.g. those with a cubic-like reaction term. In the following,
we shall exploit some of the results on the influence of general (partially reflective or albedo)
boundary conditions found in [30], as well as previous studies of the NEP [17] and of SR
[8,9,11,15]. The particular dimensionless form of the deterministic model we start with is [30,8,9]
− φ+ φhΘ(φ− φc), (1)
where Θ(z) is the Heaviside step function, φc is the value at which the piecewise-linear “reaction
term” has the jump, and D is a phenomenological “diffusion coefficient”, not necessarily related
to γ in Sec. 2.2 [in Graham’s approach, Eq. (15), D/(∆y)2 would be the matrix elements Qνν±1
in a discretization of the Laplacian]. All the effects of the parameters that keep the system away
from equilibrium (such as the electric current in electrothermal devices like the ballast resistor
[23,30], or some external reactant concentration in chemical models) are now included in φc.
For the system to display a bistable behavior, it must be 0 < φc < φh.
We consider here the class of static structures φ(y) studied in [30]. They are even solutions
to the stationary (∂tφ = 0) version of Eq. (1) in the bounded domain [−yL, yL], with equal
albedo boundary conditions (b.c.) at both ends
∂φ(y, t)
y=±yL
= ∓k φ(±yL, t).
k > 0 is called the albedo parameter: the limit k → 0 yieds Neumann’s b.c. ∂yφ(y, t)|y=±yL = 0
and the k →∞ one, Dirichlet’s b.c. φ(±yL, t) = 0.
The explicit form of these static patterns is
φ(y) = φh

sinh(yc)ρ′
k, yL+y√
k, yL√
−yL ≤ y ≤ −yc,
1− cosh(y)ρ
yL−y±c√
k, yL√
−yc ≤ y ≤ yc,
sinh(yc)ρ′
k, yL−y√
k, yL√
yc ≤ y ≤ yL,
where ρ(k, ζ) = sinh(ζ) + k cosh(ζ) and ρ′(k, ζ) = ∂ρ/∂ζ. The coordinate values y±c at which
φ(yc) = φc are
y±c =
yL − ln
Z2 + 1− k2
1 + k
with Z =
. (3)
Each real solution y±c < yL to Eq. (3) represents a structure with a central “activated” zone
(φ > φc) and two lateral “resting” regions (φ < φc). Figure 5 in [30] displays the relation yc/yL
vs k, for several values of φc/φh.
Typical shapes of the arising patterns are shown in Fig. 1. Through a linear stability analysis
it has been shown [30] that the structure with the smallest “excited” region [that is with
yc = y+c , denoted by φu(y)] is unstable, whereas the other one [with yc = y
c , denoted by φ1(y)]
is linearly stable. The trivial homogeneous solution φ0(y) = 0 exists and is linearly stable for
any parameter set. These two linearly stable solutions (φ0 and φ1) are the only stable static
structures under albedo b.c. We will concentrate on the region of values of φc/φh, yL and k
where φ1 exists.
For the finite system with albedo b.c., the NEP is a functional of φ and a function of k, yL
and φc/φh. It has the expression [17]
F([φ], φc/φh, k, yL) =∫ yL
∫ φ(y,t)
[−φ′ + φhΘ(φ′ − φc)] dφ′ +
∂φ(y, t)
φ2(y, t)
4 Will be inserted by the editor
Fig. 1. Inhomogeneous static solutions to Eq. (1) for yL = D = φh = 1.0 and φc = 0.193 (highlighted
by the dashed horizontal line). Dash-dotted upper curve: φ1(y) for k = 1.0 (in this case, yc ≡ y−c > yL).
Lower curves: φ1(y) (solid line, for which +y
c is highlighted) and φu(y) (dotted line, for which −y+c is
highlighted) for k = 7.0.
Fig. 2. NEP F([φ], k, yL) evaluated at the inhomogeneous stationary solutions φ1(y) (lower branch)
and φu(y) (upper branch), as a function of: (a) system’s size yL, with k = 3.0; (b) albedo parameter k,
with yL = 1.2. The remaining parameters are D = 1.0, φc/φh = 0.193. The NEP for the homogeneous
stationary solution φ0(y) coincides with the horizontal axis.
When the NEP is evaluated at the inhomogeneous static solutions of Eq. (1) [Eqs. (2) and
(3)] it takes the explicit form [8,17]
Fu,1(φc/φh, k, yL) = F([φu,1], φc/φh, k, yL) (4)
= φ2h
−y±c (1− 2φcφh
+ sinh
y±c /
) ρ(k, (yL − y±c )/√D)
k, yL/
while at the trivial solution φ0 ≡ 0 it is F([φ0], yL) = F0 = 0.
Figure 2a depicts the nonequilibrium potential F([φ], yL) as a function of the system’s size
yL, keeping the albedo parameter k and the ratio φc/φh fixed. The curves correspond to the NEP
Fu,1(yL), whereas F0 coincides with the x–axis. Our focus is the “bistable zone” yL ' 0.72,
where φ1(y) exists. The unstable structure φu(y) is a saddle point for F [φ] (in the medium’s
infinite-dimensional phase space), so its NEP Fu(yL) > 0 (upper branch in Fig. 2a). On the
other hand, both φ1(y) (lower branch) and φ0 (x–axis) are local minima of the NEP.
One immediately notices that ∆F1(yL) ≡ Fu(yL)−F1(yL) is an (almost linearly) increasing
function of yL (this has a profound implication for SSSR, as we shall see below). Equation (3) has
real solutions only for yL ' 0.72. This corresponds to a supercritical saddle–node bifurcation,
at which both inhomogeneous structures pop up. Now, the most important feature in Fig. 2a is
that F1(yL) vanishes at a certain system’s size y∗L (≈ 1.0 for the given values of k and φc/φh).
At that point, the stable inhomogeneous structure φ1(y) and the trivial solution φ0(y) exchange
their relative stabilities.
Will be inserted by the editor 5
For completeness, in Fig. 2b we plot F([φ], k) for the same value of φc/φh as in Fig. 2a.
For the chosen value of yL it is always F1(k) < 0, and correspondingly there is no “stability
exchange” as a function of k. Also, the initially large ∆F1(k) decreases with k as Fu(k) and
F1(k) tend (for k →∞) to the values corresponding to Dirichlet’s b.c. [17].
2.2 Results for SSSR
By including an additive spatiotemporal noise source ξ(y, t) [15,33], Eq. (1) becomes a stochastic
partial differential equation for the random field φ(y, t). The simplest assumptions about ξ(y, t)
are that it is Gaussian, with zero mean and a correlation function given by 〈ξ(y, t)ξ(y′, t′)〉 =
2γ δ(t− t′)δ(y − y′), where γ denotes the noise strength.
As discussed in [8,9,11,15], known results for activation processes in multidimensional sys-
tems [31] allow us to estimate the activation rate using the following Kramers’-like expression
for the mean first-passage time for the transitions between attractors
〈τi〉 = τ0 exp
∆F i(yL)
where ∆F i(yL) = Fu(yL) − F i(yL), i = 0, 1. The prefactor τ0 is usually determined by the
curvature of F [φ] at its extrema. On one hand, it is typically several orders of magnitude smaller
than the average time 〈τ〉, while on the other it does not change significatively when varying the
system’s parameters around the “bistable point” y∗L, where F([φ0], y
L) = F([φ1], y
L). Hence we
may simplify the analysis by assuming here that τ0 is constant, and scale it out of our results.
The behavior of 〈τ〉 as a function of k and φc/φh has been shown in [8,9,17].
As done in [8], we now assume that the system is (adiabatically) subject to an external
harmonic variation of the parameter φc: φc(t) = φ0c + δφc cos(ωt) [9,15], and exploit the “two-
state approximation” [1] as in [9,11,15]. Such approximation reduces the whole dynamics on
the bistable potential landscape to one where the transitions occur only between the states
associated to the bottom of each well, hence the only dynamical contents resides in the transition
rates. Up to first order in the amplitude δφc (assumed to be small in order that the periodic
input be sub-threshold) the transition rates Wi adopt the form
µi ∓ αi
cos(ωt)
, (5)
where (at constant φh) µi ≈ exp[−∆F i(φ0c , yL)/γ] and αi ≈ µi (∂∆F i/∂φc|φ0c ), i = 0, 1. The
quantity inside parentheses can be obtained analytically using Eq. (4). These results allow to
calculate the autocorrelation function, the power spectrum density and finally the SNR, that
we indicate by R. The detailed calculation can be found in the appendix of [11]. Up to the
relevant order (the second) in the signal amplitude δφc, we obtain
4µ0µ1
(α0µ1 + α1µ0)2
µ0 + µ1
µ0 + µ1
Φ, (6)
where we have used the form of the αi to reduce the expression, and defined Φ = [2φh yc(yL)]2.
Figure 3 (left) is a plot of R as a function of the noise intensity γ for a fixed system’s length
yL, displaying the typical maximum that has become the fingerprint of the SR phenomenon.
In Fig. 3 (right), the roles of γ and yL are exchanged (R is plotted as a function of yL for fixed
γ). Such a response is the expected one for a system exhibiting SSSR. In both cases, the values
of k and φ0c/φh are kept fixed.
Within the NEP context and in this kind of systems, the phenomenon arises due to the
breakdown of the NEP’s symmetry. This means that (as shown in Fig. 2) when varying yL,
both attractors can exchange their relative stability. For yL = y∗L ≈ 1 both stable structures—
the inhomogeneous one φ1(y) and the trivial one φ0—have the same value for the NEP. For
yL < y
L, φ1(y) becomes less stable than φ0 so transitions from φ1(y) to φ0 are more frequent
6 Will be inserted by the editor
Fig. 3. Left: SNR vs the noise intensity γ, with system’s size yL = 1.1. Right: SNR vs yL, for γ = 0.1.
The remaining parameters are k = 3.0, D = 1.0, and φ0c/φh = 0.193.
(the barrier is lower) than in the reverse direction, thus reducing the system’s response. When
yL ∼ 0.72, φ1(y) and φu(y) coalesce and disappear, and the response is strictly zero (within
the linear response scheme implicit in the two-state approximation). When yL > y∗L, φ1(y)
becomes more stable than φ0, making now transitions from φ0 to φ1(y) more frequent than in
the reverse direction, and reducing again the system’s response. Clearly, the system’s response
has a maximum when both attractors have the same stability (yL = y∗L), and decays when
departing from that situation. Hence, for this system and within this framework, SSSR arises
as a particular case of the more general discussion done in [11]. It should not come as a surprise
to find an analogy with the mechanism of double stochastic coherence described in [32], where
the NEP’s symmetry is induced by (an additional, multiplicative) noise.
By comparing figures 2a and 3 it becomes apparent that the value of yL at which the SNR
has its maximum differs slightly from y∗L (where the crossing between F
1 and F0 takes place).
The origin of this discrepancy is the following: whereas on qualitative grounds we have argued
that the maximum of the SNR should be related to the potential being symmetric (both wells
having the same “energy”) [11], the exact condition is that the transition rates between both
wells be equal. In general, due to small differences between the curvatures at the bottom of each
well, those rates become equal for values of yL slightly different from the one at the symmetric
case. Although by adopting here a constant value of τ0 we have assumed equal curvatures, there
is still a difference in the values of the αi, since the ∂∆F i/∂φc|φ0c , i = 0, 1 are slightly different
(a fact reflected in the dependence of Φ on yL).
Additional light can be shed on the phenomenon when viewed from a different angle. Figure
4 is a plot of R as a function of k at fixed values of γ, yL and φ0c/φh. It exhibits a broad resonance
since for k not too large (indicating a high reflectiveness at the boundary or a reduced exchange
with the environment) R increases with k, whereas it slightly decreases for larger k values (where
the system’s boundaries become absorbent). An explanation of this behavior in terms of the
NEP has been given in [29]: as already observed, the NEP’s symmetry is broken for this value of
yL; moreover, whereas the lower branch in Fig. 2b goes rapidly towards the value corresponding
to Dirichlet’s b.c. the upper branch keeps increasing, thus degrading the SNR. In any case, the
fact that the resonance is broad indicates the robustness of the system’s response with regard
to k, a parameter that (together with γ) encodes the coupling with the environment.
We stress the fact that the NEP framework put forward in this review allows to study
SSSR between whole patterns. The explanation offered in [27] to the phenomenon resorted to a
collective variable X ≈ (1/N)
j=1 xj , and to the fact that the noise in the effective stochastic
differential equation for X scaled with size. In [29] it was shown that all the cases discussed in
[27] can be put within the same NEP framework than the above studied scalar model. In fact,
the aforementioned almost linear increasing dependence of ∆F1(yL) on yL can be interpreted
as a noise scaling with size. There are however situations where the NEP’s symmetry is retained
as the system’s size is varied. We may then speak of a genuinely noise-scaling SSSR, in contrast
to the cases that could be called NEP symmetry breaking SSSR [29].
Will be inserted by the editor 7
Fig. 4. SNR vs k for yL = 1.2, γ = 0.1 and φ
c/φh = 0.193.
3 Case of selective coupling
In this section we analyze SR in two extended systems with density-dependent diffusive-like cou-
pling: an extension of the scalar RD model considered in Sec. 2 [12], and an array of FitzHugh–
Nagumo [23] units.
3.1 Scalar model
Here we extend the one-component RD model discussed in Sec. 2 by letting the diffusive pa-
rameter D in Eq. (1) depend on the field φ(x, t). As a matter of fact, since in the ballast resistor
[23,30] the thermal conductivity is a function of the energy density, the resulting equation for
the temperature field includes a temperature-dependent diffusion coefficient in a natural way.
The form of the governing equation is now
∂tφ(x, t) = ∂x [D(φ)∂xφ] + f(φ) + ξ(x, t), (7)
with ξ(x, t) and f(φ) as in Sec. 2.
As it was done for the reaction term, a simple choice (that retains however the qualitative
features of the system) is to consider the following dependence of the diffusion term on the field
variable
D(φ) = D0[1 + hΘ(φ− φc)].
For simplicity, here we choose the same threshold φc for the reaction term and the diffusion
coefficient.
We assume the system to lie in a bounded domain [−L,L], with Dirichlet b.c. at both ends:
φ(±L, t) = 0. The form of the patterns is analogous to what has been obtained in Sec. 2, the
only difference being that in the present case dφ/dx|xc is discontinuous and the area of the
“activated” central zone depends on h.
As before, the indicated patterns are extrema of the NEP: the unstable pattern φu(x) is
a saddle-point of this functional, separating the attractors φ0(x) and φs(x). For the case of a
field-dependent diffusion coefficient D(φ(x, t)) as described by Eq. (7), the NEP reads [12]
F [φ] =
D(φ′)f(φ′) dφ′ +
Given that ∂tφ = −[1/D(φ)](δF/δφ), one finds dF/dt = −
(δF/δφ)2 dx ≤ 0, thus warranting
the Lyapunov’s functional property.
Whereas in Sec. 2 we kept φ0c constant and varied yL, we now vary instead φ
c at constant
L. Similarly as before, both linearly stable states have the same value of the NEP (i.e., they
are equally stable) at some value φ∗c of the threshold. The way F [φ] depends on φ0c resembles
the dependence on yL shown in Sec. 2, but now ∆Fs is an (almost linearly) decreasing function
8 Will be inserted by the editor
0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14
γ (arb. units)
-0.2 -0.1 0.0 0.1 0.2
h (arb.units)
 (arb. units)
0.90 0.95 1.00 1.05 1.10
Fig. 5. Left: SNR vs the noise intensity γ, for D0 = 1.0 and h = 0.0 (full line), −0.25 (dashed line)
and 0.25 (dotted line). Right: the maximum Rmax of the SNR curve as a function of the selectiveness
h of the coupling, for D0 = 0.9 (dashed line), 1.0 (full line) and 1.1 (dotted line). The arrows a and b
indicate the response gain due respectively to a homogeneous increase of the coupling and to a selective
one. The larger gain in the second case is apparent. The inset shows the dependence of Rmax on D0 for
h = −0.25 (lower line), 0.0 and 0.25 (upper line). The remaining parameters are L = 1.0, δφc = 0.01
and Ω = 0.01.
of φ0c and both inhomogeneous structures coalesce and disappear through a subcritical saddle–
node bifurcation. As in the previous case, we analyze only the neighborhood of φ∗c . We shall
moreover consider only the neighborhood of h = 0, where the main trends of the effect can be
captured.
Figure 5 (left) depicts the dependence of R on the noise intensity γ for three values of h,
each curve displaying the typical SR maximum. Figure 5 (right) is a plot of the value Rmax
of these maxima as a function of h. The dramatic increase of Rmax (several dB for a small
positive variation of h) is apparent, and shows the strong effect that the selective coupling (or
field-dependent diffusivity) has on the system’s response.
It must be noted that the only two approximations made in order to derive our results—
namely the Kramers-like expression and the two-level approximation used for the evaluation of
the correlation function [12]—break down for large positive values of h because for increasing
selectivity the curves of F [φ] vs φ0c shift towards the left [12], which in turn means that the
barrier separating the attractors at φ∗c tends to zero. This effect is basically the same as the one
discussed in Refs. [9,11] in connection with global diffusivity D0. It is also worth noting that
except for the two aforementioned approximations, all the previous results (e.g. the profiles
of the stationary patterns and the corresponding values of the nonequilibrium potential) are
analytically exact.
3.2 FitzHugh–Nagumo model
Here we study an array of FitzHugh–Nagumo [23] units, with a density-dependent (diffusive-
like) coupling. The NEP for this system was found within the excitable regime and for particular
values of the coupling strength [19]. In the general case, however, the form of the NEP has not
been found yet. Hence, we have resorted to a study based on numerical simulations, analyzing
the influence of different parameters on the system’s response. Nevertheless, the idea of the
existence of such a NEP has always underlied this study. The results show that the enhancement
of the SNR found for the scalar system [12] is robust, and that the indicated non-homogeneous
coupling could clearly contribute to enhance the SR phenomenon in more general situations.
We consider a simplified version of the FitzHugh–Nagumo model [12,17,23], which has been
useful for gaining qualitative insight into the excitable and oscillatory dynamics in neural and
chemical systems. It consist of two variables:
– a (fast) activator field u, that in the case of neural systems represents the voltage variable,
while in chemical systems represents the concentration of an autocatalytic species.
Will be inserted by the editor 9
– an inhibitor field v, associated (within a neural context) to the concentration of potassium
ions in the medium, and that in a general chemical reaction inhibits the generation of the
u species.
Instead of considering the usual cubic-like nonlinear form, we use a piecewise linear version
∂u(x, t)
Du(u)
+ f(u)− v + ξ(x, t) (8)
∂v(x, t)
Dv(v)
+ β u− α v, (9)
where f(u) = −u + Θ(u − φc), and ξ(x, t) is a δ–correlated white Gaussian noise, as before.
γ indicates the noise intensity and φc is the “discontinuity” point, at which the piecewise
linearized function f(u) presents a jump. The parameter � = τu/τv � 1 indicates the timescale
ratio between the (fast) activator and the inhibitor. Together with α and β, it is chosen to
correspond to the excitable regime. We consider Dirichlet b.c. at x = ±L. Although the results
are qualitatively the same as those that could appear considering the usual FitzHugh–Nagumo
equations, this simplified version allows us to compare directly with the previous analytical
results for this system [12,15].
As in [12], we assume that the diffusion coefficient Du(u) is not constant, but depends on
the field u according to Du(u) = D0u [1 + hΘ(u− φc)]. This form implies that the value of
Du(u) depends “selectively” on whether the field u fulfills u > φc or u < φc. D0u is the value of
the diffusion constant without such “selective” term, and h indicates the size of the difference
between the diffusion constants in both regions [if h = 0 then Du(u) = D0u constant]. Dv(v) is
the diffusion for the inhibitor v, that we assume to be homogeneously constant.
This system is known to exhibit two stable stationary patterns. One of them is u(x) = 0,
v(x) = 0, while the other is one with nonzero values. Further, we consider, as before, that
an external, periodic, signal enters into the system through the value of the threshold φc,
φc(t) = φ0c + δφ cos(ωt), where ω is the signal frequency, and δφ its intensity.
All the results were obtained through numerical simulations of the system. The continuous
version of the system indicated by Eqs. (8), (9) was transformed into a second-order spatially
discrete one
u̇i = Dui∆ui + (Dui+1 −Dui−1)(ui+1 + ui−1) + f(ui)− vi + ξi(t),
v̇i = Dv∆vi + β ui − α vi,
with ∆φi ≡ (φi−1 + φi+1 − 2φi). The extensive numerical simulations performed for a set of
equations were done exploiting Heun’s algorithm [33].
In this spatially-extended system there are different ways of measuring the overall system’s
response to the external signal. In particular, we show the evaluated output SNR in two different
ways (the units being given in dB):
– SNR for the middle element of the chain evaluated over the dynamical evolution of uN/2,
that we call SNR2 (however, having Dirichlet b.c., the local response depends on the distance
to the boundaries).
– In order to measure the overall response of the system to the external signal, we computed
the SNR as follows: We digitized the system’s dynamics to a dichotomic process s(t): At
time t the system has an associated value of s(t) = 1 (0) if the Hilbert distance to pattern
1 (0) is lower than to the other pattern. Stated in mathematical terms, we computed the
distance D2[·, ·] defined by
D2[f, g] =
dx [f(x)− g(x)]2
in L2([−L,L]), the Hilbert space of the real-valued functions in that interval. At time t, a
digitized process is computed by means of
s(t) =
1 if D2 [Pu1 (x), u(x, t)] < D2 [Pu0 (x), u(x, t)]
0 if D2 [Pu1 (x), u(x, t)] ≥ D2 [Pu0 (x), u(x, t)]
10 Will be inserted by the editor
Model Numerical
α β � φ0c δφ ω Dv ∆t N
0.3 0.4 0.03 0.52 0.4 5π/8 1.0 10−3 51
Table 1. Fixed parameters for the FitzHugh–Nagumo model with Du(u)
Fig. 6. SNR2 and SNRp vs the noise intensity γ, for D
u = 0.3 and h = −2.0 (+), −1.0 (4), 0.0 (�),
1.0 (�), and 2.0 (©). Both measures reveal a systematic enhancement of the SNR as h increases.
-2 0 2 4 6 8 10
-2 0 2 4 6 8 10
Fig. 7. SNR2 and SNRp vs the selectiveness of coupling h, for D
u = 0.3 and γ = 0.032 (©), 0.32 (�),
0.6 (�), 1.2 (4), and 3.2 (×).
We call this “global-like” measure SNRp.
The parameters kept fixed have been summarized in Table 1. The simulation was repeated 250
times for each parameter set, and the SNR was computed by recourse of the average power
spectral density.
Figure 6 depicts the results for the different SNR measures we have previously defined, as
functions of the noise intensity γ. For both measures it is apparent that there is an enhancement
of the response for h > 0, when compared with the h = 0 case, while for h < 0 the response is
smaller.
In Fig. 7 we show again the response’s measures, but now as functions of h. We have plotted
the maximum of each SNR curve for D0u = 0.3, and γ = 0.01, 0.1, and 0.3. It is clear that there
exists an optimal value of γ for which the response is largest. The rapid fall in the response for
h < 0 is also apparent.
In Fig. 8 we show the dependance of SNR on h, for different values of the diffusion which
depends on the activator density D0u. It is apparent that the response becomes larger when
the value of D0u is larger. However, as was discussed in [8,15], it is clear that for still larger
values of D0u, the symmetry of the underlying potential (that is the relative stability between
the attractors) is broken and the response finally falls down.
The previous figures clearly show that the response to the external signal grows with the
“selectiveness” of the coupling, showing the robustness of the phenomenon presented in [12,15].
Will be inserted by the editor 11
-2 0 2 4 6 8 10
-2 0 2 4 6 8 10
Fig. 8. SNR2 and SNRp vs h for γ = 0.032 and D
u = 0.0 (©), 0.1 (�), and 0.3 (�).
4 Nonlocal Interaction
Let us consider again a system like the one described by Eqs. (8), (9), but now assume that
Du and Dv are constant. In Ref. [18] it was assumed that the inhibitor-like field has a diffusive
transport behavior, and is fast enough that can be adiabatically eliminated, thus yielding an
effective scalar RD equation with a nonlocal term, characterized by a diffusive kernel G(x, x′).
After briefly reviewing the derivation of the NEP for that situation, we shall assume in this
section that the transport mechanism of the adiabatically eliminated inhibitor-like field is of
nondiffusive character, thus yielding a kernel H(x, x′) that is more localized in space, and with
a controllable interaction range.
Following Ref. [18], let the system be defined by
∂tu(x, t) = Du∂
xu(x, t)− u(x, t) +Θ[u(x, t)− a]− v(x, t),
�−1 ∂tv(x, t) = Dv∂
xv(x, t) + βu(x, t)− αv(x, t), (10)
where � was defined after Eqs. (8), (9), and let it be confined to the domain [−L,L], with
Dirichlet b.c. at both ends: u(±L, t) = v(±L, t) = 0. Contrarily to the standard hypothesis,
we now assume that the inhibitor is much faster than the activator (i.e. τv � τu). In the limit
�→∞, we can rewrite Eq. (10) as
∂tu(x, t) = Du∂
xu(x, t)− u(x, t) +Θ[u(x, t)− a]− v(x, t),
0 = Dv∂
xv(x, t) + βu(x, t)− αv(x, t).
In the last pair of equations we can eliminate the inhibitor (which is now slaved to the activator)
by solving the second equation using the Green’s function method
[−Dv∂2x + α]G(x, x
′) = δ(x− x′),
v(x) = β
dx′G(x, x′)u(x′),
where the Green’s function G(x, x′) is given by
G(x, x′) =
[sinh k(L− x′)/ sinh 2kL] sinh k(L+ x) x < x′,
[sinh k(L+ x′)/ sinh 2kL] sinh k(L− x) x > x′,
with k = (α/Dv)1/2. This slaving procedure reduces our system to a nonlocal equation for the
activator only, that has the form
∂u(x, t)
∂2u(x, t)
+ f(u)− β
G(x, x′)u(x′) dx′. (11)
12 Will be inserted by the editor
From this equation, and taking into account the symmetry of the Green’s function G(x, x′), we
can obtain the Lyapunov functional for this system, which has the form
F [u] =
f(w) dw +
dx′G(x, x′)u(x′)u(x)
. (12)
This spatial nonlocal term in the NEP takes into account the repulsion between activated
zones. When two activated zones come near each other, the exponential tails of the inhibitor
concentration overlap, increasing its concentration between both activated zones and creating
an effective repulsion between them. Hence the Green’s function plays the role of an exponential
screening between the activated zones. In Ref. [9] the knowledge of such NEP was exploited to
study SR on the system indicated by Eq. (11).
The starting point of our present analysis will be the effective, nonlocal and stochastic RD
equation for the real (activator-like) field φ(x, t), analogous to Eq. (11), defined in the one
dimensional domain x ∈ [−L,L] by
+ f(φ)− β
H(x, x′)φ(x′) dx′ + ξ(x, t),
where the diffusivityD is constant, and we assume a cubic nonlinear term f(φ) = φ (φ−b) (2−φ).
Here ξ(x, t) is an additive Gaussian white noise, as in the previous cases. As before, the system
is subject to Dirichlet b.c. φ(±L, t) = 0.
Similarly to Eq. (11), this system could be written in a variational form, with the func-
tional F [φ] given by Eq. (12). As anticipated, here we consider a nondiffusive kernel, with a
controllable interaction range. In order to keep our analysis simple we propose the following
H(x, x′) =
1/2; |x− x′| ≤ l
0; |x− x′| > l, (13)
which allows the analysis by just varying the interaction range 2l.
The new effective RD equation contains local and nonlocal couplings (corresponding to the
diffusive and the nonlocal contribution, respectively). The last one contains the nonlocal kernel
given by Eq. (13), with a variable range 2l, that corresponds to the interaction of the field at
points x′ ∈ [x − l, x + l]. However, such points will contribute if and only if they are inside
the domain [−L,L]. We are now in position to study the role played by the nonlocal kernel
(particularly by its range 2l) on the SR phenomenon.
As before, the SR between stationary solutions was investigated in terms of the two-state
approach (all the details about the procedure and the evaluation of the SNR can be found
in Refs. [11,15]). As usual, we subject our system to a weak external signal b = b0 + A(t) =
b0 +∆b cos(ωst), rocking the NEP [15]. In order to have a subthreshold signal, the amplitude
∆b should satisfy ∆b� b. We have chosen b0 as the value of b at which F [φs] = F [φ = 0] = 0
when β = 0.
Up to first order in the small amplitude ∆b, the transition rates Wi and the functions αi have
the form indicated in Eqs. (5), (6). But now Φ, that depends on the inhomogeneous attractor
φs, has the form
After fixing the length of the system L and the kernel range 2l we can use the above indicated
expressions to find the SNR, that shows the usual bell-shaped form of stochastic resonance as
a function of the noise intensity.
In Fig. 9 we show the dependence of the SNR on the kernel range 2l for a fixed value of 2L.
There is a nonmonotonic behavior in the system’s response against variation of 2l, that can be
explained by the following facts:
– On one hand, the transition rates are decreasing functions of the range 2l for fixed 2L.
Therefore, the ratio µ1µ2/(µ1 +µ2) in the expression for the SNR also reflects this behavior.
Will be inserted by the editor 13
– The other factor in this expression has a maximum for a kernel range that corresponds to
“first neighbor” sites, namely x = x′ ± l.
Fig. 9. SNR as a function of 2l for 2L fixed. The parameters are D = 0.6, 2L = 6.25, β = 0.02 and
b0 = 0.719123.
The maximum in the system’s response as a function of the kernel range is due to the interplay
between these two factors. From this analysis of the comparative weights of the local (diffusive)
and nonlocal terms contributing to the SR response, it is apparent that the range of such
nonlocal kernel has an optimum value yielding a maximum for the SNR [36].
5 Conclusions
We have discussed three different aspects of the phenomenon of stochastic resonance in reaction–
diffusion systems, within the nonequilibrium potential’s framework. In first place we have dis-
cussed system-size SR in a scalar model. Even though we have not shown the details here, it
has been also possible to also study other cases [29]. In particular, a model of globally coupled
nonlinear oscillators discussed in [27], showing that it can also be described within the NEP
framework, with SSSR arising through an “effective” scaling of the noise intensity with the
system’s size.
In second place we presented a study of SR in systems with a density-dependent (diffusive-
like) coupling. We initially discuss the case of a scalar system [12], and afterwards extent the
analysis to an array of FitzHugh–Nagumo units, with a field-dependent activator diffusion [22].
For the second system, when both diffusions are constant (that is: Du > 0 and Dv = 0), has
a known form of the NEP [15]. However, in the general case we have not been able to find
the form of the NEP (but the idea of such a NEP is always underlying our analysis) and have
to resort to an analysis based on numerical simulations. The result shows that the system’s
response is enhanced due to the particular form of the non-homogeneous coupling. From such
results, we can conclude that the phenomenon of enhancement of the SNR, due to a selectivity
in the coupling, initially found for a scalar system [12] is robust, and that the indicated non-
homogeneous coupling could clearly contribute to enhance the SR phenomenon in very general
systems.
Finally, we analyzed an activator-like field including a nonlocal contribution that arise
through an effective adiabatic elimination of an auxiliary (inhibitor-like) field. By exploiting the
knowledge of the nonequilibrium potential in such a case, we have analyzed the dependence of
the SNR on the nonlocal interaction kernel range, founding that there is an optimal value of the
kernel range, yielding a maximum in the system’s response, corresponding to a very localized
interaction.
The indicated results clearly show that the “nonequilibrium potential” (even if not known
in detail [10]) offers a very useful framework to analyze a wide spectrum of characteristics
associated to SR in spatially extended or coupled systems. For instance, within this framework,
14 Will be inserted by the editor
the phenomenon of SSSR looks—as other aspects of SR in extended systems [11]—as a natural
consequence of a breaking of the symmetry of the NEP [27].
Appendix: Brief review of the nonequilibrium potential scheme
Loosely speaking, the notion of NEP is an extension to nonequilibrium situations of that of
equilibrium thermodynamic potential. In order to introduce it, we consider a general system of
nonlinear stochastic equations (admitting the possibility of multiplicative noises)
q̇ν = Kν(q) + gνi (q) ξi(t), ν = 1, . . . , n; (14)
where repeated indices are summed over. Equation (14) is stated in the sense of Itô. The {ξi(t)},
i = 1, . . . ,m ≤ n are mutually independent sources of Gaussian white noise with typical strength
Graham’s approach
The Fokker–Planck equation corresponding to Eq. (14) takes the form
Kν(q)P +
∂qν ∂qµ
Qνµ(q)P (15)
where P (q, t; γ) is the probability density of observing q = (q1, . . . , qn) at time t for noise
intensity γ, and Qνµ(q) = gνi (q) g
i (q) is the matrix of transport coefficients of the system,
which is symmetric and non-negative. In the long time limit (t→∞), the solution of Eq. (15)
tends to the stationary distribution Pst(q). According to [16], the NEP Φ(q) associated to Eq.
(15) is defined by
Φ(q) = − lim
γ lnPst(q, γ). (16)
In other words
Pst(q) d
nq = Z(q) exp
+O(γ)
where Φ(q) is the NEP of the system and Z(q) is defined as the limit
lnZ(q) = lim
lnPst(q, γ) +
Here dΩq = dnq/
G(q) is the invariant volume element in the q-space and G(q) is the deter-
minant of the contravariant metric tensor (for the Euclidean metric it is G = 1). It was shown
[16] that Φ(q) is the solution of a Hamilton–Jacobi-like equation (HJE)
Kν(q)
Qνµ(q)
and Z(q) is the solution of a linear first-order partial differential equation depending on Φ(q)
(not shown here).
Equation (16) and the normalization condition ensure that Φ is bounded from below. Fur-
thermore, from the separation of the streaming velocity of the probability flow in the steady
state into conservative and dissipative parts, it follows that
dΦ(q)
= Kν(q)
∂Φ(q)
Qνµ(q)
i.e. Φ is a LF for the dynamics of the system when fluctuations are neglected. Under the
deterministic dynamics, q̇ν = Kν(q), Φ decreases monotonically and takes a minimum value on
attractors. In particular, Φ must be constant on all extended attractors (such as limit cycles or
strange attractors) [16].
Will be inserted by the editor 15
Ao’s approach
An alternative way to look into this problem is due to Ao [37]. Let us refer again to the system
in Eq. (14). Following [37], we introduce now the auxiliary matrix
H−1(q) = S(q) + A(q),
with S(q) a symmetric matrix while A(a) is an antisymmetric one. H−1(q) is now used to
rewrite the initial system as
H−1(q) q̇ = H−1(q) K(q, t) + H−1(q) ξ(q, t) = −∇Φ(q, t) + η(q, t),
where −∇Φ(q, t) = H−1(q) K(q, t), η(q, t) = H−1(q) ξ(q, t) and H−1(q) q̇ = −∇Φ(q, t) +
η(q, t). The new stochastic variables, η(q, t), fulfill
〈η(q, t)ηT (q, t′)〉 = 2 S(q) δ(t− t′) = 2 G−1(q)Q(q)[G−1(q)]T δ(t− t′),
that imposes a condition on the arbitrary definition of S(q) as we have
[S(q) + A(q)] Q(q) [S(q)−A(q)] = S(q).
As shown by Ao, the last equation also implies H(q) + HT (q) = 2 Q(q). We also have
H−1(q) K(q, t)
= ∇× [−∇Φ(q, t)] = 0. From the previous equations, we have in principle
all the needed conditions to determine H−1(q), and to obtain from it the potential Φ(q, t). The
interesting feature of this approach is that it resorts neither to Pst(q) nor to the small-noise
limit, thus being applicable in principle to more general situations.
Acknowledgements
The authors acknowledge the collaboration of B. von Haeften, S. Bouzat, C. J. Tessone, G. G.
Izús, M. Kuperman, S. Mangioni, A. Sánchez, F. Castelpoggi, in different aspects and/or stages
of this work. HSW thanks the European Commission for the award of a Marie Curie Chair at
the Universidad de Cantabria, Spain.
References
1. L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni, Rev. Mod. Phys. 70, 223 (1998).
2. T. Wellens, V. Shatokin, and A. Buchleitner, Rep. Prog. Phys. 67, 45 (2004).
3. A. S. Pikovsky and J. Kurths, Phys. Rev. Lett. 78, 775 (1997).
4. J. K. Douglas et al., Nature 365, 337 (1993); J. J. Collins et al., Nature 376, 236 (1995); S. M.
Bezrukov and I. Vodyanoy, Nature 378, 362 (1995).
5. A. Guderian, G. Dechert, K. Zeyer, and F. Schneider, J. Phys. Chem. 100, 4437 (1996); A. Förster,
M. Merget, and F. Schneider, J. Phys. Chem. 100, 4442 (1996); W. Hohmann, J. Müller, and F.
W. Schneider; J. Phys. Chem. 100, 5388 (1996).
6. J. F. Lindner et al., Phys. Rev. E 53, 2081 (1996).
7. A. R. Bulsara and G. Schmera, Phys. Rev. E 47, 3734 (1993); P. Jung, U. Behn, E. Pantazelou,
and F. Moss, Phys. Rev. A 46, R1709 (1992); P. Jung and G. Mayer-Kress, Phys. Rev. Lett. 74,
2130 (1995); J. F. Lindner et al., Phys. Rev. Lett. 75, 3 (1995); F. Marchesoni, L. Gammaitoni,
and A. R. Bulsara, Phys. Rev. Lett. 76, 2609 (1996).
8. H. S. Wio, Phys. Rev. E 54, R3075 (1996); H. S. Wio and F. Castelpoggi, Unsolved Problems of
Noise, Proc. Conf. UPoN’96, C. R. Doering, L. B. Kiss, and M. Schlesinger, Eds. (World Scientific,
Singapore, 1997), pg. 229; F. Castelpoggi and H. S. Wio, Europhys. Lett. 38, 91 (1997).
9. F. Castelpoggi and H. S. Wio, Phys. Rev. E 57, 5112 (1998).
10. H. S. Wio et al., Physica A 257, 275 (1998); M. Kuperman, H. S. Wio, G. Izús, and R. Deza, Phys.
Rev. E 57, 5122 (1998).
11. S. Bouzat and H. S. Wio, Phys. Rev. E 59, 5142 (1999).
16 Will be inserted by the editor
12. B. von Haeften, R. Deza, and H. S. Wio, Phys. Rev. Lett. 84, 404 (2000).
13. C. J. Tessone, H. S. Wio, and P. Hänggi, Phys. Rev. E 62, 4623 (2000).
14. M. A. Fuentes, R. Toral, and H. S. Wio, Physica A 295, 114 (2001).
15. H. S. Wio, S. Bouzat and B. von Haeften, in Proc. 21st IUPAP Int. Conf. on Statistical Physics,
STATPHYS21, A. Robledo and M. Barbosa (Eds.), Physica A 306C, 140-156 (2002).
16. R. Graham, in Instabilities and Nonequilibrium Structures, Eds. E. Tirapegui and D. Villaroel (D.
Reidel, Dordrecht, 1987); R. Graham and T. Tél, Phys. Rev. A 42, 4661 (1990); R. Graham and T.
Tél, in Instabilities and Nonequilibrium Structures III, E. Tirapegui and W. Zeller, eds. (Kluwert,
1991); O. Descalzi and R. Graham, Phys. Lett. A 170, 84 (1992); O. Descalzi and R. Graham, Z.
Phys. B 93, 509 (1994); H. S. Wio, in 4th Granada Seminar in Computational Physics, Eds. P.
Garrido and J. Marro (Springer-Verlag, Berlin, 1997), pg.135.
17. G. Izús et al., Phys. Rev. E 52, 129 (1995); G. Izús et al., Int. J. Mod. Phys. B 10, 1273 (1996).
18. G. Drazer and H. S. Wio, Physica A 240, 571 (1997).
19. D. H. Zanette, H. S. Wio and R. Deza, Phys. Rev. E 53, 353 (1996); F. Castelpoggi, H. S. Wio
and D. H. Zanette, Int. J. Mod. Phys. B 11, 1717 (1997); G. Izús, R. Deza and H. S. Wio, Phys.
Rev. E 58, 93 (1998); S. Bouzat and H. S. Wio, Phys. Lett. A 247, 297 (1998); G. G. Izús, R. R.
Deza and H. S. Wio, Comp. Phys. Comm. 121-122, 406 (1999).
20. B. von Haeften et al., Phys. Rev. E 69, 021107 (2004); B. von Haeften et al., in Noise in Complex
Systems and Stochastic Dynamics, Z. Gingl, J. M. Sancho, L. Schimansky-Geier, and J. Kerstez
(Eds.), Proc. SPIE 5471, 258-265 (2004).
21. H. S. Wio, in Noise and Fluctuations, Proc. 18th Int. Conf. on Noise and Fluctuations-ICNF2005,
Eds. T. Gonzalez, J. Mateos, D. Pardo, AIP 780, 55 (2005).
22. C. J. Tessone and H. S. Wio, Physica A 374, 46 (2006).
23. H. S. Wio, An Introduction to Stochastic Processes and Nonequilibrium Statistical Physics (World
Scientific, Singapore, 1994); A. S. Mikhailov, Foundations of Synergetics I (Springer-Verlag, Berlin,
1990).
24. G. Schmid, I. Goychuk, and P. Hänggi, Europhys. Lett. 56, 22 (2001); ibid., Phys. Biol. 1, 61
(2004).
25. P. Jung and J. W. Shuai, Europhys. Lett. 56, 29 (2001); J. W. Shuai and P. Jung, Phys. Rev. Lett.
88, 068102 (2003).
26. R. Toral, C. Mirasso, and J. Gunton, Europhys. Lett. 61, 162 (2003).
27. A. Pikovsky, A. Zaikin, and M. A. de la Casa, Phys. Rev. Lett. 88, 050601 (2002).
28. C. J. Tessone and R. Toral, Physica A 351, 106-116 (2005).
29. B. von Haeften, G. G. Izús, and H. S. Wio, Phys. Rev. E 72, 021101 (2005).
30. C. L. Schat and H. S. Wio, Physica A 180, 295 (1992).
31. P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
32. A. Zaikin et al., Phys. Rev. Lett. 90, 030601 (2003).
33. J. Garćıa-Ojalvo and J. M. Sancho, Noise in Spatially Extended Systems (Springer-Verlag, New
York, 1999).
34. It is worthwhile noting that when the parameter h is large enough, under some circumstances the
coupling term might become negative, in what is known as “inhibitory coupling” [35]. This is a
very interesting kind of coupling that has attracted much attention, both in neural and chemical
context, that we will not discuss here.
35. P. Dayan and L. F. Abbott, Theoretical neuroscience: Computational and mathematical modeling
of neural systems (MIT Press, Cambridge, 2001).
36. B. von Haeften and H. S. Wio, Physica A 376, 199 (2007).
37. P. Ao, cond-mat/0302081 (2003); P. Ao, J. Phys. A 37, L25-L30 (2004).
http://arxiv.org/abs/cond-mat/0302081
	Introduction
	System-size stochastic resonance
	Case of selective coupling
	Nonlocal Interaction
	Conclusions
