Given a set of reference spike trains \(x_(t)\), the corresponding set of burst spike trains \(^_(t)\), and a corresponding set of stimuli \(s_(t)\) for \(i = 1, \ldots , N_}\), we would like to relate the spectral statistics introduced in Section 2 for the original (reference) spike train and for the burst spike train (burst spikes added according to the algorithm introduced in the preceding section).

We start with the burst-induced change in the linear susceptibility, then continue with the relation for the second-order susceptibility, and finally derive the relation between the spike train power spectra with and without bursting.

The linear susceptibility describes how a weak time-dependent stimulus affects the time-dependent firing rate (the instantaneous mean value of the spike train). The magnitude of the susceptibility \(|\chi _1(\omega )|\) at a certain frequency can be interpreted as scaling the amplitude of the firing rate modulation, \(r(t)=r_0+\varepsilon |\chi _1(\omega )| \sin (\omega t-\phi )\), in response to a very weak sinusoidal stimulus, \(\varepsilon \sin (\omega t)\); of interest is here, whether stimulus-unrelated burst spikes may boost (\(|^_(\omega )|>|\chi _1(\omega )|\)) or merely diminish (\(|^_(\omega )|<|\chi _1(\omega )|\)) the linear response.

The next-order (nonlinear) response is characterized by the susceptibility \(^_(\omega _1,\omega _2)\) (or, for the reference spike train, \(\chi _2(\omega _1,\omega _2)\)) that depends on two frequency arguments and would describe the response up to second order in a small signal amplitude \(\varepsilon \). Also here, we would like to know the effect of additional burst spikes on the response.

Last but not least, we also aim at the power spectrum for the case without stimulus (spontaneous activity). This statistics (in combination with the linear susceptibility) is useful to compute a lower bound on the neural information transmission (with and without burst spikes).

First, we want to study the effect of burst spikes on the linear response function \(^_(\omega )\). Only the spike train is modified by the bursts (but not the signal) and hence only the cross-spectrum changes, \(S_ \rightarrow ^_\), yielding

$$\begin ^_(\omega )&= \frac^_(\omega )}(\omega )} \ . \end$$

(14)

The Fourier transform of the burst spike train, Eq. (12), is given by

$$\begin }^(\omega )&= \tilde(\omega ) + \sum \limits _^ \sum \limits _^^_} e^ + \sum \limits _^ ^_ \right) } \ , \end$$

(15)

which inserted in Eq. (5) yields the burst cross-spectrum:

$$\begin ^_(\omega )&= \lim \limits _ \frac(\omega ) \tilde^(\omega ) \right\rangle } + \nonumber \\&\quad \lim \limits _ \frac \left\langle \sum \limits _^ \sum \limits _^^_} e^ + \sum \limits _^ ^_ \right) } \tilde^(\omega ) \right\rangle \ . \end$$

(16)

Due to the additional randomness associated with the jittered IBIs and the burst-spike distribution, the brackets imply now not only an average over the intrinsic noise \(\xi \) and the stochastic signal s but also over \(^\) and \(^\): \(\left\langle \cdot \right\rangle = \left\langle \cdot \right\rangle _^,^}\). The first term in Eq. (16) is the cross-spectrum between stimulus and reference spike train. This leaves only the calculation of the second term. First, we carry out the average with respect to the IBIs. Since the random numbers \(^_\) are drawn independently for each burst spike (they form the local renewal process \(y_(t)\)), the average factorizes:

$$\begin \left\langle e^^ ^_} \right\rangle _^}&= \prod \limits _^ \left\langle e^^_} \right\rangle _^} \nonumber \\&= \prod \limits _^ \int \text ^_} \ e^^_} \rho \left( ^_\right) \nonumber \\&= \prod \limits _^ \varphi _^_} (\omega ) = \varphi ^_^}(\omega ) \ . \end$$

(17)

Here in the last line we used the definition of the characteristic function \(\varphi (\omega )\), the Fourier transform of \(\rho (^_)\), and that the \(^_\) are all drawn from the same distribution \(\rho \left( ^\right) \). For the example of the Gaussian, the characteristic function is well known:

$$\begin \varphi _}(\omega )&=e^\sigma ^}} \ . \end$$

(18)

To average over the burst spikes we again use that the number \(^_\) of burst spikes is drawn independently for each reference spike \(t_\) in each realization (trial). For the total number of realizations \(N_}\) we have for the k-th reference spike (assuming that it exists in every realization):

$$\begin \left\langle \sum \limits _^^_} \varphi _^}^(\omega ) \right\rangle _^}&= \frac}} \Bigl [ \left( \varphi _^} + \varphi _^}^ + \ldots + \varphi _^}^^_} \right) \nonumber \\&\qquad + \left( \varphi _^} + \varphi _^}^ + \ldots + \varphi _^}^^_} \right) \Bigr . \nonumber \\&\qquad + \ldots \nonumber \\&\qquad + \Bigl . \left( \varphi _^} + \varphi _^}^ + \ldots + \varphi _^}^^_}}} \right) \Bigr ] \nonumber \\&= \frac}} \left( p_N_}\varphi _^} + p_N_}\varphi _^}^ + \ldots \right) \nonumber \\&= \sum \limits _^ p_ \varphi _^}^(\omega ) \ . \end$$

(19)

Here we suppressed, for the ease of notation, the limit \(N_} \rightarrow \infty \) that should be taken for a proper ensemble average. We arranged the terms to illustrate how the probability \(p_\) to have at least n burst spikes, emerges. The latter probability can be calculated from the probability \(P_\) to have exactly j burst spikes as follows:

$$\begin p_ = \sum \limits _^ P_ \ , \end$$

(20)

and we may also easily invert this relation and write

$$\begin P_ = p_-p_ . \end$$

(21)

The remaining averages for the second term in Eq. (16) are now the same as for the cross-spectrum \(\left\langle \cdot \right\rangle = \left\langle \cdot \right\rangle _\), and we obtain for the burst cross-spectrum:

$$\begin ^_(\omega )&= S_(\omega ) + \lim \limits _ \frac \times \nonumber \\&\quad \left\langle \sum \limits _^ e^} \sum \limits _^ p_ \varphi _^}^(\omega ) \tilde^(\omega ) \right\rangle \nonumber \\&= S_(\omega ) + \lim \limits _ \frac(\omega ) \tilde^(\omega ) \right\rangle } \sum \limits _^ p_ \varphi _^}^(\omega ) \nonumber \\&= S_(\omega ) \left( 1 + \sum \limits _^ p_ \varphi _^}^(\omega ) \right) . \end$$

(22)

Therefore, we find for the linear response function with burst spikes Eq. (14)

$$\begin ^_(\omega )&= \chi _ (\omega ) \left( 1 + \sum \limits _^ p_ \varphi _^}^(\omega ) \right) \nonumber \\&= \chi _(\omega ) f (\omega ) , \end$$

(23)

which is the linear response function given by Eq. (7) multiplied by a frequency-dependent factor. The latter depends on \(\omega \) solely through the characteristic function, i.e. \(f(\omega ) = F(\varphi (\omega ))\), where \(F(\varphi )\) is a function of the characteristic function \(\varphi \). Using the Gaussian approximation for the IBI distribution, the factor reads

$$\begin f_}(\omega )&= 1 + \sum \limits _^ p_ e^\omega ^\sigma ^\right) } , \end$$

(24)

which assumes the form of a complex-valued damped oscillation with respect to the frequency argument \(\omega \). The “frequency” of this undulation has the physical dimension of a time, corresponding to multiples of the delay \(\tau \); the damping in turn is determined by the standard deviation of the jitter, \(\sigma \).

We finally note that there is also a different interpretation for the terms in \(f(\omega )\) that should also become apparent from our derivation above. For the nontrivial sum term we can write

$$\begin \sum \limits _^ p_ \varphi ^_^} (\omega )&= \sum \limits _^ \sum _^ P_ \int \limits _^ \text \ e^ \rho _(T) \\&= \sum \limits _^ P_ \sum \limits _^ \int \limits _^ \text \ e^ \rho _(T) \end$$

In the first line of the above equation, we have used the n-th order interval density \(\rho _(T)\); in the second line we have exchanged the sums and obtained for a given burst count n a sum over all n-th order intervals giving us the probability to obtain any spike after the reference spike. The outer sum then averages this over all possible total numbers of burst spikes.

We can further simplify the right hand side by exploiting the fact that the probability density of spiking after the (arbitrarily chosen) k-th reference spike is the ensemble average of the renewal train \(y_k(t_ + T)\), leading to

$$\begin \sum \limits _^ p_ \varphi ^_^} (\omega )&= \int \limits _^ \text \ e^ \left\langle y_(t_ + T) \right\rangle \nonumber \\&= \int \limits _^ \text \ e^ m_(T) = \widetilde_(\omega ) , \end$$

(25)

where \(m_(T)\) is the conditional firing rate for a burst spike in the k-th burst at time \(t_ + T\). In the last step of Eq. (25) it also becomes clear that for \(\omega = 0\) the Fourier transform, turning into a pure integral over the conditional rate, yields the full mean number of burst spikes (for one burst and without counting the reference spike). Furthermore, the factor \(f(\omega )\) can then be interpreted as the Fourier transform of

$$\begin \delta (T) + m_(T) , \end$$

(26)

i.e. the conditional firing rate within a burst which includes (by the delta function) the reference spike itself.

Next, we would like to study the effect of burst spikes on the nonlinear response function Eq. (8):

$$\begin ^_(\omega _, \omega _)&= \frac^_(\omega _,\omega _)}(\omega _) S_(\omega _)} . \end$$

(27)

As before, the power spectrum of the signal \(S_(\omega )\) is unaffected by the burst spikes, which leaves only the calculation of the third-order burst cross-spectrum \(^_(\omega _, \omega _)\). We obtain \(^_(\omega _,\omega _)\) by inserting Eq. (15) now evaluated at \(\omega \rightarrow \omega _ + \omega _\) in Eq. (6):

$$\begin ^_(\omega _,\omega _)&= \lim \limits _ \frac}^(\omega _ + \omega _) \tilde^(\omega _) \tilde^(\omega _) \right\rangle } \nonumber \\&= S_(\omega _,\omega _) \left( 1 + \sum \limits _^ p_ \varphi _^}^(\omega _ + \omega _) \right) . \end$$

(28)

We can directly write down the result for the third-order burst cross-spectrum, because only the spike train is affected by the burst spikes, and all steps from the calculation of the second-order cross-spectrum Eqs. (16)-(22) apply in the same manner. For the nonlinear response function with burst spikes we then obtain:

$$\begin ^_(\omega _,\omega _)&= \chi _(\omega _, \omega _) \left( 1 + \sum \limits _^ p_ \varphi _^}^(\omega _ + \omega _) \right) \nonumber \\&= \chi _(\omega _,\omega _) f(\omega _ + \omega _) \ , \end$$

(29)

which is the nonlinear response function given by Eq. (8) multiplied with the same frequency-dependent factor as in Eq. (23), evaluated now at \(\omega \rightarrow \omega _+\omega _\).

We note that for the Gaussian approximation the same applies as for the linear response: the factor \(f_}(\omega _ + \omega _)\) introduces a damped oscillation into the nonlinear response function. Furthermore, the factor will be constant along the antidiagonal \(\omega =\omega _ + \omega _=\text \), which suggests to consider a nonlinear response averaged over the anti-diagonal:

$$\begin \mathbb _} (\omega )\!&=\! \displaystyle \frac^ \text } \left| \chi _ (\omega _, \omega -\omega _) \right| }^ \text }}, & 0< \omega \le \omega _} \\ \displaystyle \frac}}^}}\!\!\!\!\!\! \text } \left| \chi _ (\omega _, \omega -\omega _) \right| }}}^}}\!\!\!\!\!\! \text }}, & \omega _}< \omega < 2\omega _} \end\right. } \end$$

(30)

Because of the projection on the summed frequencies, this function is considered in the interval \(\left( 0,2\omega _}\right) \), i.e. up to twice the cut-off frequency \(\omega _}\). Figure 3 illustrates how Eq. (30) comes about.

Projection of the nonlinear response function. The domain of the nonlinear response functions is limited by \(\omega _}\). We integrate over \(\left| \chi _(\omega _,\omega _) \right| \) along the anti-diagonals \(\omega = \omega _ + \omega _ =\) const and normalize these values by the length of the corresponding anti-diagonal. In the lower triangle we evaluate the projection for projection-frequencies \(0 < \omega \le \omega _}\) (red), and the upper triangle gives us the evaluation for the projection-frequencies \(\omega _}< \omega < 2\omega _}\) (blue)

The effect of burst spikes on the power spectrum of the spike train is complicated due to the fact that it involves second-order statistics of the spike train itself. Inserting Eq. (15) in Eq. (3) yields

$$\begin ^_(\omega ) =&\lim \limits _ \frac(\omega )\tilde^(\omega ) \right\rangle } \nonumber \\&+ \lim \limits _ \frac \left[ \left\langle \tilde(\omega ) \sum \limits _^ \tilde_^(\omega ) \right\rangle + \text \right] \nonumber \\&+ \lim \limits _ \frac \left\langle \sum \limits _=1}^\sum \limits _=1}^ \tilde_}(\omega ) \tilde_}^(\omega ) \right\rangle \nonumber \\ =&\lim \limits _ \frac(\omega )\tilde^(\omega ) \right\rangle } \nonumber \\&\hspace+ \lim \limits _ \frac \left[ \left\langle \tilde(\omega ) \sum \limits _^\sum \limits _^^_} e^ + \sum \limits _^ ^_\right) } \right\rangle + \text \right] \nonumber \\&\hspace+ \lim \limits _ \frac \left\langle \sum \limits _=1}^ \sum \limits _ = 1}^^_}} e^} + \sum \limits _=1}^} ^_,k_}\right) } \times \right. \nonumber \\&\left. \hspace \sum \limits _=1}^ \sum \limits _=1}^^_}} e^} + \sum \limits _=1}^} ^_,k_}\right) }\right\rangle \ . \end$$

(31)

The brackets indicate an average over the intrinsic noise \(\xi \), the IBIs \(^\) and the burst-spike distribution \(\left\langle \cdot \right\rangle = \left\langle \cdot \right\rangle _^,^}\). The first term is the power spectrum of the reference spike train. The averages over the IBIs and burst-spike distribution in the second term can be calculated in the same way as for the second-order burst cross-spectrum in Eqs. (17) and (19). Therefore, we obtain for the first and second term in Eq. (31):

$$\begin S_&(\omega ) \left( 1 + \sum \limits _^ p_ \left[ \bigl ( \varphi _^}(\omega ) \bigr )^ + \bigl ( \varphi _^}^(\omega ) \bigr )^ \right] \right) \nonumber \\&= S_(\omega ) \left( 1 + 2 \sum \limits _^ p_\, \text \^}^ (\omega ) \} \right) . \end$$

(32)

To evaluate the third term, we distinguish between the terms of the sum referring to the same reference spike (\(\sum _=k_}\)) or not (\(\sum _\ne k_}\)). First, we want to focus on the second case: because we are referring to different reference spikes, \(t_}\) and \(t_}\), the random numbers \(^_,k_}\) and \(^_,k_}\) are independent. Therefore, the average over the IBI distribution can be calculated independently:

$$\begin \left\langle e^=1}^} ^_,k_}} \right\rangle _^}&\left\langle e^=1}^} ^_,k_}} \right\rangle _^} \nonumber \\&= \bigl ( \varphi _^} (\omega ) \bigr )^} \bigl ( \varphi _^}^ (\omega ) \bigr )^} . \end$$

(33)

Furthermore, this also allows us to compute the average over the burst-spike distribution independently (Eq. (19)):

$$\begin&\left\langle \sum \limits _=1}^^_}} \bigl ( \varphi _^} (\omega ) \bigr )^} \right\rangle _^} \left\langle \sum \limits _=1}^^_}} \bigl ( \varphi _^}^ (\omega ) \bigr )^} \right\rangle _^} \nonumber \\&\qquad = \sum \limits _=1}^ p_} \bigl ( \varphi _^} (\omega ) \bigr )^} \sum \limits _=1}^ p_} \bigl ( \varphi _^}^ (\omega ) \bigr )^} . \end$$

(34)

It remains to calculate the exponential containing the spike times:

$$\begin \left\langle \sum \limits _\ne k_}^ e^} - t_})} \right\rangle&= \left\langle \sum \limits _^ (N-k) \left[ e^} + \text \right] \right\rangle \nonumber \\&= \left\langle \sum \limits _^ (N - k) \bigl [ \tilde_(\omega ) + \tilde_^(\omega ) \bigr ] \right\rangle \ . \end$$

(35)

Here we have rewritten the differences of the spike times by means of the k-th order interval \(T_ = t_ - t_\) and assumed that the average over \(\left\langle e^} \right\rangle \) does not depend on the spike-time index i; this is reflected by the suppression of the index in the notation of \(T_\) and by the prefactor \(N-k\) of the number of identical terms appearing in the sum. Furthermore, in the second line we used the fact, that the average of the phase factor \(\left\langle e^} \right\rangle \) for fixed k over different realizations of the intrinsic noise \(\xi \) will result in the Fourier transform of the k-th order-interval density \(\tilde_(\omega )\) (Holden, 1976). We keep the averaging brackets because the total spike count N is still a stochastic variable; in the limit \(T\rightarrow \infty \), the prefactor \((N-k)/T\) approaches the firing rate \(r_\) and we may omit the averaging brackets. Combining all steps, Eqs. (33)-(35), in consideration of the last sum in Eq. (31), we obtain for the terms referring to different spike times:

$$\begin&\lim \limits _ \frac\left\langle \sum _\ne k_} \sum \limits _=1}^^_}} \sum \limits _=1}^^_}} A^_} \right\rangle =\nonumber \\&r_\!\!\sum _ \bigl [ \tilde_(\omega )\! +\! \tilde_^ (\omega ) \bigr ] \!\!\sum \limits _=1}^ \!\!p_} \!\bigl ( \varphi _^} (\omega ) \bigr )^} \!\!\sum \limits _=1}^ \!\!p_} \!\bigl ( \varphi _^}^ (\omega ) \bigr )^}\!\!, \end$$

(36)

where \(A_ = \exp \left[ i\omega \left( t_ + \sum _^ ^_\right) \right] \). For the terms of the sum referring to the same reference spikes, we have to distinguish additionally between the terms of the sum referring to the same burst spikes (\(\sum _=n_}\)) or not (\(\sum _\ne n_}\)), yielding:

$$\begin \left\langle N \sum _=n_}^^_}} 1 + N \underbrace\ne n_}^^_}} e^=1}^} ^_,k_}} e^=1}^} ^_,k_}}}_} \right\rangle \end$$

(37)

The second sum \(\mathcal \) can be rewritten as follows:

$$\begin \mathcal = \left\langle \sum \limits _^^_}} \left( ^_} - n \right) \left[ e^^ ^_}} + \text \right] \right\rangle _^} \end$$

(38)

Note, that we used here the fact, that the IBI’s are drawn from the same distribution and only the length of the sequence and not the explicit index \(m_\) or \(m_\) is important. Evaluating the averages \(\left\langle \cdot \right\rangle _^,^}\) yields:

$$\begin&\left\langle \sum \limits _^^_}} \left( ^_} - n \right) \left[ e^^ ^_}} + \text \right] \right\rangle _^,^} \nonumber \\&= \left\langle 2\, \text \left\^^_}} \left( ^_} - n \right) \varphi _^}^(\omega )\right\} \right\rangle _^} \nonumber \\&= \left\langle 2\, \text \left\^} \left( \varphi _^}^^_}} - ^_} \varphi _^} + ^_} - 1 \right) }^} \right) ^}\right\} \right\rangle _^} \nonumber \\&= 2 \sum \limits _^ P_\, \text \left\^} \left( \varphi _^}^ - n \varphi _^} + n - 1 \right) }^} \right) ^}\right\} \ , \end$$

(39)

where we used in the second last step the result of the finite series

$$\begin \sum \limits _^ (M-n) a^ = \frac - aM + M - 1 \right) }} \ . \end$$

(40)

As for the other terms, the ratio N/T approaches the firing rate in the limit of an infinite time window T. With the above, we obtain for the terms of the sum referring to the same spike times:

$$\begin \lim \limits _ \frac&\left\langle \sum \limits _=k_} \sum \limits _=1}^^_}} \sum \limits _=1}^^_}}A^_} \right\rangle = r_ \sum \limits _^ p_ + \nonumber \\&\quad 2 r_ \sum \limits _^ P_\, \text \left\^} \left( \varphi _^}^ - n \varphi _^} + n - 1 \right) }^} \right) ^}\right\} \ . \end$$

(41)

Using the general result for a stationary spike train (Holden, 1976)

$$\begin S_ (\omega ) = r_ \left( 1 + \sum \limits _ \bigl [ \tilde_(\omega ) + \tilde_^ (\omega ) \bigr ] \right) \ , \end$$

(42)

and our results in Eqs. (32), (36) and (41), we obtain for the burst-spike-train power spectrum:

$$\begin ^_(\omega ) =&S_ (\omega ) \left[ 1 + 2 \sum \limits _^ p_\, \text \left\^}^ (\omega ) \right\} + \right. \nonumber \\&\left. \sum \limits _=1}^ p_} \bigl ( \varphi _^} (\omega ) \bigr )^} \sum \limits _=1}^ p_} \bigl ( \varphi _^}^ (\omega ) \bigr )^} \right] \nonumber \\&+ r_ \left[ \sum \limits _^ p_ \!-\! \sum \limits _=1}^ p_} \bigl ( \varphi _^} (\omega ) \bigr )^} \sum \limits _=1}^ p_} \bigl ( \varphi _^}^ (\omega ) \bigr )^} \right. \nonumber \\&+ \left. 2 \sum \limits _^ P_\, \text \left\^} \left( \varphi _^}^ - n \varphi _^} + n - 1 \right) }^} \right) ^}\right\} \right] \end$$

(43)

This can be further simplified to yield

$$\begin&^_(\omega )= S_(\omega ) \left| f(\omega ) \right| ^ + \nonumber \\&r_ \left[ \sum \limits _^ p_ \left( 1 + 2\, \text \left\^} \left( 1 - \varphi _^}^ \right) }^}}\right\} \right) - \left| f(\omega ) - 1 \right| ^ \right] \nonumber \\&= S_(\omega ) \left| f(\omega )\right| ^ + r_ g(\omega ) \ . \end$$

(44)

Here we have used Eq. (21), performed a few of algebraic manipulations, and expressed the sum term with the factor \(f(\omega )\) via

$$\begin \sum \limits _=1}^ p_} \bigl ( \varphi _^} (\omega ) \bigr )^}&= f(\omega ) - 1 \ . \end$$

Unlike the response functions we do not obtain in Eq. (44) the burst-spike-train power spectrum by a pure product of the reference spike-train spectrum and a frequency-dependent factor. Besides the original spectrum (first line) being multiplied with the squared absolute value of the factor introduced in Eq. (23), we have now also an additional term that includes the firing rate \(r_\) multiplied by a function \(g(\omega )\) of the burst characteristics.

When we inspect the high-frequency limit of the spectrum, it is useful to know the limits for the characteristic function \(\varphi _^}(\omega )\) and the factor \(f(\omega )\). In the case of jittered IBI’s with a smooth probability density, we can assume that

$$\begin \lim \limits _ \varphi _^}(\omega )&= 0 \quad \Rightarrow \quad \lim \limits _ f(\omega ) = 1 \ , \end$$

(45)

i.e. the probability density does not change in an infinitely fast manner, and hence its Fourier transform decays for very high frequencies.

The additive term in Eq. (44) ensures the saturation of \(^_\) in the large frequency limit at an increased firing rate

$$\begin \lim \limits _ ^_ (\omega )&= \lim \limits _ S_(\omega ) + r_ \sum \limits _^ p_ \nonumber \\&= r_ \left( 1 + \left\langle ^ \right\rangle \right) = ^_ \ , \end$$

(46)

which is then just the burst firing rate \(^_\). Notably, we recognize from Eq. (44) that without a jitter (when \(\varphi _}(\omega ) = e^\)) the burst-spike-train power spectrum will not saturate at \(^_\) but oscillate around this value.

From our derivation it is also evident that the additive term in Eq. (44) can never be negative,

$$\begin g(\omega )&= \sum \limits _^ p_ \left( 1 + 2\,\text \left\^} \left( 1 - \varphi _^}^ \right) }^}}\right\} \right) \nonumber \\&\quad - \left| f(\omega ) - 1 \right| ^ \nonumber \\&\ge 0 \ . \end$$

(47)

This inequality is not obvious but can be understood by considering the burst \(y_(t) = \sum \limits _^^_} \delta (t-t_)\), i.e. the finite (stochastic) number of burst spikes added to the k-th spike. Then its mean value is given by

$$\begin \left\langle y_(t) \right\rangle&= m_(t - t_) \ , \quad t > t_ \ . \end$$

(48)

The (one-sided) Fourier transform of this function is given by \(f(\omega ) - 1\). From our derivation it has become clear that the infinite sum in the first line of Eq. (47) is equal to the second moment of \(\tilde_(\omega )\), cf. the first line in Eq. (31) and in particular the terms in the last double sum with \(k_=k_\). Hence, the left hand side of Eq. (47) can be regarded as a variance

$$\begin g(\omega )&= \left\langle \tilde_(\omega ) \tilde_^(\omega ) \right\rangle - \left\langle \tilde_(\omega ) \right\rangle \left\langle \tilde_^(\omega ) \right\rangle \ge 0 \ , \end$$

(49)

which cannot be negative.

The seemingly exotic non-negativity of the additional term in the power spectrum is consistent with the insight that by adding burst spikes in a signal-unrelated manner to a spike train, we can only degrade the information that the spike train carries about the stimulus. At the level of a linear approximation this becomes apparent in terms of the coherence function between burst spike train and stimulus

$$\begin ^(\omega )&= \frac^_(\omega ) \right| ^}(\omega ) ^_(\omega )} = \frac^_(\omega ) \right| ^ S_(\omega )}^_(\omega )} \nonumber \\&= \frac(\omega ) \right| ^ \left| f(\omega ) \right| ^ S_(\omega )} \left| f(\omega ) \right| ^ + r_ g(\omega )} \nonumber \\&= \frac(\omega )}(\omega ) + r_ g(\omega )/\left| f(\omega ) \right| ^} C(\omega ) \ , \end$$

(50)

where \(C(\omega )\) is the coherence between the stimulus and the reference spike train. From the structure of the prefactor it is clear that as long as \(g(\omega ) \ge 0\), which was shown in the relations Eqs. (47)-(49), we have

$$\begin ^(\omega ) \le C(\omega ) \ , \end$$

(51)

i.e. we never increase the correlation coefficient between input and output by adding burst spikes to the spike train in a signal-unrelated manner. Likewise, we can conclude that by adding burst spikes in this way we cannot increase the lower bound of the mutual information rate, Eq. (10), which is a monotonic function of the coherence function. We will discuss these information-theoretic measures for the P-units in Section 6.3.