Table Of ContentCorrelation transfer in stochastically driven oscillators over long and short time scales
Aushra Abouzeid and Bard Ermentrout
University of Pittsburgh
(Dated: January 11, 2011)
In the absence of synaptic coupling, two or more neural oscillators may become synchronized by
virtue of the statistical correlations in their noisy input streams. Recent work has shown that the
degree of correlation transfer from input currents to output spikes depends not only on intrinsic
oscillator dynamics, but also depends on the length of the observation window over which the
correlationiscalculated. Inthispaperweusestochasticphasereductionandregularperturbations
toderivethecorrelationofthetotalphaseelapsedoverlongtimescales,aquantitywhichprovides
aconvenientproxyforthespikecountcorrelation. Overshorttimescales,wederivethespikecount
correlationdirectlyusingstraightforwardprobabilisticreasoningappliedtothedensityofthephase
difference. Our approximations show that output correlation scales with the autocorrelation of the
1 phase resetting curve over long time scales. We also find a concise expression for the influence of
1 the shape of the phase resetting curve on the initial slope of the output correlation over short time
0
scales. These analytic results together with numerical simulations provide new intuitions for the
2
recent counterintuitive finding that type I oscillators transfer correlations more faithfully than do
n type II over long time scales, while the reverse holds true for the better understood case of short
a time scales.
J
0
While the jury is still out on the functional role of cording to the bifurcations that occur as the dynamical
1
synchrony and correlations in neural firing, the ubiquity system goes from a stable rest state to a stable limit
] of these phenomena in the nervous system is suggestive. cycle. Furthermore, the oscillator’s bifurcation class has
S
Onelong-standinghypothesisholdsthatcorrelatedactiv- beenshowntodeterminetheshapeofit’sPRCandthere-
D
ity in the visual system underlies feature binding. Syn- foreit’sabilitytosynchronize. TypeIoscillatorsundergo
. chronous oscillations may also play a role in amplifying the saddle-node-on-an-invariant-circle, or SNIC, bifurca-
h
t signals [1], transmitting information from one layer to tion and the resulting PRC is strictly positive, indicat-
a
another [2–4], or such oscillations may encode informa- ing that perturbations can only advance the oscillator’s
m
tiondirectly[5–12]. Ontheotherhand,correlationsmay phase. Type II cells undergo the Andronov-Hopf bifur-
[ negatively impact the signal-to-noise ratio [13–16], and cation, which produces a PRC with both negative and
1 excessive synchrony is a hallmark of neurological disor- positive regions; typically, inputs occurring early in the
v ders such as epilepsy and Parkinson’s disease. cycle can delay the phase while later inputs advance it.
9 See Fig.(1).
To understand the function of oscillatory correlations,
1
9 or one day achieve clinically relevant control over them, An expanding body of work has demonstrated that
1 we must first understand the underlying biophysical over short time scales of less than one period, type II
. mechanisms. While synchrony can arise as the result of oscillators are more susceptible to stochastic synchrony
1
anatomical connectivity between neurons, much recent than type I. This has been shown via simulations and
0
1 work [17–22] has brought to light ways in which corre- in vivo [17, 25], by deriving the probability distribution
1 lated activity develops from the inherent stochastisicity ofthephasedifference[26], byminimizingtheLyapunov
: of neural systems. Thus, in the absence of direct cou- exponent of the phase difference [27], and most recently
v
pling, two or more neural oscillators may become syn- by calculating the spike count correlation over a range
i
X chronized by virtue of the statistical correlations in their of time windows [28]. The latter study further reports
r noisy input streams – a phenomenon we will refer to as that this finding reverses over long timescales, namely
a stochastic synchrony. that type I oscillators transmit correlations more faith-
For our analysis of stochastic synchrony, we appeal to fully than type II when observed over lengths of time
thetheoryofweakcoupling,whichholdsinthestochastic much greater than one period.
contextprovidedtheamplitudeofthenoiseissufficiently
In Section I we provide a brief introduction to the
small. Inparticular,anumberofgroups[18–20,23]have
phase reduction technique in a stochastic setting. Next
proved that the phase reduction technique [24] can be
inSectionIIweuseregularperturbationstogiveanovel
applied to oscillators receiving additive noise. Thus, we
and straightforward analysis of correlation transfer over
reduce a noisily driven oscillator to a scalar differential
long time scales. To facilitate our derivation, we use the
equation describing the evolution of the phase. This so-
total elapsed phase as a proxy for the spike count. Note
called phase equation depends only on the properties of
that the total phase (modulo the period) and the spike
thenoiseandtheoscillator’sphaseresettingcurve(PRC)
countdifferbyatmostone,whichisanegligiblequantity
which characterizes how small perturbations influence
when many spikes have been observed over a long time
the oscillator’s subsequent timing or phase.
window. Theexpressionwederiveforthecorrelationco-
Neural oscillators can be classified into two types ac- efficient of the total phase agrees both qualitatively and
2
(cid:68)(cid:72)Θ(cid:76)
along the limit cycle. In the case of a neural oscillator,
we assume the noisy perturbations arise as the result of
2
stochasticsynapticinput,whichinfluencesonlythevolt-
age variable. Hence Z(θ) has only one nonzero compo-
1 nent, which is proportional to the phase resetting curve
∆(θ).
Θ Thusfar,wehaveusedtheconventionalchangeofvari-
Π 2Π ables to obtain Eq.(1), which therefore must be under-
stood as a stochastic differential equation (SDE) in the
(cid:45)1 Stratonovich sense. In order to eliminate the correlation
betweenθ andξ wemustusetheItˆochangeofvariables,
FIG. 1. We use the parametrization ∆(θ) = −sin(θ+α)+ which will introduce an additional drift term:
sin(α)tovarythePRCsmoothlyfromtypeI(red),whereα=
π2 and ∆(θ)=1−cos(θ), to type II (blue), where α=0 and θ˙ =1+σ∆(θ)ξ+ σ2∆(cid:48)(θ)∆(θ).
∆(θ)=−sin(θ). Note that intermediate values of α produce 2
PRCshapes(dashedpurple)thatmorecloselyresemblethose
found empirically in vivo. Here (cid:48) denotes differentiation with respect to θ. For a
detaileddiscussionofphasereductioninnoisyoscillators
see [29].
quantitatively with the results found in [28].
In Section III we consider short time scales less than
II. CORRELATION TRANSFER OVER LONG
or equal to the period of the oscillation. In this case, TIME SCALES
the total phase cannot be used to approximate the spike
count. We therefore derive the spike count correlation
We now consider the transfer of correlations over time
directly, using simple probabilistic reasoning applied to
scales much larger than the natural period of the oscil-
the density of the phase difference. Our analytic results
lators. Given the level of correlation between the noisy
together with Monte Carlo simulations corroborate ear-
inputs,wewishtoknowwhatlevelofcorrelationremains
lier work showing type II oscillators transfer correlations
between the spike count of two oscillators after some
more readily than type I over short time windows.
time. For analytic convenience, however, we will use
the total phase that has elapsed as a proxy for the spike
count. Since these quantities differ by at most one, the
I. NOISY OSCILLATORS
discrepancy will be negligible for the large spike counts
that accrue over long time scales.
Letusbeginwithaneuraloscillatorreceivingadditive
Our system will consist of two identical phase oscilla-
noise with equations of motion given by
tors receiving weak, correlated, but not identical, addi-
tive white noise. Keeping only terms up to order σ, we
dX =F(X)dt+σξ,
have
whereX ∈Rn andξ isawhitenoiseprocess. Whenσ =
0, we assume the noiseless system has an asymptotically θ˙ =1+σ∆(θ )ξ (t)
1 1 1
stableperiodicsolutionX0(t)=X0(t+τ)withperiodτ. θ˙ =1+σ∆(θ )ξ (t). (2)
2 2 2
As in the deterministic case, we can reduce this high-
dimensionalsystemtoascalarequationfortheevolution
The noise takes the form
of the phase θ around the limit cycle. Let φ : Rn → S1
√ √
map a neighborhood of the limit cycle to the phase on ξ = c ξ + 1−c ξ
1 C A
a circle. That is, θ = φ(X), with θ ∈ [0,1). Then θ √ √
ξ = c ξ + 1−c ξ , (3)
satisfies 2 C B
dθ where ξ , ξ and ξ are mutually independent, zero
=1+σ∇ φ(X)·ξ, A B C
dt X mean white noise processes, and c ∈ [0,1] is the cor-
relation between ξ and ξ , which we will refer to as the
where we have normalized the unperturbed period to be 1 2
input correlation.
one. Next we can close the equation by assuming the
NextletusrewriteEq.(2)intheformofintegralequa-
noise amplitude σ is sufficiently small, so that the sys-
tions:
temtrajectorycanbeapproximatedbythenoiselesslimit
cycle X0: (cid:90) t
θ (t)=t+θ (0)+σ ∆(θ (s))ξ (s)ds
θ˙ ≈1+σZ(θ)·ξ, (1) 1 1 1 1
0
(cid:90) t
where Z(θ) = ∇Xφ(X0(θ)) is the adjoint, or phase- θ2(t)=t+θ2(0)+σ ∆(θ2(s))ξ2(s)ds.
dependent sensitivity of the trajectory to perturbation 0
3
A! ! ! ! ! ! ! ! ! B
P Φ
0.5
0.5
0.4
0.4
0.3 0.3
0.2 0.2
0.1 0.1
Φ
3 2 1 0 1 2 3 3 2 1 0 1 2 3
FIG. 2. The steady state distribution P(φ) of phase differences φ is shown for type I (red) and type II (blue) as well as for
intermediate PRCs (dashed purple). Note that the unperturbed period of the oscillators is 2π. (A) Input correlation c=0.4.
(B) Input correlation c=0.8.
Let T be length of the window of time over which ditions, say θ (0), to be distributed uniformly on the in-
1
we will observe the system. Throughout this discus- terval [0,2π]. However, at equilibrium the phases obey
sion we will assume that our system has reached equi- the steady state probability distribution P(φ) derived in
librium, and that time has been reparametrized so that [26]and[30],whichdependsonlyonthephasedifference
our observation takes place on the interval t ∈ [0,T]. φ(t) = θ (t)−θ (t). Therefore, the average of Eq.(5) is
2 1
In order to quantify the total phase traversed during computed as
this time, we subtract the initial phases by defining
q (T) = θ (T)−θ (0) for i = 1,2. Thus the total phase
i i i (cid:34) (cid:35)
(cid:90) T
traversed over a time window of length T is given by:
E[q (T)]=E T +σ ∆(s+x)ξ (s)ds
i i
0
(cid:90) T
q (T)=T +σ ∆(θ (s))ξ (s)ds. 1 (cid:90) 2π(cid:90) 2π
i i i = P(y−x)×
0 2π
0 0
(cid:34) (cid:35)
(cid:90) T
with q (0) = 0 for i = 1,2. Finally, since we assume σ
i T +σ ∆(s+x)(cid:104)ξ (s)(cid:105)ds dxdy
i
is small, let us simplify the integrands by expanding the
0
phase to lowest order:
σ (cid:90) 2π(cid:90) 2π
=T + P(y−x)×
2π
0 0
(cid:90) T
θ (t)=t+θ (0)+O(σ). (4)
i i ∆(θ (s))(cid:104)ξ (s)(cid:105)dsdxdy
i i
0
Then we have ∆(θ (s))=∆(t+θ (0)), and thus =T, (6)
i i
where 2π is the unperturbed period of the oscillators,
(cid:90) T
q (T)=T +σ ∆(s+θ (0))ξ (s)ds (5) P(φ) is the steady state probability distribution of the
i i i
0 phase difference, and x and y represent the initial condi-
tions θ (0) and θ (0), respectively. The last line follows
1 2
When taking expectations of the quantities in Eq.(5), because the white noises have zero mean.
wemustkeepinmindthattherearefourstochasticvari- Our goal is to compute the correlation of the total
ables over which averaging must take place. In partic- phase traversed by the two oscillators:
ular, we must average over the white noise signals ξ (t)
1
Cov[q ,q ]
and ξ2(t) and the initial conditions θ1(0) and θ2(0). Cor[q ,q ]= 1 2 . (7)
1 2 (cid:112)
Assuming we begin observation after the system has Var[q1]Var[q2]
reached equilibrium, we can take one of the initial con-
4
First, we derive the covariance as follows
Cov[q ,q ](T)=E[(q (T)−E[q (T)])(q (T)−E[q (T))]]
1 2 1 1 2 2
=E[(q (T)−T)(q (T)−T)]
1 2
(cid:34) (cid:35)
(cid:90) T (cid:90) T
=E σ2 ∆(s+θ (0))ξ (s)ds ∆(s(cid:48)+θ (0))ξ (s(cid:48))ds(cid:48)
1 1 2 2
0 0
1 (cid:90) 2π(cid:90) 2π (cid:90) T (cid:90) T
=σ2 P(y−x) ∆(s+x)∆(s(cid:48)+y)(cid:104)ξ (s)ξ (s(cid:48))(cid:105)dsds(cid:48)dxdy
2π 1 2
0 0 0 0
c (cid:90) 2π(cid:90) 2π (cid:90) T (cid:90) T
=σ2 in P(y−x) ∆(s+x)∆(s(cid:48)+y)δ(s−s(cid:48))dsds(cid:48)dxdy
2π
0 0 0 0
c (cid:90) 2π(cid:90) 2π (cid:90) T
=σ2 in P(y−x) ∆(s+x)∆(s+y)dsdxdy.
2π
0 0 0
Similarly, we find the variance to be c
out
Var[q ](T)=E[(q (T)−E[q (T)])2] 1
1 1 1
1 (cid:90) 2π(cid:90) 2π (cid:90) T
=σ2 P(y−x) ∆(s+x)2dsdxdy.
2π 0.8
0 0 0
Note that we therefore have Var[q ] = Var[q ], and
1 2
hence the denominator of Eq.(7) can be simplified: 0.6
(cid:112)
Var[q ]Var[q ] = Var[q ]. This gives the total phase
1 2 1
correlation as
0.4
Cor[q ,q ](T)
1 2
(cid:82)2π(cid:82)2πP(y−x)(cid:82)T ∆(s+x)∆(s+y)dsdxdy 0.2
=c 0 0 0 .
(cid:82)2π(cid:82)2πP(y−x)(cid:82)T ∆(s+x)2dsdxdy
0 0 0 !
(8) 0
0 !/4 !/2
(cid:82)2π
Nowleth(x)= ∆(y)∆(y+x)dybetheautocorrela-
0 FIG. 3. Output correlation for large time windows is shown
tionofthePRC,andletφ(t)=θ (t)−θ (t)representthe
2 1
asafunctionofthePRCshapeparameterα. Notethatwhen
phase difference as before. Then we can rewrite Eq.(8)
α=0thePRCisapuresinusoidandthereforetheoscillator
as
istypeII;whenα=π/2,theoscillatoristypeI(seeEq.(10)).
Theoreticalcurves(solid)areagoodmatchforboththesimu-
(cid:82)2π latedtotalphasecorrelation(dotted)andthesimulatedspike
P(φ)h(φ)dφ
c :=Cor[q ,q ](T)=c 0 . countcorrelation(starred). Colorsindicatethelevelofinput
out 1 2 (cid:82)2πP(φ)h(0)dφ correlation: 0.2(blue),0.4(green),0.6(red),0.8(cyan),0.99
0
(purple). In all cases, noise amplitude σ=0.05.
NotethattherighthandsidenolongerdependsonT af-
terweswitchedtheorderofintegrationandcanceledthe
resulting factors of T in both numerator and denomina-
tor. Nextwecandoawaywiththedenominatorentirely, where G(x)=1−c(h(x)/h(0)), and N is a normalizing
(cid:82)2π
since h(0) does not depend on φ, which leaves simply constant, N = 1/ 1/G(x)dx. Let us further define
0
the PRC to be
(cid:90) 2π h(φ)
c = P(φ)c dφ. (9)
out h(0)
0 ∆(θ;α)=−sin(θ+α)−sin(α), (10)
An expression for the steady-state probability density
of the phase difference P(x) was derived by Marella and
where α is a parameter that allows us to vary the PRC
Ermentrout in [26]. Specifically, we have
shape smoothly between type I (α = π/2) and type II
(α = 0). See Fig.(1). Using this, the phase distribu-
N tion over long time scales becomes a function of input
P(φ)= ,
G(φ) correlation and the PRC shape parameter:
5
c Our intuition for this finding can be honed by per-
out
forming a further perturbation expansion, now assuming
small input correlation. For sufficiently small c, we can
0.06 make the approximation
1 1 h(x)
0.04 = ≈1+c .
G(x) 1−ch(x) h(0)
h(0)
0.02 When we substitute this into Eq.(9) we find
Α N˜ (cid:90) 2π
c =c h(φ)dφ+O(c2), (14)
0 Π Π out h(0)
4 2 0
where N˜ = 1/(cid:82)2π(1+c h(x)/h(0))dx is likewise ap-
FIG. 4. The perturbation expansion of cout for small input 0 in
correlation (dashed) agrees well with the full output corre- proximated to lowest order in c.
lation (solid). Note that, to lowest order in c , the output TheformofEq.(14)demonstratesthatoutputcorrela-
in
correlation goes to zero as the PRC shape parameter α goes tion scales with the integral of the PRC autocorrelation,
to zero, that is, as the PRC shape approaches the pure type and for the parametrized PRC in Eq.(10) we have
II. Colors indicate the level of input correlation: 0.01 (blue),
0.05 (green), 0.1 (red).
(cid:90) 2π
h(φ)dφ=4π2sin(α)2.
0
In particular, α = 0 for the type II PRC, and hence
(cid:112) c = 0 to lowest order. Clearly, we have nonzero au-
(c−1)(cos(2α)−2)(2+(c−1)cos(2α)) out
P(φ;c,α)= . tocorrelation for nonzero α ≤ π, and hence PRCs that
2π(2−c+(c−1)cos(2α)−ccos(φ)) 2
deviatefrompuretypeIIwillproducehigheroutputcor-
(11)
relation over the long timescales considered here.
In the special cases where α=π/2 and α=0, Eq.(10)
ExpandingtheremainingtermsinEq.(14),wefindthe
and Eq.(11), together with Eq.(8), yield
approximated output correlation takes the form
2csin(α)2
Type I c = . (15)
out 2+c−(1+c)cos(2α)
∆ (x)=1−cos(x) In Fig.(4) we see that this approximation agrees with
I
√ √ Eq.(8) for c = 0.01 and 0.05 but diverges for c = 0.1.
3 c2−4c+3
P (φ;c)= (12) Note that these curves would all lie below the lowest
I 2π (3−2c−ccos(φ)) curve plotted in Fig.(3) if shown on the same scale.
1(cid:112) We verify the preceding analysis by simulating two
c =1− 3(c−3)(c−1)
out,I 3 phase oscillators perturbed by additive white noise as
described in Eq.(2) and Eq.(3). The simulations used
noiseamplitudeσ =0.05, andtheinputcorrelationtook
Type II
the values c∈{0.2,0.4,0.6,0.8,0.99}.
We computed the correlation coefficient of both the
∆ (x)=−sin(x) total phase and the spike count, using a range of obser-
II
√
vation windows T. As shown in Fig.(3), the total phase
1 1−c2
P (φ;c)= (13) correlation and the spike count correlation agree closely
II 2π(1−ccos(φ))
both with each other and with the theoretical curves as
(cid:112)
c =1− 1−c2 a function of the PRC shape parameter α.
out,II
III. SHORT TIME SCALES
As in [28], we see in Fig.(3) that type I oscillators dis-
play greater output correlation than type II oscillators Now we will calculate the spike count correlation di-
foranyfixedvalueoftheinputcorrelationc,asurprising rectly for observation windows T that are shorter than
finding in light of earlier results that demonstrated the or equal to the natural period, which we will assume to
opposite relationship over short windows of observation. be 2π. First let us consider the probability that a spike
6
functions of T:
P spike
f (T):=P[θ ≤2π−T,θ ≤2π−T]
1 f f 00 1 2
00 11 1 (cid:90) 2π−T (cid:90) 2π−T
= P(y−x)dxdy
2π
0 0
f (T):=P[θ >2π−T,θ ≤2π−T]
01 1 2
0.5 f f 1 (cid:90) 2π (cid:90) 2π−T
01, 10 = P(y−x)dxdy
2π
2π−T 0
f (T):=P[θ ≤2π−T,θ >2π−T]
10 1 2
1 (cid:90) 2π−T (cid:90) 2π
= P(y−x)dxdy
T
2π
0 2π−T
0 Π 2Π
f (T):=P[θ >2π−T,θ >2π−T]
11 1 2
1 (cid:90) 2π (cid:90) 2π
FIG.5. Jointspikingprobabilityfortwooscillatorsreceiving = P(y−x)dxdy.
2π
partially correlated noise is shown for observations windows 2π−T 2π−T
T ≤ 1, where 1 is the natural frequency of the oscillation.
Let X be the random variable such that X = 1 if θ
The subscripts ij indicate the probability that oscillator i or 1
j does (1) or does not (0) spike. spikes during the observation period T, and X = 0 if θ1
doesnotspike. Similarly, letY representthepresenceor
absence of a spike in oscillator θ . Then the covariance
2
is given by Cov[X,Y]=E[XY]−E[X]E[Y]. In terms of
the functions defined above we have
occurs in [0,T]. We say that oscillator i spikes when its
phase θi reaches 2π, which is to say θi(T) ≥ 2π. As- E[X]=0·(f00+f01)+1·(f10+f11)
suming as usual that the noise amplitude σ is small, we =(f +f )=E[X2]
10 11
expand the phase to lowest order as in Eq.(4), that is
E[Y]=0·(f +f )+1·(f +f )
θ (T)=θ (0)+T+O(σ). Thereforetheprobabilitythat 00 10 01 11
i i
oscillator i spikes is simply =(f01+f11)=E[Y2]
E[XY]=0·0·f +1·0·f +0·1·f +1·1·f
00 10 01 11
=f .
P[θ spikes]=P[θ +T ≥2π] 11
i i
P[θ does not spike]=P[θ +T <2π]. A few simplifications are possible. In particular, the
i i
sumf (T)+f (T)isjustthemarginalprobabilitythat
10 11
θ spikeswithintimeT. Sinceθ isuniformlydistributed,
1 1
For two oscillators, there are four possibilities for the this probability is simply T . Furthermore, we also have
2π
joint spike count: f = f by the symmetry of the density P, and hence
10 01
(cid:112)
Var[X]Var[Y] = Var[X]. Therefore the spike count
correlation over short time windows is
P[θ does not spike,θ does not spike]
1 2
=P[θ +T <2π,θ +T <2π]
1 2
Cor[X,Y](T;c) (16)
P[θ spikes,θ does not spike]
1 2
E[XY]−E[X]E[Y]
=P[θ1+T ≥2π,θ2+T <2π] =
Var[X]
P[θ does not spike,θ spikes]
1 2
f −(f +f )2
=P[θ +T <2π,θ +T ≥2π] = 11 10 11
1 2
(f +f )(1−(f +f ))
10 11 10 11
P[θ spikes,θ spikes]
1 2 f −(cid:0)T (cid:1)2
=P[θ1+T ≥2π,θ2+T ≥2π]. = T11(cid:0)1−2πT (cid:1)
2π 2π
1 (cid:20) (cid:90) 2π (cid:90) 2π (cid:21)
= 2π P(y−x)dxdy−T2 .
2πT −T2
2π−T 2π−T
These probabilities can be obtained directly by in- (17)
tegrating the density of the phase difference, Eq.(11),
over the appropriate domain. Note that this gives four Fig.(6A,B)showshowthisanalyticallyderivedoutput
discrete joint probabilities for each observation window correlationcompareswithnumericalsimulationsfortype
T ∈ [0,2π]. For convenience, let us define the following I and type II oscillators, respectively.
7
A! B! ! ! ! ! ! ! ! !
c
out
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2 0.2
0 0
0 ! 2! 0 ! 2!
Cc! D
out
0.5
0.5
0.4
0.4
0.3 0.3
0.2 0.2
0.1 0.1
T
0 Π Π 0 Π Π
2 2
FIG. 6. (A,B) Theoretical (solid) and simulated (dotted) output correlation curves are shown as a function of the observation
windowT ≤2π. (A)TypeIoscillators. (B)TypeIIoscillators. (C,D)Theinitialslope(dashed)ofthespikecountcorrelation
(solid)isthelinearapproximationofEq.(17)atT =0,whichisgiveninEq.(19). (C)TypeIoscillators. (D)TypeIIoscillators.
For all plots, noise amplitude σ=0.04, and colors indicate the level of input correlation: 0.2 (blue), 0.4 (green), 0.6 (red), 0.8
(cyan), 0.99 (purple). In all cases, noise amplitude σ=0.05.
We can make a further simplification by considering
the linear part of Eq.(17) for T close to zero:
(cid:32) (cid:33)
T c
c ≈
out,I π 3(1−c)+(cid:112)3(c−1)(c−3)
c
=T +O(c2)
6π
(cid:18) 1 (cid:19) T (cid:18) 1+c (cid:19)
cout =T P(0)− 2π +O(T2) cout,II ≈ 2π √1−c2 −1 (18)
c
=T +O(c2). (19)
2π
Thus, the initial slope of the output correlation is pro- From here, it is clear that the initial slope of c is
out
portional to the peak of the stationary phase difference greater for type II than for type I oscillators; in fact the
distribution, P(φ)| . Substituting P (0) and P (0) type II output correlation rises three times faster than
φ=0 I II
from Eq.(12) and Eq.(13), we obtain: the type I, to lowest order in c. See Fig.(6C,D).
8
A!c DISCUSSION
out
We have demonstrated a novel approach to approxi-
1
mating the spike count correlation of noisy neural oscil-
lators over both long and short time scales. In the case
0.8 of long windows of observation T much greater than the
natural period of oscillation, we used the total elapsed
phase,modulotheperiod,asaproxyforthespikecount.
0.6
The difference between these quantities is at most one
and hence is negligible for when many spikes are ob-
0.4 served over large time windows T. In our perturbation
expansion to lowest order in the noise amplitude, σ, the
correlationbetweenoscillatorsdependsonlyonthePRC
0.2 andthestationarydistributionofthephasedifference. A
further approximation assuming small input correlation
crevealsthatoutputcorrelationscaleswiththeautocor-
0
0 20! 40! 60! 80! 100! relationofthePRC,whichisanonnegativequantitythat
equals zero precisely when the PRC is a pure sinusoid,
i.e., when the oscillator displays type II dynamics. This
B! observationshedssomelightonthecounterintuitivefind-
ing, first reported by Barreiro, et al. [28], whereby type
1
I oscillators transfer correlations more faithfully than do
type II over long time scales, although the reverse holds
0.8 true for the better understood case of short time scales.
Usingstraightforwardprobabilisticreasoning,wecom-
puted the spike count correlation directly for short time
0.6
scales. In the limit of small T and small c, we obtain
an expression for the initial slope of the output correla-
0.4 tion, also known as the correlation susceptibility [8]. In
[8], de la Rocha, et al. use a phenomenological model to
explore the complex relationship between susceptibility,
0.2 firingrateandthresholdnonlinearities. Thepresentanal-
ysis illustrates the contribution of bifurcation structure
T viaphaseresettingdynamics. Inparticular,thesuscepti-
0
0 ! 2! 3! 4! bility is proportional to the peak of the stationary phase
difference distribution, P(φ)| , which in turn depends
φ=0
on the shape of the PRC.
FIG. 7. Output correlation is shown as a function of
intermediate-length observation windows T. Colors indicate Our analytic expressions in the limit of small noise
thelevelofinputcorrelation: 0.2(blue),0.4(green),0.6(red), agree well with spike count correlations computed from
0.8 (cyan), 0.99 (purple). (A) Type II oscillators (solid) ex- simulated oscillators. However, for tractability we in-
hibit higher output correlations over short time scales than cluded only terms of order one in the perturbation ex-
do type I (dashed). (B) This result reverses over long time pansion of the phase given in Eq.(4). As a result, the
scales. In all cases, noise amplitude σ=0.2. present analysis cannot account for the slow drift of the
correlation due to noise, which is visible for values of T
near 2π in Fig.(6), and is even more apparent for the in-
termediate values of T shown in Fig.(7). New analytic
methodscapableofaddressingnon-extremalcaseswould
shed light on this and many other questions in mathe-
matical biology.
[1] P. H. E. Tiesinga, Phys. Rev. E 69, 031912 (2004). [4] T. Tetzlaff, S. Rotter, E. Stark, M. Abeles, A. Aertsen,
[2] E. Salinas and T. J. Sejnowski, J. Neurosci. 20, 6193 and M. Diesmann, Neural Comput. 20, 2133 (2008).
(2000). [5] R. C. deCharms and M. M. Merzenich, Nature 381, 610
[3] A. Kuhn, A. Aertsen, and S. Rotter, Neural Comput. (1996).
15, 67 (2003). [6] J. M. Samonds, J. D. Allison, H. A. Brown, and A. B.
9
Bonds, J. Neurosci. 23, 2416 (2003). 045201(R) (2005).
[7] A.KohnandM.A.Smith,J.Neurosci.25,3661(2005). [19] J. N. Teramae and D. Tanaka, Phys. Rev. Lett. 93,
[8] J.delaRocha,B.Doiron,E.Shea-Brown,K.Josic, and 204103 (2004).
A. Reyes, Nature 448, 802 (2007). [20] H. Nakao, K. S. Arai, K. Nagai, Y. Tsubo, and Y. Ku-
[9] C.M.Gray,P.Knig,A.K.Engel, andW.Singer,Nature ramoto, Phys. Rev. E 72, 026220 (2005).
338, 334 (1989). [21] G.B.Ermentrout,R.F.Gala´n, andN.N.Urban,Phys.
[10] J. Biederlack, M. Castelo-Branco, S. Neuenschwander, Rev. Lett. 99, 248103 (2007).
D. W. Wheeler, W. Singer, and D. Nikoli, Neuron 52, [22] R.F.Gala´n,G.B.Ermentrout, andN.N.Urban,Phys.
1073 (2006). Rev. E 76, 056110 (2007).
[11] M. J. Chacron and J. Bastian, Journal of Neurophysiol- [23] K.YoshimuraandK.Arai,Phys.Rev.Lett.101,154101
ogy 99, 1825 (2008). (2008).
[12] K. J, E. Shea-Brown, B. Doiron, and J. de la Rocha, [24] Y. Kuramoto, Chemical Oscillation, Waves and Turbu-
Neural Comput. 21, 2774 (2009). lence (Springer-Verlag, 1984).
[13] E.Zohary,M.N.Shadlen, andW.T.Newsome,Nature [25] R. F. Gala´n, G. B. Ermentrout, and N. N. Urban, J.
370, 140 (1994). Neurophysiol. 99, 277 (2008).
[14] K. O. Johnson, Journal of Neurophysiology 43, 1793 [26] S. Marella and G. B. Ermentrout, Phys. Rev. E 77,
(1980). 041918 (2008).
[15] K.H.Britten,M.N.Shadlen,W.T.Newsome, andJ.A. [27] A. Abouzeid and G. B. Ermentrout, Phys. Rev. E 80,
Movshon, J. Neurosci. 12, 4745 (1992). 011911 (2009).
[16] W. Bair, E. Zohary, and W. T. Newsome, J. Neurosci. [28] A. K. Barreiro, E. Shea-Brown, and E. L. Thilo, Phys.
21, 1676 (2001). Rev. E 81, 011916 (2010).
[17] R. F. Gal´an, N. Fourcaud-Trocme, G. B. Ermentrout, [29] J. N. Teramae, H. Nakao, and G. B. Ermentrout, Phys.
and N. N. Urban, J. Neurosci. 26, 3646 (2006). Rev. Lett. 102, 194102 (2009).
[18] D. S. Goldobin and A. Pikovsky, Phys. Rev. E 71, [30] H. Nakao, K. S. Arai, and Y. Kawamura, Phys. Rev.
Lett. 98, 184101 (2007).
