Data from auditory brainstem neurons in barn owls, chicks and American alligators were used to assess the effect of sampling on the calculation of VS. Animal husbandry and experimental protocols were approved by the Animal Care and Use Committee of the University of Maryland, the Regierung von Oberbayern (Germany), the University of Sydney Animal Ethics Committee, and/or the Marine Biological Laboratory (Woods Hole, MA, USA). Detailed procedures for surgery, stereotaxis, acoustic stimulus generation, and data collection have been provided by Carr and Köppl (2004) for owls, Köppl and Carr (2008) for chicks, and Carr et al. (2009) for alligators. In brief, animals were anesthetized and placed in a sound-attenuating chamber. Body temperature was maintained by a feedback-controlled heating blanket. An electrocardiogram was recorded via needle electrodes placed in the muscles of legs and/or wings to monitor muscle potentials and the heart beat. The head was held in a constant position by gluing a stainless steel head post and the skull was opened to expose the cerebellum. If necessary, a portion of the cerebellum was aspirated to expose the dorsal surface of the brainstem. Recordings were made with tungsten (2–20 MΩ) or glass electrodes (5–100 MΩ).

Custom-written software (xdphys, Caltech, CA, USA) was used for controlling acoustic stimuli and collecting data together with the TDT2 signal-processing system (Tucker Davis Technology, TDT, Gainesville, FL, USA). Acoustic stimuli were passed through a D/A converter (TDT DD1), filtered (TDT FT6-2), attenuated (TDT PA4), impedance-matched (TDT HB4) and delivered to the animal by earphones placed into the ear canals. Sound pressure levels were calibrated before recordings using built-in miniature microphones (Knowles EM3068, Itasca, IL, USA). Responses to acoustic stimuli were continuously monitored until the electrode reached the cochlear nuclei in the auditory brainstem (nucleus magnocellularis, NM; or nucleus laminaris, NL). After isolating a single unit, characteristic frequency (CF) and response threshold at CF were determined (Köppl and Carr, 2003). To measure the degree of phase-locking, continuous tones at or near the CF were presented with an intensity of 20 dB above the threshold. Signals from the electrode were amplified and filtered by a custom-built headstage and amplifier and passed through an A/D converter (TDT DD1), a threshold discriminator (TDT SD1) with an event timer (TDT ET1) and fed to the computer. In about half of the recordings, extracellular potential waveforms were stored to the computer and later analyzed. In other cases, only spike timing data generated by the level detector (TDT SD1) were stored. Both the potential waveforms and the spike timing data were digitized and stored at a sampling rate of 48077 Hz.

Custom-written Matlab (MathWorks, Natick, MA, USA) scripts were used for data analysis. For units with potential waveform data, spike timings t_j were calculated by peak detection (Figure 1) and VS was calculated according to Eqs. 1–3. For units without potential waveforms, stored spike timing data (which was generated by the threshold discriminator) was used to calculate VS. Note that no significant difference between data with and without potential waveforms was found in the results shown in Section “Examples From In Vivo Recording.” For each single unit, timing data from 400 to 10000 action potentials were stored. For Figure 6, we used timing data of 400 spikes from each unit recording to calculate VS.

To quantify the effect of sampling rate on VS calculation, potential waveforms or spike timing data were down-sampled with various sampling frequencies f_sample. Peaks tj′ of each downsampled waveform were detected and VS of the spike train was computed. For a unit without a stored waveform, downsampled spike timing tj′ was assigned by shifting each spike time t_j to the nearest sampling point after t_j and VS was calculated. In order to test significance of the phase preference, we calculated the significance probability for VS of each spike train by P = exp(−N(VS)²) with N being the number of spikes (Fisher, 1993). All the single unit data used in our analysis satisfy VS > 0.2 and N > 400, yielding P < 1.1 × 10⁻⁷.

In this subsection, we derive the theoretical upper and lower bounds of VS values with temporal sampling errors. We assume, for theoretical simplicity, that a sufficiently large number of spikes are collected and that the von Mises distribution (Fisher, 1993) can properly approximate the phase histogram of the spike trains.

Let g(x) be a periodic function with a period of 2π and be normalized as ∫−ππg(x)dx=1. The mean vector (X,Y) of the function g(x) is defined as: (4)X=∫−ππg(x) cosxdx, (5)Y=∫−ππg(x) sinxdx and the VS is: (6)VS=X2+Y2. The von Mises distribution is defined as: (7)g(x)=12π I0exp(kcos(x−m)), where k and m are the parameters determining the concentration and the mean phase, respectively. I₀ is the modified Bessel function of order zero satisfying I0=(1/2π)∫−ππexp(kcosx)dx and thus ∫−ππg(x)dx=1. By assuming m = 0 without any loss of generality, VS with the von Mises distribution can simply be calculated as: (8)VSexact=∫−ππg(x)cosxdx=12πI0∫−ππexp(kcosx)cosxdx. The subscript “exact” means that no temporal sampling error is incorporated in this calculation. An example is given in Figure 2A.

As discussed in the previous subsection, collected spike timing can be shifted within the length of the sampling window T = 1/f_sample. This temporal sampling error corresponds to a maximum phase error of ±πR. In the following text, R = f_signal/f_sample is referred to as the “sampling ratio.” The theoretical upper bound of the VS is obtained by assuming that all the spike timings are shifted in a biased fashion toward the direction of the mean phase of the original distribution to increase the value of VS (Figure 2B). In this case, the length of the mean vector of the shifted spike train is calculated as: (9)LU=∫−π+θ0g(x−θ)cosxdx+∫0π−θg(x+θ)cosxdx+∫−θθg(x)dx, where θ = πR = πf_signal/f_sample. The first, second, and third terms denote the contribution of the probability distributions on (−π,0), the distribution on (0,π) and the distribution concentrated at phase 0, respectively. The upper bound of VS is: (10)VSU=LU. The lower bound of VS can be obtained similarly but assumes that all the spike timings are shifted toward the opposite direction of the mean phase of the original distribution to decrease the value of VS (Figure 2C). In this case, the length of the mean vector of the shifted spike train is calculated as: (11)LL=∫−π−θg(x+θ)cosxdx+∫θπg(x−θ)cosxdx −∫−π−π+θg(x)dx−∫π−θπg(x)dx. The first, second, third, and fourth terms denote the contribution of the probability distributions on (−π,0), the distribution on (0,π), the distribution concentrated at phase −π, and the distribution concentrated at phase π, respectively. In contrast to the upper bound L_U, the value of L_L can be less than 0, since the “length” here is calculated with respect to the direction of the mean phase of the original distribution. A negative value of L_L means that the mean vector of the shifted spike train lies in the opposite direction of the original direction and in such a case VS can take an arbitrary value between 0 and VS_exact. Therefore we obtain the lower bound of VS as: (12)VSL=max⁡{0,LL}. The upper and lower bounds for five VS values ranging from 0.1 to 0.9 are shown in Figure 3 (dashed lines). The horizontal axis is the sampling ratio R = f_signal/f_sample. When the sampling ratio increases to 1, the upper bound of VS approaches to 1 and the lower bound to 0. This means that we cannot obtain a good estimate of VS if the sampling rate is as low as the reference stimulus frequency. Since the upper and lower bounds depend on VS_exact, we calculated the theoretical “maximum error” as max⁡0≤VSexact≤1{VSU−VSL}. Maximum VS errors calculated for several sampling rates are shown in Table 1. For R < 0.1, the maximum error is almost linear with R.

In the previous section, we obtained the upper and lower bounds of VS, assuming the von Mises distribution. Although these upper and lower bounds are of theoretical importance, it is practically unlikely that sampling is totally biased toward the direction where these limits are attained. In this section, we derive another estimate for error in VS by adopting the more natural assumption that collected spike timing is jittered randomly within the sampling window. Generally, this random sampling jitter flattens the spike distribution. Figure 4 shows examples of narrow (A), wide (B) and extremely wide sampling windows (C). Note that the length of sampling window (=1/f_sample) is converted to the length of the window function (=2πf_signal/f_sample, see next paragraph for detail). If the sampling window is small (or equivalently, if the sampling rate is high) compared to the reference signal, the effect of temporal sampling error is limited (Figure 4A). If the sampling rate is equal to the signal frequency, the temporal sampling error totally hides the temporal structure of the spike trains (Figure 4C).

Let g(x) be a periodic function with a period of 2π and be normalized as ∫−ππg(x)dx=1. In the following derivation, we do not need to assume any particular shape for g(x). Only a sufficiently large number of spikes are assumed to be collected to form the distribution function g(x). Since a spike occurred at phase x is assumed to be randomly shifted within the range of ±θ (θ= πR = πf_signal/f_sample), the distribution function h(x) of sampled spikes (Figure 4, gray areas) can be obtained as a convolution of the original distribution function g(x) (Figure 4, dashed lines) and a window function w(x) (Figure 4, insets). Precisely, (13)h(x)=(w∗g)(x)=∫−∞∞w(x−t) g(t)dt=12θ∫x−θx+θg(t)dt. The window function w(x) = 1/2θ (−θ < x < θ) and = 0 (otherwise). Since the Fourier transform of a convolution is the product of the Fourier transforms of the two functions, the mean vector (X_sampled, Y_sampled) of the function h(x) can be calculated as: (14)Xsampled=∫−ππ(w∗g)(x)cosxdx           =∫−ππw(x) cosxdx∫−ππg(x) cosxdx=(sinθθ)Xexact, (15)Ysampled=∫−ππ(w∗g)(x)sinxdx                        =∫−ππw(x)sin⁡xdx∫−ππg(x)sinxdx=(sinθθ)Yexact. Thus VS of sampled spike train is: (16)VSsampled=(sinθθ)VSexact=(sinπRπR)VSexact. Note that VS_sampled obtained here does not depend on a specific shape of the spike distribution g(x) whereas the upper and lower bounds discussed in the previous section were obtained only with the von Mises distribution.

We calculated VS_sampled for five VS_exact values ranging from 0.1 to 0.9 (Figure 3, solid lines). Although VS_sampled approaches to 0 when the sampling ratio R = f_signal/f_sample increases to 1, it is much more robust to R than the lower bound VS_L (Figure 3, dashed lines). Since VS_sampled = (sinπR/πR) VS_exact, the “expected error” of VS, defined as e_expected = (VS_exact − VS_sampled)/VS_exact can be calculated as: (17)eexpected=1−sin⁡πRπR. Expected errors with several sampling rates are shown in Table 1. Expected error is much smaller than the theoretically calculated maximum error (see also Figure 2), and is less than 2% if the sampling frequency f_sample is only 10 times greater than the signal frequency f_signal.

Figure 3 and Table 1 imply that the expected error increases quite slowly with the sampling ratio R for small R values. Using the Taylor expansion sinπR = (πR) − (πR)³/3! + O(R⁵), the expected error can be calculated as: (18)eexpected=1−sinπRπR=(πR)26+O(R4) The approximation e_expected = (πR)²/6 is 99.5% accurate for R < 0.1. This approximation explains the slow increase in the expected error to the sampling ratio.

In this section, we compare the expected VS errors obtained in the previous subsection with spiking data recorded in vivo. We use data from neurons in the nucleus magnocellularis (NM) and the nucleus laminaris (NL) in the auditory brainstem of owls, chicks, and alligators. These neurons show phase-locked spiking activity and play a key role in sound localization (Carr and Konishi, 1990; Köppl, 1997; Köppl and Carr, 2008; Carr et al., 2009). In our original data set, spike timing was collected with a sampling frequency of 48077 Hz. We downsampled the data with various sampling frequencies and re-calculated VS values (see Materials and Methods). Figure 5 shows the phase-locked activity of eight neurons with best frequencies ranging from 350 to 7000 Hz and with VS ranging from 0.27 to 0.82. In all the neurons shown, VS values decay according to the estimation given as VS_sampled= (sinπR/πR) VS_exact (Eq. 16), where the sampling ratio R = f_signal/f_sample.

The above result was entirely consistent with much larger data sets we have tested (Figure 6). Since VS_sampled = (sinπR/πR) VS_exact, we can estimate VS_exact = (πR/sinπR) VS_sampled. We use the data recorded at 48 kHz (original sampling frequency) to obtain the estimate value of VS_exact. In Figure 6, we plotted VS_sampled from downsampled spike data divided by estimated VS_exact. Decay of VS_sampled with the sampling ratio R is accurately predicted by the equation VS_sampled= (sinπR/πR) VS_exact. When the sampling rate f_sample is 20 times as large as the signal frequency f_signal (i.e., R = 0.05), VS_sampled can be predicted with a root mean square error of about 1%.

In this section, we examine the sampling effect on several circular statistics other than VS.