Revealing Spectrum Features of Stochastic Neuron Spike Trains

: Power spectra of spike trains reveal important properties of neuronal behavior. They exhibit several peaks, whose shape and position depend on applied stimuli and intrinsic biophysical properties, such as input current density and channel noise. The position of the spectral peaks in the frequency domain is not straightforwardly predictable from statistical averages of the interspike intervals, especially when stochastic behavior prevails. In this work, we provide a model for the neuronal power spectrum, obtained from Discrete Fourier Transform and expressed as a series of expected value of sinusoidal terms. The ﬁrst term of the series allows us to estimate the frequencies of the spectral peaks to a maximum error of a few Hz, and to interpret why they are not harmonics of the ﬁrst peak frequency. Thus, the simple expression of the proposed power spectral density (PSD) model makes it a powerful interpretative tool of PSD shape, and also useful for neurophysiological studies aimed at extracting information on neuronal behavior from spike train spectra.


Introduction
A sequence of action potentials is the main way neurons communicate each other. Neuroscientists agree that information is conveyed by the positions in the time of the action potentials, rather than by their waveform, so that often a sequence of unitary impulses located at the instant of the spikes onset, here referred to as neuronal spike train, is used as idealized neuron output [1].
Oscillations in neuronal spike train firing rate reflect encoding properties of single neurons and/or neuronal networks [2][3][4][5][6][7][8][9] and reveal different physiological and pathological conditions [10][11][12]. Indeed, this code reflects not only the neuron's input stimuli, but also the neuron's biophysical properties [1], so that one of the open questions in neuroscience is the evaluation of the extent to which a neuron's spiking response is due to the applied stimulus and to its intrinsic biophysical properties [13].
Estimation of the spike train power spectrum is a powerful method for analyzing neuronal outputs, used on a variety of neurophysiological recordings, such as multiunit activity (MUA), single unit spike trains and cellular membrane potentials [14][15][16][17][18][19][20], as well as on the output of different neuronal models to characterize the response of the modeled neuron to input stimuli [5,[21][22][23][24][25] Recently, spectral analysis was applied also to experimental recordings from Globus pallidum neurons of Parkinson's patients [18] and the general form of the power spectrum shape was suggested as a classification criterion of neural systems and their activity [18][19][20]. Power spectral analysis is also at the basis of several techniques aimed at the understanding of neuronal encoding properties of weak input signals [4,5,23,26,27]. Among them we may mention the classical signal to noise ratio (SNR), calculated as the ratio of the output power spectrum at the input frequency and the background noise spectrum [4,21,28,29], the coherence function between the input signal s(t) and the output y(t), calculated from second-order spectral measures of input and output [30], and advanced techniques based on homomorphic deconvolution [23,31].
Spectral estimates based on the Discrete Time Fourier Transform (DFT) and its fast implementation, the Fast Fourier Transform (FFT) are versatile and efficient techniques, largely used for evaluating the power spectral density (PSD) of neuronal spike trains coming from model simulations or experimental recordings. The typical shape of a neuronal spectrum presents a main peak which was interpreted as the mode of firing rate and less well-defined peaks approximately centered around the harmonics of the first one [32]. Our previous studies on spike trains coming from stochastic Hodgkin and Huxley (HH) neuron models [22,23,31,32] have shown that positions of such peaks depend, in a not trivial way, on the level of stimulation current and on the noise level, related to the size of membrane patch. Moreover, neuronal spectral features cannot be directly obtainable from statistical averages of the interspike interval (ISI).
To reveal such spectral features, we derive a model that expresses the PSD of a neuronal spike train as a series of expected values of cosines, function of frequency and multiples of the ISI random variable. The model, valid under the well-assessed assumptions that the ISI random variables are independent and identically distributed (i.i.d.), allows us to obtain a PSD model, which perfectly matches the FFT estimate obtained over one thousand neuronal traces. Moreover, the first term only of the series allows an accurate and straightforward estimation of the spectral peaks positions, giving a useful interpretation of the relationship between neuronal spectral components and ISI distribution.
In this paper, PSDs have been calculated and tested on the outputs of a stochastic HH model, under several conditions of background stimulation (stimulation current density) and channel noise (size of the membrane patch), but it is applicable to any kind of spike trains coming from different neuronal models, such as integrate-and-fire [33,34], multi-compartment models [35], or experimental recordings. The paper is organized as follows. In Section 2 the PSD model is derived. Results are shown and commented in Section 3. In particular, in Section 3.1, neuronal spectral features coming from the classical FFT based estimation of PSD are recalled and their dependence on the HH model parameters is analyzed. Section 3.2 presents the PSD model estimation of neuronal spectral features. Finally, in Section 4, conclusions are drawn.

Spike Train PSD Model
Substituting each action potential instance with a unitary impulse, located at the instant of the spikes onset, is a common practice to extract neuronal encoding from the neuron output. These instants are not deterministic but they can be considered as random variable (r.v.), therefore the spike train is a stochastic process, very often considered as a stochastic point process [1]. In particular, since data come from simulations or sampled neuronal recordings, we can consider the point process as a discrete time process, provided that at most one event occurs in each sampling period T s ; this condition is also called orderly. The point process represents the ordered times of occurrence of spikes as integer multiple of the sampling interval [36]. We can define the difference between two spike occurrences, called interspike interval (ISI), as where θ i is the time of a spike occurrence. Let us assume the r.v. Θ i as independent and identically distributed (i.i.d.). A point process (p.p.) may be defined formally as a random non-negative integer-valued measure [37]. For a p.p., the counting variate, u(0, t] counts the number of events in the interval (0, t]. An important elementary function of a p.p. is the differential increment. The differential increment, denoted by du(t), is defined as du(t) = u(t, t + dt], where in discrete time processing dt can be assumed equal to the sampling period, T s [36]. In the following we will use the current notation for discrete time signals, i.e., the normalised time, n = (t/T s ) and frequency, ω k = 2π k N , where N is the observation interval and is the floor operator that implements the integer division between t and T s .
For an orderly p.p., du[n] will take the value 0 or 1 depending on the occurrence of a spike in the sampling interval [36]; thus, for a given observation interval, the differential increment of the p.p. can be also written as where the r.v. θ i represents the normalized time occurrence of i-th spike, M is the number of spikes occurring in N, and δ is the discrete time unitary impulse. The DFT is an important means for the study of point processes. The DFT of a weakly stationary p.p. can be defined as [37] where ω k = 2π N k with k = 0 . . . N/2 are the discrete frequencies. The normalized discrete frequencies ω k can be expressed in hertz by doing f k = ω k 2πT s , like the normalized spike occurrence or ISI in second by multiplying them by T s . By substituting (2) in (3), we obtain (see also (25) in [36]) Finally, the power spectrum of the p.p. can be defined as This spectrum has a characteristic shape with one dominant peak around the mode of rate of firing and other less well-defined peaks centered around the harmonics of the first one [22,23,31,32,36].
If one wants to describe |U[k]| 2 as a function of ISI, the (4) can be further developed: From this formula one can notice that the final spectrum is the combination of cosine functions, dependent on the ISI distribution and the number of spikes of the realizations.
To clarify which terms make it up, the (6) can be written as The spectral representation in (7) is a stochastic process consisting of the sum of many components. Let us consider a generic term of (7).
Given a finite and different from zero value of ω k , the terms of the series are i.i.d. r.v., being a deterministic function of the i.i.d. random variables Θ i . In addition, the integer M is a r.v., but not independent from Θ i . Therefore M can also be written as a stopping time for the sum of Θ i : Under this assumption the validity of Wald's equation [38] holds, and we can write the expected value of Υ p (ω k ) as where , respectively the expected value and the variance of χ p (ω k ), Υ p (ω k ) properly shifted and normalized is which, considering a sequence with an infinite number of spikes, M, for the central limit theorem, tends toward a standard normal distribution N (0, 1). Eventually, since the sum of normal r.v. is a normal r.v. with expected value equal to the sum of the expected values, we can approximate the PSD of the spike train, taking only the first P terms of (7), obtaining the desired PSD model.

Test Data
The spike trains (du[t]) used as test were obtained from the complete stochastic HH model described in detail in [23,31]. Stochasticity of sodium and potassium currents was accounted for by solving a channel-state-tracking algorithm [39,40] with the Monte Carlo approach [41]. The bias current density I 0 , describing the excitability level of the neuron, was set to constant values (7,10,20) µA/cm 2 . All these values are considered supra-threshold, i.e., they induce a sustained firing mode even in the absence of internal noise [42]. The higher I 0 the more frequent the induced firing activity will be. Three values for the membrane patch size were considered: (11, 100, 300) µm 2 . Assuming constant spatial densities for the sodium (60 N Na /µm 2 ) and potassium (18 N K /µm 2 ) channels, patch area determines the channel noise level, which is inversely proportional to its square root [43,44]. The model was run in the C ++ environment using the explicit Euler integration method with time step T s equal to 10 µs. In our simulation conditions this numerical method is sufficiently accurate, although higher order methods may be necessary when considering time dependent I 0 [45]. The observation interval NT s is equal to 2 s, resulting in a spectrum frequency resolution of 0.5 Hz.

Spectral Features of Stochastic Neuron Spike Trains
The power spectra, obtained from the simulated spike trains by using the FFT, are shown in Figure 1 for the three current densities: I 0 = 20 µA/cm 2 (panel a), I 0 = 10 µA/cm 2 (panel b), I 0 = 7 µA/cm 2 (panel c). In each panel, spectra obtained for the three patch areas are compared with the deterministic behavior (purple line). As expected, except for the artifact induced by the spectral leakage, deterministic spectra were a train of pulses centered at the firing frequency and its higher order harmonics. The firing frequency increased with the input current density as reported in literature for an HH neuron model [34]. In the stochastic case, the main peak was not exactly centered at the deterministic firing rate and the other peaks did not correspond to the higher order harmonics of the first one. Moreover, as the frequency increased, they got wider and smaller and shortly disappeared for a high channel noise level, e.g., for the smallest patch area of 11 µm 2 . Figure 2 shows the differences of peak positions between stochastic and deterministic spectra, for all considered simulation conditions. For patches of 100 µm 2 and 300 µm 2 results for the first five peaks are reported in Figure 2, whereas, in the case of 11 µm 2 , only three peaks were faithfully detectable. Deviation of the main spectral peak from the deterministic firing rate increased with the channel noise and decreased with the input current density, i.e., the more the stochastic behavior prevailed on the deterministic input. It is worth noting that this deviation increased for the higher order peaks, faster and faster as the current density and the patch area decreased.

Estimation of Neuronal Spectral Features
The main results of Section 2.1 (Equation (12)) asserts that the frequency positions of spectral peaks depend on a linear combination of µ p (ω k ), the mean values of the function χ p (ω k ) = cos pω k Θ. By using the well known formula of the probability density of a non invertible function of random variable, the density of χ p (ω k ) can be written as where the TH q are the inverse images of x, and the summation is carried out over the Q k intervals in which the function χ p is monotone. The expression of f χ p (x) in (13) shows that it is very difficult to get analytically the trend of µ p (w), as well as to connect f Θ with the frequencies of spectral peaks.
So, in the following, we will investigate this connection starting from the analysis of f Θ (TH) and ending to deepen the PSD model (12). Figure 3 shows the densities f Θ per different values of I 0 : 20 µA/cm 2 (panel a), 10 µA/cm 2 (panel b), and 7 µA/cm 2 (panel c), and three patch areas: 11 µm 2 (blue lines), 100 µm 2 (red lines), and 300 µm 2 (yellow lines). We can notice that for higher I 0 values the main peaks were narrower and taller and centered at lower ISI values. Even the patch area affected the ISI distribution: the smaller the area the wider and smaller the main peak, centered at lower ISIs. For I 0 = 7 µA/cm 2 and patch area 300 µm 2 (yellow line in Figure 3c), the current density was inside the bifurcation range, (6.3-9.8) µA/cm 2 , and a reduced noise level determined the typical multimodal distribution, in line with results presented in [46]. To investigate the relationship between spectral features and the ISI distribution f Θ , in Table 1 we compared the first peak position in the PSD, f 1 PSD , with the mean firing rate λ = 1/E{Θ} [47] and with the inverse of the mode of Θ. Looking at Table 1, we can observe that they were almost coincident for the highest current density value (I 0 = 20 µA/cm 2 ) and the highest patch area (300) µm 2 , i.e., when the deterministic behavior dominated. In this case, the frequency of the first peak of µ 1 f 1 µ 1 , the mean firing rate λ, and 1/mode{Θ} could all be used as good estimates of the position of the main peak in the spike train spectrum. In the other cases, when the stochastic behavior became more relevant, λ underestimated the frequency of the main peak and 1/mode{Θ} overestimated it.
Except for the simplest case of deterministic models, the frequency of the main peak of the spectrum did not exactly correspond either to the inverse of the mean or the mode of the ISI distribution and the secondary peaks were not centered at multiple frequencies of the main one. Therefore, neuronal spectral features were hardly predictable and interpretable and could not be directly obtained from statistical averages of the ISI distributions.
The PSD model derived in Section 2 allowed us to interpret and evaluate in a very simple and direct way the position of the spectrum peaks.
Equation (12), derived in Section 2, furnishes an approximate model of the PSD of a spike train (S P u [k]), based on DFT and obtained as a linear combination of µ p (ω k ), expected value of χ p (ω k ), for p = 1, . . . , P. Figure 4, left column, shows probability density functions of χ p (ω k ), as a function of frequency, for p = 1 (panel a), p = 2 (panel c), p = 3 (panel e), in the case I 0 = 20 µA/cm 2 and patch area = 100 µm 2 . Panels on the right side of Figure 4 show the expected values of such random variables versus frequency, for the same values of p. As evidenced by looking at panels b, d, and f, µ p (ω k ) presents cosinusoidal terms which oscillated and decreased faster and faster with frequency as p index increased. As an example, for p = 3, µ 5 was almost completely damped for f > 300 Hz. Therefore, for high p values, the contribution to the summation of (12) was significant only in the low frequency range. Furthermore the sum of these µ p , which oscillated faster at low frequencies for high p values, contributed mainly to the almost constant plateau typical of the spike train spectra at lower frequencies than the main peak, whereas at highest frequency tended to zero.
If one considers only the first term (P = 1) and compares the PSD model S 1 (12) with the PSD calculated by FFT, as shown in Figure 5 for patch area = 100 µm 2 and I 0 = 20 µA/cm 2 , 10 µA/cm 2 , and 7 µA/cm 2 , significant differences were evident only in the low frequency range. Indeed, as shown in the Figure 4, in the low frequency range the µ p (ω k ) with high p were dominant. However, the peak positions almost coincided, as quantitatively shown in Table 2 for all simulated conditions and the first 5 peaks. For the smallest patch area (11 µm 2 ), only the first three peaks were reported as already mentioned in the comment to Figure 2. Table 2 reports the errors, in terms of difference between peak frequencies of FFT based PSD and µ 1 , as well as an average value calculated over all the peaks. It is worth noting that the difference was always below 2 Hz, except for the third peak in the case of the smallest area, i.e., when random fluctuations in the PSD may affect the automatic peak identification. (Looking only at the first peak, which was the most significant from a physiological point of view, the error was equal or lower than 2 Hz, also for the smallest patch).
If P increased, S P u [k] progressively captured all the shape of the PSD. Figure 6 shows comparison between the FFT based PSD and S P u [k] for P = 5 (panel a) and P = [M] − 1 (panel b), in the case of I 0 = 10 µA/cm 2 and patch = 100 µm 2 . When summation was limited to the first 5 terms (P = 5), except for fast oscillations below 60 Hz, the PSD model S P u [k] was a good approximation of the FFT spectrum; when all terms of summation were considered (P = [M] − 1), the FFT spectrum was very accurately captured in the whole frequency range. Furthermore, recalling, with the help of the graphics of µ p (ω k ) reported in Figure 4, that µ p (ω k ) tended to zero as frequency increased, the PSD model S P u [k] showed that the value of PSD at high frequencies was equal to M, the square root of the PSD DC value M 2 .

Conclusions
The PSD of a neuronal spike train has a characteristic shape consisting of peaks and valleys, whose shape and position depend on neuron features and input stimuli. This behavior can be observed in Figure 1, showing the spectra calculated on the outputs of a stochastic HH model where the current density takes into account a deterministic input, whereas the patch area determines the channel noise level. Similar PSD trends were obtained in [24,25] for simple leaky integrate-and-fire (LIF) neurons and in [48] for a sparse recurrent network of LIF neurons driven by stochastic inputs.
When the neuron exhibits a quasi deterministic behaviour, we may consider that the main peak is centered at a frequency equal to the mean firing rate (Figure 3). However, when the stochastic behavior of the neuron becomes significant with respect to the deterministic one, this is no longer true and the characteristic frequencies are difficult to be estimated and interpreted.
Moreover, as evidenced in [49], the detection of spectral peaks in spike trains could be impeded by several biological and experimental factors, so that analytical or semi-analytical approaches are advisable [48,49].
In this work, under the general and well assessed hypothesis that the ISI random variables Θ i , are i.i.d., we demonstrated that the S P u [k], given by (12) can successfully approximate the PSD of neuronal spike trains as shown in Figure 6b, which compares the S P u [k] with the PSD numerically calculated using FFT. S P u [k] consists of the linear combination of µ p (ω k ), the expected values of the random variables χ p , p = 1, . . . , P. These terms have a damped cosinusoidal trend versus ω that decays more and more rapidly as p increases. This means that terms corresponding to high p values can be neglected if the exact PSD shape is not needed and that the value of PSD at high frequency is equal to M, the number of spikes in the observation period.
Indeed, the maxima of the first term alone, µ 1 (ω k ), are able to estimate the positions of the spectrum peaks with an average error from 0.6 Hz to 1.5 Hz, depending on the model parameters, i.e., the input current density and the membrane patch area and with an absolute error averaged over all simulation conditions equal to 0.94 Hz. This allows us, not only to accurately calculate the frequencies at which the spectral peaks of a spike train are centered, but also to provide a theoretical interpretation of them.
Compared to literature results obtained using PSD analytic formulations on LIF neurons and networks driven by stochastic inputs [24,25,48], our PSDs of the neuronal spontaneous activity exhibit a damped quasi-periodic structure as in [24,48]. Recently, Droste and colleagues [25] showed that the damped PSD structure was strongly related to the features of the considered noise that, in our work, was realistically modeled as the endogenous channel noise. Similarly to the other analytic models present in the literature [24,25,48], even our PSD model gives an exact formulation of the spike-train power spectrum. In addition, it allows us to estimate the PSD features from the ISI statistics, up to the second order, and to give an interpretation of the typical PSD quasi-periodic structure. Indeed, by truncating the expansion of (12) to p = 1, it is possible an accurate, fast and straightforward calculation of the frequencies corresponding to the main and the higher order spectral peaks. Our PSD model demonstrates that such frequencies do not coincide with the mean firing rate and its multiples, as usually reported in the literature [48].
Although S P u [k] has been validated on the ISIs obtained from a stochastic HH neuron model, yet it is applicable to any measured or simulated trace as a quick and accurate estimate of the peak positions. Thus, the simple analytical expression of the proposed PSD model, with respect those available in literature [24,25,48,49], makes it a powerful interpretative tool of PSD shape, also useful for neurophysiological studies aimed at extracting information on neuronal behavior from spike train spectra. Moreover, given the general hypotheses underlying our PSD model, the proposed analysis may be expanded to the outputs of higher levels of biological complexity, such as neuronal networks, and may be applied even to experimental traces recorded through the patch clamp technique or multielectrode array (MEA).

Conflicts of Interest:
The authors declare no conflict of interest.