PSD at harmonics of unconstrained stimuli

−π −π/2 0

(b)

π/2 π

χ 3,0(ψ,0|T o, M)

ψ 0

0.2 0.4 0.6 0.8

Figure 4.9: Joint spike phase densities for (a) the first and (b) the third spike in a train.

Colored histograms are count-conditional densitiesχk,0(ψ,0|To, M) for spike countsM = 5 (red),M = 8 (blue), andM = 12 (yellow). The black line is the plain jSPDχk,0(ψ,0). All parameters are identical to Fig. 4.7.

4.6.2 PSD at harmonics of unconstrained stimuli

For a sinusoidal stimulus of frequency Ω, the spike-count-conditional spectrum at the harmonics reduces to

S_To,M(nΩ) = M

πT_o[1 + 2Reh_M(nΩ)] (4.50) with

h_M(nΩ) = 1 M

M−1

X

j=1 M−j

X

k=1

e^{−in(ψj+k−ψj)}

(4.51) and the spike phases ψ_j as defined by Eq. (4.16). This average may be eval-uated within the framework of the Markov chain in phase, offering essential advantages. First of all, the Markov chain in phase is asymptotically stable and the stationary density χ^(s) easily obtained as shown above. The spectrum

72 Spike Trains

can thus be computed in the stationary state. Second, as demonstrated by a comparison of Figs. 4.9 and 4.7, the count-conditional joint phase density

χ_j+k,j(ψ, φ|T_o, M) = Prob

The j +k^th spike occurs at ψ and the j^th atφ, provided M spikes occur in (0, T_o].

(4.52) is approximated better by its “plain” counterpart χ_j+k,j(ψ, φ) than the corre-sponding spike time densities. The reason for this is that the global constraint on the spike count primarily alters the distribution of spikes between stimu-lus periods, but affects their distribution within periods but little. Plotting is facilited by χ_k,0(ψ,0|T_o, M) =χ_k(ψ).

Exploiting thatχ_j+k,j(ψ, φ) =χ_k(ψ|φ)χ^(s)(φ) in the stationary state, one may thus write to a good degree of approximation

e^{−in(ψj+k−ψj)} With an eye on numerical evaluation, the remainder of the derivation is given for the discrete Markov chain; the generalization to the continuous-phase Markov process should pose no principal difficulties. The discrete equivalent to Eq. (4.53) is

e^{−in(ψj+k−ψj)}

=ˆa(n)^tr·T^k·b(n)ˆ (4.54) with vectors (∆ψ = 2π/L is the phase discretization)

ˆ

a(n) = 1,e^−in∆ψ, e^−2in∆ψ, . . . , e^{−(L−1)in∆ψ}tr

ˆb(n) = χ^(s)₀ , χ^(s)₁ e^in∆ψ, . . . , χ^(s)_L−1e^{(L−1)in∆ψ}tr

. Upon inserting Eq. (4.54) into Eq. (4.51), observe that

e^{in(ψj+k−ψj)} de-pends only on the “distance”k of the spikes, but not on their “position”j due to stationarity. Hence the outer summation overj may be performed to obtain

h_M(nΩ) =ˆa(n)^tr where S^(M) is a diagonal matrix with elements¹¹

S^(M_mm⁾ =

11The corresponding Eq. (A5) ofPlesser and Geisel (1999) contains a typographical error.

4.6 Power spectral density 73

The λ_m are the ordered eigenvalues ofT (λ₁ = 1>|λ2| ≥. . .≥ |λL|), and the vectors a(n) and b(n) are given by

a(n) =C^trˆa(n), b(n) =C⁻¹ˆb(n) (4.57a) with the matrix decomposition

T=C·L·C⁻¹ , L= diag(λ₁ = 1, λ₂, . . . , λ_L). (4.57b) Note that the term for m = 1 has no counterpart in the constrained case, while the terms for m > 1 are very similar to the second and third terms of the power spectral density for the constrained case, cf. Eq. (4.47). Inserting Eq. (4.55) into Eq. (4.50) yields

S_To,M(nΩ) = M πT_o

1 + 2Rea^tr(n)S^(M)b(n)

. (4.58)

The matrix element S^(M₁₁ ⁾ is linear in the number M of spikes, and, since T_o/M → hτiasymptotically, this term gives rise to the singular spectral com-ponent that grows linearly with the length of the spike train. It is convenient to split the spectrum into its singular and bounded parts

S_To,M(nΩ) = M and the singular periodic contribution

B(n) = Rea1(n)b1(n) =

The last equality follows directly from Eq. (4.57), because the first column of Cis an eigenvector ofT to eigenvalueλ₁ = 1 and thus equal to the stationary phase distribution χ^(s) up to a factor c 6= 0. The first row of C⁻¹ as the corresponding right eigenvector has all elements equal to 1/c, yielding the result above. Note that A(n, M) is bounded for M → ∞, since |λm| <1 for m >1 from Section 4.4.2 so that

M→∞lim

74 Spike Trains

Figure 4.10: Power spectral density for unconstrained stimulation with the same sinusoidal stimuli used in Fig.4.8and observation timeTo= 200. Results at harmonics from Markov chain approach in blue, and spectra from simulated trains in red. Circles mark the estimated spectrum ST_o, triangles mark bounded part [1 +A(n, To/hτi)]/πhτi, and the filled square the limit of the spectrum at the origin from Eq. (4.49). The abscissa is in units of the stimulus frequency Ω and spectral power is in dB.

4.6 Power spectral density 75

Figure 4.11: Power spectral density at stimulus frequency S(Ω) from Markov chain (Eq.4.62) and simulated trains are in excellent agreement. (a) Scatter plot of results from 245 parameter sets with Markov chain result on the abscissa and simulation result on or-dinate; error bars for simulation results are shorter than the symbol size. The correlation coefficient is > 0.999. Observation time is To = 200. (b) Distribution of difference ∆S between simulation and Markov chain results, and fitted Gaussian distribution with mean h∆Si=−0.01 dB and standard deviationp

h∆S²i= 0.11 dB. The thick black bar marks two standard errors of the mean of the simulation results.

In analogy to the constrained case, the spike count is now estimated as¹² M = T_o/hτi ≈ hMi_To to obtain the time-limited PSD at the stimulus

Figure 4.10 indicates excellent agreement between the power spectral density at the harmonics given by this approximation with spectra obtained from simulated spike trains via discrete Fourier transform.

To quantify the quality of the approximation, S(Ω) was computed for a large number of stimuli from Eq. (4.62) and from simulated trains of 50000 spikes.¹³ The results are shown in Fig. 4.11 and indicate near-perfect agree-ment. The standard deviation of the difference ∆S between simulation and Markov chain results is only 0.11 dB, which is on the order of the statistical error of the simulation results, while the mean of ∆Sis nearly zero (−0.01 dB).

12Plesser and Geisel (1999) usedM =bTo/hτicto approximate the spike count averaging;

the precise value of the fraction improves the approximation slightly.

13Stimulus parameters were an evenly spaced subset of those used in Fig.5.17, p.104.

76 Spike Trains

0 2Ω 4Ω 6Ω

ω S T o(ω) [dB]

(a)

−25

−20

−15

−10

−5 0

0 2Ω 4Ω 6Ω

ω

(b)

−25

−20

−15

−10

−5 0

Figure 4.12: PSD at the stimulus harmonics increases with observation time for un-constrained stimulation, while the background is independent of T_o: (a) T_o = 250 and (b) T_o = 1000. Filter ringing is reduced at larger observation times. The stimulus is the same as in Fig. 4.10(c). The abscissa is in units of the stimulus frequency Ω and spectral power in dB.

The resulting spectra invite some comments:

• The first term of the power spectral density is again the Poisson train PSD.

• The stimulus is clearly discernible in the spectrum for unconstrained stimulation even for those stimuli where it has become invisible in the constrained case, compare Fig. 4.8.

• For long timesT_oand consequently large spike numbersM, the spectrum is dominated by the singular part∼M B(n), compare Fig.4.12. It arises from the persistent eigenmode to eigenvalue λ₁ = 1, and represents the long-range correlations between spikes induced by the synchronization of the spike train to the stimulus. One might call it the PSD of the stationary phase densityχ^(s), see Eq. (4.61). The peaks at the harmonics thus decay the faster, the wider the phase density. For perfect phase-locking, i.e. χ^(s)(ψ) =δ(ψ−ψ₀), the peaks would not decay at all.

• Since the Markov chain in phase does not describe the distribution of spikes across periods, it yields no information about the noise background of the spectrum. In particular, the continuous part of the spectrum

∼ [1 +A(n, M)] (triangles in Fig. 4.10) does not reliably estimate the noise level unless the neuron rarely skips a stimulus period without firing.

4.6 Power spectral density 77

• The dip of the spectrum at lowest frequencies is due to refractory effects as for the constrained stimulation. Since refractoriness should depend little on stimulation details, the limiting value S_∞(ω → 0) given by Eq. (4.49) applies to the spectra of unconstrained spike trains as well, as shown in Fig. 4.10.

• The derivation above applies to non-sinusoidal periodic stimuli if Ω is replaced by the pseudo-frequency ˆΩ. The resulting spectral values at multiples of ˆΩ will usually be of little interest.

• The full spectra displayed in Fig. 4.10 for comparison show noticeable filter ringing particularly for low stimulus frequencies. This is a conse-quence of the rectangular window [ 0, T_o] imposed in the definition of S_To. It is thus an artefact of the definition, not a numerical one. Ringing might be reduced by using windows with softer flanks (Press et al. 1992), but spectra would no longer be in accordance with the definition ofS_To, cf. Eq. (4.39). Since the evaluation of the power spectral densities at the harmonics in closed form would be severely complicated if, e.g., a Han-ning window were incorporated in the definition, only the rectangular window is used here. Errors at the spectral peaks are negligible provided that the observation timeT_o 2π/Ω. This condition is met throughout.

Im Dokument Aspects of Signal Processing in Noisy Neurons (Seite 85-91)