A bias in the traditional method for calculating PRCs
Somatic whole-cell patch-clamp recordings were made in current-clamp mode from spontaneously firing Purkinje cells in mouse cerebellar slices. To construct PRCs, a single brief depolarizing current pulse (amplitude, 0.05 nA; duration, 0.5 ms) was injected after a 100–150 ms baseline period (see ). Repeating this protocol many times should result in a homogenous sampling of phase space in spontaneously spiking neurons. The resulting change in interspike interval (ISI) relative to the mean ISI corresponds to the PRC value denoted by
. Plotting, for each trial,
as a function of the phase,
, at which the pulse arrived shows the overall ISI shortening corresponding to a positive PRC (, neuron firing at 180 Hz; see Materials and Methods
). Three observations can be made. First, at late phases there is a triangular region entirely void of data points (outlined in green) which we call the “Bermuda Triangle”. This causes a negative bias of the running average at late phases (, dashed red line). Second, the intrinsic variability in the ISI 
of spontaneously firing Purkinje cells acts as a source of noise, giving rise to data points with
. However, removing all points beyond 1 does not eliminate the negative bias (, solid red line). Finally, many trials (typically more than 5000) were required to calculate the Purkinje cell PRC reliably, while PRCs in other cell types are normally obtained from 100–200 trials 
. The ISI variability in Purkinje cells 
results in PRCs with low signal-to-noise ratio, increasing the bias at late phases and leading to a miscalculation of the PRC when this traditional method is used. Thus, a robust and unbiased method for the calculation of Purkinje cell PRCs in the presence of noise is required.
Purkinje cell PRCs determined using the traditional method.
To better understand the negative deflection of the PRC at late phases, a control PRC (cPRC; see Materials and Methods
) was calculated from the unperturbed voltage traces prior to the current pulse. The cPRC should be zero throughout all phases. However, the negative bias of the PRC at the late phases persisted in the cPRC (). We conclude that it is not the result of the brief current pulse injection. Rather, it results from the inhomogeneous sampling of the phase in the presence of noise. Indeed, the phase histogram (, lower panel) indicates that late phases are sampled less frequently.
To reproduce the effect of noise, PRCs were obtained from a Purkinje cell model 
in which Gaussian current noise was added to reproduce the irregularity of real Purkinje cell spiking (example model neuron firing at its spontaneous firing rate of 27 Hz; see Materials and Methods
). The model PRC exhibited the same negative deflection at late phases as observed in the experimental PRC (, dashed red line). As before, removing all points for which the phase exceeds 1 did not eliminate the negative deflection (, solid red line). Similarly, the cPRC in the model exhibited the same negative bias and the same inhomogeneous phase distribution () as the experimental cPRC. Therefore, the negative bias at late phases is a general feature of the traditional method for calculating PRCs, and must be due to the intrinsic ISI variability.
In order to explain how the ISI variability might affect the PRC calculation, we sketch twelve representative scenarios in which spike jitter due to noise causes misclassification of the phase and/or the PRC value. In these scenarios (shown in ), we jittered either the first or the second AP (, black lines) with respect to a perfectly periodic cycle of firing (, grey lines). We divided the sketches into three blocks depending on the phase of the current pulse within the cycle (: early phase; B, E, H, K: middle phase; C, F, I, L: late phase). The misclassification of phase and/or PRC value (arrows) becomes clear when comparing them against their deterministic counterparts. The jitter of the spike preceding or following the brief current pulse can lead to a loss of causality and hence to a drastic miscalculation of the PRC. The most serious consequences of the ISI variability due to noise occur in the scenarios illustrated in , where the jitter causes the current pulse to fall into a different cycle of firing, resulting in a significant bias at the late and early phases of the PRC, respectively. Specifically, the “Bermuda Triangle” effect present in both model and experiment can be explained by means of the sketch in : when the pulse arrives at late phases, and the AP jitter results in the pulse falling into the subsequent ISI as compared to the deterministic case, the resulting phase is small according to the new ISI boundaries. Due to causality, it is impossible for a PRC point to fall into the green “Bermuda Triangle” in , since for all points in the triangle the shortening of the ISI would be larger than the actual phase difference of the pulse to the end of the ISI. This explains the observation that phases are sampled less frequently in the late part of the ISI, and thus the PRC values are underestimated and the effect of ISI noise is not averaged out.
Interspike interval variability causes a bias in the traditional method to calculate PRCs.
To visualize the resulting phase and PRC misclassification, we translated each of these twelve sketches onto a corresponding phase plot (). This allows the resulting phase and PRC values of each of the twelve cases to be compared against their deterministic counterparts. More specifically, regularly spaced spike times were defined and jittered independently by noise taken from a Gaussian distribution. The known actual phase without noise was plotted against the sampled phase. The assumption that the process underlying spiking is perfectly periodic and that the presence of a spike does not reset this underlying process is made only for generating the data in (and subsequently ), and only for purposes of illustration. In a purely deterministic scenario, the sampled phase is linearly dependent on the actual phase (, points on the diagonal). This is also the case for occurrences in which the noise has no effect on the phase (e.g. the scenarios in or ; yellow points in ). For any deviations of the sampled phase from the actual phase due to noise, the points are scattered across the plot (, color coding as in A–L). Based on the same principles, the effect of noise in each of the twelve scenarios on the PRC plot is shown in (color coding as in A–L).
A new approach for generating corrected PRCs.
To summarize, the bias at late phases of the PRC calculated using the traditional method is due to erroneous phase sampling, which results from the substantial ISI variability present in spontaneously firing Purkinje cells, and the loss of causality between the current pulse and the jitter in the times of either of its two surrounding APs.
Improving the traditional method to obtain PRCs in the presence of noise
Our new method to correct for the bias in the traditional PRC and obtain a homogeneous phase histogram is illustrated in . The red spike immediately preceding the pulse is the one used as a reference (
) in the traditional method. In our new method, instead of using just the spike immediately preceding the current injection, each spike in the spike train is taken as a reference one at a time and the corresponding phase values (indicated under the arrows in ) are all taken into account (see also Materials and Methods
). In this case, the two spikes prior to the stimulation pulse (red and black in ) predominantly contribute to the phase interval [0,1] of the PRC (, red and black points). The impact of the pulse on the subsequent ISI, the PRC2
, is then determined by the two spikes following the current pulse ( blue and cyan spikes; also compare 
) and so on. It is worth emphasizing that even though more than one spike is included in the PRC calculation, the presence of each reference spike resets the phase to zero (
). Our method restores periodicity in the spiking jitter as can be seen in (all points, in analogy to ). By taking only the points according to the traditional calculation of the PRC, a sharp boundary is drawn (, red) resulting in an inhomogeneous distribution of sampled phases (, red). In contrast, by including the second spike prior to the pulse, spike jitter effects are averaged out (). The bias at the late phases of the PRC plot observed when taking points according to the traditional calculation of the PRC (, red) is thereby eliminated (, all points), as is the bias in the cPRC (not shown).
In order to validate our new method, we applied it to neuronal models for which the PRC can be calculated analytically (from the adjoint 
). PRCs of the Morris-Lecar model (parameters from 
), in the presence and absence of noise, were compared with the analytically derived PRC (Fig. S1A
). The PRCs calculated using both the traditional and our corrected method overlap perfectly (except near
, due to the finite time step and finite amplitude of the current pulse in the simulations), and match the analytically derived PRC. In the presence of noise, the PRC calculated by the traditional method is biased at late phases, as described above. Our new method eliminates most of this bias.
However, it has been shown that noise can directly affect the dynamics of neurons underlying the PRC, leading to changes in the PRC which are not due to measurement errors (e.g. in the Morris-Lecar model 
). We therefore used an additional model, the non-leaky integrate-and-fire model, in which noise-dependent changes of dynamics can be excluded. When noise was introduced in this model, the traditional method resulted in a biased PRC, as compared to the analytically derived PRC and the PRC in the absence of noise. Again, our corrected method removed most of this bias (Fig. S1B
). The same analysis was repeated in a leaky integrate-and-fire model. This shows that the “Bermuda Triangle” and its consequences on the PRC are the result of the traditional calculation of PRCs, separate from the effect of noise on the dynamics of the system (Fig. S1C
Next, using the Purkinje cell model 
, we compare the result of our method (, black) to the PRC obtained with the traditional method (, red) and the deterministic PRC without noise (, green). When the noise is increased, reflected by an increased coefficient of variation (CV) of ISIs, the traditional PRC deviates from the deterministic one and the bias becomes more pronounced (, dashed red line). In contrast, our corrected method performs as well as with low CV (, dashed black line). The strong bias at late phases is eliminated. In order to evaluate the performance of our method in comparison to the traditional method, we calculated the integral of the differences between PRCs and their deterministic counterparts (PRC error; ). As the CV increases, the PRC error shows larger increases using the traditional (, red line) compared to our corrected method (, black line).
Validation of the corrected method for obtaining PRCs.
In conclusion, the “Bermuda Triangle“ present in PRCs is due to shortcomings of the traditional method for calculating PRCs. The bias can for the most part be compensated for by taking the two spikes preceding the pulse as a reference, one at a time, instead of just the spike immediately preceding the pulse as in the traditional method.
A frequency-dependent switch in Purkinje cell dynamics
Spontaneous firing frequencies of Purkinje cells range from 10–180 Hz both in vitro 
and in vivo 
. To test how the dynamics of Purkinje cells change according to the firing frequency, we recorded from cells firing spontaneously at low (15–40 Hz, n
10) and high (55–180 Hz, n
6) rates and calculated their PRCs using our corrected method.
A representative corrected PRC is shown in (the same example of a rapidly firing (180 Hz) Purkinje neuron as in ). The PRC is positive, indicating that the brief current pulse causes an advance of the following spike (shortening of the ISI relative to the mean) with maximum displacement when the pulse arrives near the middle of the ISI. It is worth noting that the phase histogram is homogeneous (, lower panel), suggesting that, with the corrected method, the ISI is equally sampled throughout. In order to study the effects of the brief pulse on the subsequent intervals we plotted the PRC2–5
(; see Materials and Methods
is negative, suggesting that the subsequent ISI is lengthened relative to the mean. A PRC2
with opposite sign to the PRC has been previously reported 
and it is believed to be due to a compensatory effect on the current ISI length. Indeed, as seen from PRC2–5
, these curves are negative and the effect dies out after about 4 ISIs.
Two types of PRCs depending on the Purkinje cell firing rate.
In comparison, an example of a PRC of a slowly firing (30 Hz) Purkinje neuron is shown in . The brief current pulse causes the same positive displacement of the following spike independently of its position within the ISI, resulting in a square PRC. The phase histogram is homogeneous, indicating that there is an equal probability for the pulse to arrive at each phase within the ISI (, lower panel). In order to study the effects of the brief pulse on the subsequent intervals we calculated the PRC2–5
(; see Materials and Methods
). They were negative, similar to those of cells firing at a high rate, but exhibited larger fluctuations. It is interesting to note that the PRC phase advances occur at a different scale in the slowly and rapidly firing Purkinje cells. However, when converted back into time units, the PRC values are of the same order of magnitude in both cases (see below).
The PRCs of Purkinje cells exhibiting slow (15–40 Hz; n
10) and rapid (55–180 Hz; n
6) spontaneous firing were calculated using our corrected method. The PRCs switched from square (phase-independent) for lower frequencies () to phase-dependent for higher frequencies (). The switch occurred at a frequency of approximately 50 Hz. The average PRC of all neurons firing at low rates (, thick line) is phase-independent. To our knowledge, such a square PRC has not been previously reported. A square PRC can only be obtained if the cells act as perfect non-leaky integrators. In contrast, the average PRC of all Purkinje cells firing at high rates (, thick line) exhibited a sharp peak. It is useful to compare these average PRCs (Fig. S2
, thick black and red lines) with the biased ones obtained with the traditional method (Fig. S2
, thick green lines). To quantitatively assess the switch in dynamics we plotted the peak-to-baseline ratio of the PRCs in relation to the firing rate (; see Materials and Methods
). This quantity essentially compares the extreme value in the first half of the PRC with the extreme value in the second half. The switch at a firing rate of approximately 50 Hz can be seen clearly in this representation. The switch becomes particularly apparent when both the phase and the phase shift of the PRC are plotted in units of time, and phases are aligned with respect to the second AP in the ISI (). Then, the peaks of the PRCs measured at high firing rates coincide (red), indicating that an input signal causes an effect only in a 3 ms window prior to the output spike irrespectively of the precise firing rate of the cells in that group. This peak in the PRC is shown to give way to a larger phase-independent plateau (black) at low firing rates, in which incoming signals will affect the spiking of the cell regardless of the time at which they arrive. A transitory PRC (thin solid red lines in and , indicated by arrows) showing both a plateau at early phases and a peak at late phases was observed in a cell with intermediate firing frequency (55 Hz). To summarize, the PRCs of Purkinje cells largely depend on the intrinsic firing frequency of the cells: they are phase-independent at low firing rates (15–40 Hz) and phase-dependent at high frequencies (55–180 Hz).
A frequency-dependent switch in Purkinje cell dynamics.
The firing rate of a Purkinje cell changes depending on modulation of its inputs 
. For example, during locomotion in cats the firing frequencies of Purkinje cells can increase from an average of about 40 Hz 
to more than 100 Hz 
. To test whether the switch in Purkinje cell dynamics can occur in the same cell, we recorded Purkinje cell PRCs while modulating their firing frequencies using injected current (n
3; , points labeled with two colors). We first recorded at the spontaneous firing frequency, and if the spontaneous frequency was low, we next increased the firing rate by injecting a positive constant current. Alternatively a negative constant current was injected if the spontaneous frequency was high. The PRCs for both fast and slow states were calculated (, color coding as in ). When Purkinje cell spiking was changed from slow (33 Hz) to fast (104 Hz), the originally square PRC (), exhibited a sharp peak (). This change in the PRC was reversible, as when the neuron was allowed to relax back to its intrinsic firing rate (40 Hz) the PRC returned to a square shape (). Conversely, another neuron initially firing at a high rate (71 Hz) exhibited a peaked PRC (), which was switched to a square shape by reducing its firing rate to 26 Hz via injection of hyperpolarizing current (). When the neuron was then allowed to fire at its intrinsic firing rate (84 Hz) the sharp peak in the PRC reappeared (). Therefore, the switch in Purkinje cell dynamics reflected in the switch of the PRC can also occur in the same cell.
The switch in PRC shape can occur within the same cell.