Introduction {#s1}
============

Cerebellar Purkinje cells exhibit a wide range of dynamical phenomena. They are intrinsic neural oscillators, firing spontaneously and highly rhythmically in the absence of synaptic input, at a rate of 10--180 Hz [@pcbi.1000768-Husser1]--[@pcbi.1000768-Raman1]. They also exhibit intrinsic bistability [@pcbi.1000768-Williams1], [@pcbi.1000768-Loewenstein1], which influences their responses to sensory stimulation [@pcbi.1000768-Loewenstein1]. In addition, interactions between spontaneously firing Purkinje cells can result in waves of activity travelling down the cerebellar folia [@pcbi.1000768-Watt1], or in high frequency oscillations [@pcbi.1000768-Maex1], which may contribute to the generation of precise temporal patterns in the cerebellar network [@pcbi.1000768-deSolages1]. Hence, the firing of Purkinje cells is highly time- and state-dependent, and thus they represent excellent targets for dynamical systems analysis.

The phase response curve (PRC; [@pcbi.1000768-Reyes1]--[@pcbi.1000768-Canavier1]) is a powerful tool to study neuronal dynamics at the cellular level. The PRC describes the effect of a brief perturbation on the firing phase of a neuron, and can be used to predict the response of a neuron to more complex stimulation patterns [@pcbi.1000768-Reyes1]--[@pcbi.1000768-Canavier1]. The shape of the PRC is linked to the type of neuronal excitability [@pcbi.1000768-Izhikevich2], [@pcbi.1000768-Galn1], to oscillatory stability [@pcbi.1000768-Tateno1] and to network synchronization properties [@pcbi.1000768-Ermentrout1]--[@pcbi.1000768-Achuthan1]. Studying Purkinje cell PRCs is therefore an essential step to explore their dynamic repertoire, probe their biophysical mechanisms, and to construct models of Purkinje cells to determine their role in information processing at the network level.

PRCs can be obtained by directly perturbing the membrane potential by short (infinitesimal) square current pulses [@pcbi.1000768-Reyes1]--[@pcbi.1000768-Canavier1] or synaptic conductance pulses [@pcbi.1000768-Canavier1], [@pcbi.1000768-Achuthan1]--[@pcbi.1000768-Maran1], and via indirect methods [@pcbi.1000768-Ermentrout3]--[@pcbi.1000768-Ota1]. Using the direct method, infinitesimal PRCs are obtained by repeatedly injecting brief current pulses while neurons are firing action potentials (APs). Phase and phase perturbation are measured by using the AP immediately preceding the current pulse as a reference, and we refer to these PRCs as "traditional" PRCs throughout this paper. We show using electrophysiological experiments and in simulations that the interspike interval variability present in Purkinje cells introduces a systematic bias in this traditional calculation of the PRC. The bias results from loss of causality caused by the jitter of the APs surrounding the current pulse, and gives rise to an empty triangular region in the PRC plot, which we call the "Bermuda Triangle". We introduce a method for calculating the PRC which corrects for this bias by using all spikes in the spike train as a reference, one at a time. We refer to PRCs obtained by this method as "corrected" PRCs. Note that in our study both "traditional" and "corrected" PRCs are calculated using the same experimental data: perturbation of the firing of Purkinje cells with brief square current pulses. Using the corrected method we show that the shape of the Purkinje cell PRC changes fundamentally depending on the firing rate of the neuron.

Results {#s2}
=======

A bias in the traditional method for calculating PRCs {#s2a}
-----------------------------------------------------

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 [Fig. 1A](#pcbi-1000768-g001){ref-type="fig"}). 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 ([Fig. 1B](#pcbi-1000768-g001){ref-type="fig"}, neuron firing at 180 Hz; see [Materials and Methods](#s4){ref-type="sec"}). 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 ([Fig. 1B](#pcbi-1000768-g001){ref-type="fig"}, dashed red line). Second, the intrinsic variability in the ISI [@pcbi.1000768-Husser1] 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 ([Fig. 1B](#pcbi-1000768-g001){ref-type="fig"}, 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 [@pcbi.1000768-Reyes1], [@pcbi.1000768-Tateno1]. The ISI variability in Purkinje cells [@pcbi.1000768-Husser1] 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.\
(A) Whole-cell patch-clamp recording of a spontaneously firing Purkinje cell. The injected current pulse (0.5 ms) is shown below. (B) Experimentally determined Purkinje cell PRC from the same cell as in (A) (mean firing rate 180 Hz) showing a triangular region ("Bermuda Triangle", dashed green lines) devoid of points, creating a bias at late phases of the PRC (dashed red line). Deleting the points at phases\>1 does not remove the bias (solid red line). (C) Control PRC calculated with no current pulse injection (cPRC) is flat for most of the phases, but the "Bermuda Triangle" and the negative bias at the late phases are still present. The panel below shows the phase histogram, which is inhomogeneous with smaller probabilities of sampling the ISI in the late phases. (D) When noise is introduced in a Purkinje neuron model (firing at its spontaneous firing rate of 27 Hz), the "Bermuda Triangle" and the bias at the late phases are reproduced. (E) Control PRC for the model neuron, calculated in the absence of current pulse injection. It is flat for most of the phases, but the "Bermuda triangle" and the negative deflection at the late phases are reproduced. The phase histogram below is inhomogeneous, with smaller probability of sampling the ISI in the late phases in the model.](pcbi.1000768.g001){#pcbi-1000768-g001}

To better understand the negative deflection of the PRC at late phases, a control PRC (cPRC; see [Materials and Methods](#s4){ref-type="sec"}) 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 ([Fig. 1C](#pcbi-1000768-g001){ref-type="fig"}). 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 ([Fig. 1C](#pcbi-1000768-g001){ref-type="fig"}, lower panel) indicates that late phases are sampled less frequently.

To reproduce the effect of noise, PRCs were obtained from a Purkinje cell model [@pcbi.1000768-Khaliq1] 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](#s4){ref-type="sec"}). The model PRC exhibited the same negative deflection at late phases as observed in the experimental PRC ([Fig. 1D](#pcbi-1000768-g001){ref-type="fig"}, dashed red line). As before, removing all points for which the phase exceeds 1 did not eliminate the negative deflection ([Fig. 1D](#pcbi-1000768-g001){ref-type="fig"}, solid red line). Similarly, the cPRC in the model exhibited the same negative bias and the same inhomogeneous phase distribution ([Fig. 1E](#pcbi-1000768-g001){ref-type="fig"}) 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 [Fig. 2A--L](#pcbi-1000768-g002){ref-type="fig"}), we jittered either the first or the second AP ([Fig. 2A--L](#pcbi-1000768-g002){ref-type="fig"}, black lines) with respect to a perfectly periodic cycle of firing ([Fig. 2A--L](#pcbi-1000768-g002){ref-type="fig"}, grey lines). We divided the sketches into three blocks depending on the phase of the current pulse within the cycle ([Fig. 2 A, D, G, J](#pcbi-1000768-g002){ref-type="fig"}: 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 [Fig. 2C and 2J](#pcbi-1000768-g002){ref-type="fig"}, 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 [Fig. 2C](#pcbi-1000768-g002){ref-type="fig"}: 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 [Fig. 1](#pcbi-1000768-g001){ref-type="fig"}, 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.\
Illustrations of the changes in and (arrows) introduced by noise, in the form of jittering of the spikes (black lines) surrounding the current pulse, one at a time; as compared to the deterministic case (gray lines). Noise in the PRC plots depends on the position of the pulse within the ISI: (A, D, G, J) the pulse arrives near the beginning of the ISI; (B, E, H, K) the pulse arrives near the middle of the ISI; (C, F, I, L) the pulse arrives near the end of the ISI. The color coding groups the different scenarios according to the effect of the noise on and . For example, green indicates those scenarios in which noise gave rise to an increase in and as compared to the deterministic case. (M) Actual phase plotted against the sampled phase (resulting from jittering the spikes around a mean). The color coding is the same as in (A--L). The extreme cases where noise misclassifies the stimulus in a different ISI are shown in red and purple. (N) Sampled PRC corresponding to data in (M), shown plotting sampled PRC against sampled phase. Color-coding corresponds to each of the cases in (A--L).](pcbi.1000768.g002){#pcbi-1000768-g002}

To visualize the resulting phase and PRC misclassification, we translated each of these twelve sketches onto a corresponding phase plot ([Fig. 2M,N](#pcbi-1000768-g002){ref-type="fig"}). 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 [Fig. 2M,N](#pcbi-1000768-g002){ref-type="fig"} (and subsequently [Fig. 3B--E](#pcbi-1000768-g003){ref-type="fig"}), and only for purposes of illustration. In a purely deterministic scenario, the sampled phase is linearly dependent on the actual phase ([Fig. 2M](#pcbi-1000768-g002){ref-type="fig"}, 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 [Fig. 2A](#pcbi-1000768-g002){ref-type="fig"} or [2B](#pcbi-1000768-g002){ref-type="fig"}; yellow points in [Fig. 2M](#pcbi-1000768-g002){ref-type="fig"}). For any deviations of the sampled phase from the actual phase due to noise, the points are scattered across the plot ([Fig. 2M](#pcbi-1000768-g002){ref-type="fig"}, 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 [Fig. 2N](#pcbi-1000768-g002){ref-type="fig"} (color coding as in A--L).

![A new approach for generating corrected PRCs.\
(A) Schematic illustration of a spike train, coded by color: the red spike is immediately preceding the pulse and is the one used as a reference () in the traditional method. The remaining colors correspond to the temporal sequence of spiking (green: three spikes before pulse, black: two spikes before pulse, blue: first spike after pulse, cyan: second spike after pulse). In our new approach, all spikes are used as a reference one at a time to calculate the phase of the current pulse (values under arrows) and the corresponding PRC value. The red arrow indicates the phase used to calculate the PRC using the traditional method. (B) Sampled phase vs. actual phase plot as in [Fig. 2M](#pcbi-1000768-g002){ref-type="fig"} with the color code corresponding to the colors used in (A). The two spikes preceding the pulse (red and black circles) contribute to the interval of interest \[0,1\], whereas in the traditional method only the spike immediately preceding the pulse is used as a reference (red circles). (C) Phase histograms, plotted separately for the different reference spikes in the corresponding color code. (D) Phase histogram of sampled phases using all five spikes as a reference one at a time. (E) Same sampled PRC vs. sampled phase as in [Fig. 2N](#pcbi-1000768-g002){ref-type="fig"} with the same color code as in (A). PRC~2~ and PRC~3~ are calculated by using only the spikes following the pulse as reference spikes, and correspond to the phase intervals \[−1,0\] and \[−2,−1\] respectively. The data shown in Fig. 3B--E (as well as 2M,N) assume for the purpose of illustration that the process underlying spiking is perfectly periodic and that the presence of a spike does not reset this underlying process. However, when actually applying our corrected method for calculating the PRC, the phase *is* yoked to the reference spike, i.e. phase is reset at each reference spike.](pcbi.1000768.g003){#pcbi-1000768-g003}

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 {#s2b}
------------------------------------------------------------------------

Our new method to correct for the bias in the traditional PRC and obtain a homogeneous phase histogram is illustrated in [Fig. 3A](#pcbi-1000768-g003){ref-type="fig"}. 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 [Fig. 3A](#pcbi-1000768-g003){ref-type="fig"}) are all taken into account (see also [Materials and Methods](#s4){ref-type="sec"}). In this case, the two spikes prior to the stimulation pulse (red and black in [Fig. 3A](#pcbi-1000768-g003){ref-type="fig"}) predominantly contribute to the phase interval \[0,1\] of the PRC ([Fig. 3B](#pcbi-1000768-g003){ref-type="fig"}, red and black points). The impact of the pulse on the subsequent ISI, the PRC~2~, is then determined by the two spikes following the current pulse ([Fig. 3A](#pcbi-1000768-g003){ref-type="fig"} blue and cyan spikes; also compare [@pcbi.1000768-Gutkin1]) 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 [Fig. 3B](#pcbi-1000768-g003){ref-type="fig"} (all points, in analogy to [Fig. 2M](#pcbi-1000768-g002){ref-type="fig"}). By taking only the points according to the traditional calculation of the PRC, a sharp boundary is drawn ([Fig. 3B](#pcbi-1000768-g003){ref-type="fig"}, red) resulting in an inhomogeneous distribution of sampled phases ([Fig. 3C](#pcbi-1000768-g003){ref-type="fig"}, red). In contrast, by including the second spike prior to the pulse, spike jitter effects are averaged out ([Fig. 3D](#pcbi-1000768-g003){ref-type="fig"}). The bias at the late phases of the PRC plot observed when taking points according to the traditional calculation of the PRC ([Fig. 3E](#pcbi-1000768-g003){ref-type="fig"}, red) is thereby eliminated ([Fig. 3E](#pcbi-1000768-g003){ref-type="fig"}, 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 [@pcbi.1000768-Ermentrout4]). PRCs of the Morris-Lecar model (parameters from [@pcbi.1000768-Oprisan2]), in the presence and absence of noise, were compared with the analytically derived PRC ([Fig. S1A](#pcbi.1000768.s001){ref-type="supplementary-material"}). The PRCs calculated using both the traditional and our corrected method overlap perfectly (except near and , 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 [@pcbi.1000768-Tateno2]). 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](#pcbi.1000768.s001){ref-type="supplementary-material"}). 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](#pcbi.1000768.s001){ref-type="supplementary-material"}).

Next, using the Purkinje cell model [@pcbi.1000768-Khaliq1], we compare the result of our method ([Fig. 4A](#pcbi-1000768-g004){ref-type="fig"}, black) to the PRC obtained with the traditional method ([Fig. 4A](#pcbi-1000768-g004){ref-type="fig"}, red) and the deterministic PRC without noise ([Fig. 4A](#pcbi-1000768-g004){ref-type="fig"}, 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 ([Fig. 4A](#pcbi-1000768-g004){ref-type="fig"}, dashed red line). In contrast, our corrected method performs as well as with low CV ([Fig. 4A](#pcbi-1000768-g004){ref-type="fig"}, 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; [Fig. 4B](#pcbi-1000768-g004){ref-type="fig"}). As the CV increases, the PRC error shows larger increases using the traditional ([Fig. 4B](#pcbi-1000768-g004){ref-type="fig"}, red line) compared to our corrected method ([Fig. 4B](#pcbi-1000768-g004){ref-type="fig"}, black line).

![Validation of the corrected method for obtaining PRCs.\
(A) PRCs of a model Purkinje cell, calculated using the corrected method in the stochastic model (black), in the deterministic model (green) and calculated using the traditional method in the stochastic model (red). Traces are shown at two different noise levels (increased CV indicated by dashed red and dashed black lines respectively). (B) PRC error, expressed as the integral of the differences between PRCs and their deterministic counterparts. The corrected method performs better than the traditional method in all cases, particularly at high CV.](pcbi.1000768.g004){#pcbi-1000768-g004}

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 {#s2c}
------------------------------------------------------

Spontaneous firing frequencies of Purkinje cells range from 10--180 Hz both *in vitro* [@pcbi.1000768-Husser1], [@pcbi.1000768-Williams1], [@pcbi.1000768-Watt1], [@pcbi.1000768-Raman1] and *in vivo* [@pcbi.1000768-Loewenstein1], [@pcbi.1000768-Shin1]. 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 [Fig. 5A](#pcbi-1000768-g005){ref-type="fig"} (the same example of a rapidly firing (180 Hz) Purkinje neuron as in [Fig. 1B](#pcbi-1000768-g001){ref-type="fig"}). 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 ([Fig. 5A](#pcbi-1000768-g005){ref-type="fig"}, 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 PRC~2--5~ ([Fig. 5B](#pcbi-1000768-g005){ref-type="fig"}; see [Materials and Methods](#s4){ref-type="sec"}). PRC~2~ is negative, suggesting that the subsequent ISI is lengthened relative to the mean. A PRC~2~ with opposite sign to the PRC has been previously reported [@pcbi.1000768-Gutkin1] and it is believed to be due to a compensatory effect on the current ISI length. Indeed, as seen from PRC~2--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.\
(A) PRC of a rapidly firing cell. Using the corrected method the "Bermuda triangle" (shown in green dashed lines) and the bias at the late phases are eliminated. The brief current pulses cause an ISI shortening as reflected by the positive PRC values. The panel below shows the phase histogram, which is homogeneous, suggesting an equally sampled ISI. (B) The effect of the brief current pulse is persistent in the subsequent ISIs (PRC~2--5~) and dies out after four ISIs. PRC~2--5~ are negative suggesting ISI lengthening. (C) PRC of a slowly firing cell. The PRC is square/flat; the brief current pulse causes the same effect independent of its position within the ISI. Again, the phase histogram, below, is homogeneous revealing an equally sampled ISI. (D) Similarly to (B), the brief current pulse lengthens subsequent ISIs: PRC~2--5~ are negative.](pcbi.1000768.g005){#pcbi-1000768-g005}

In comparison, an example of a PRC of a slowly firing (30 Hz) Purkinje neuron is shown in [Fig. 5C](#pcbi-1000768-g005){ref-type="fig"}. 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 ([Fig. 5C](#pcbi-1000768-g005){ref-type="fig"}, lower panel). In order to study the effects of the brief pulse on the subsequent intervals we calculated the PRC~2--5~ ([Fig. 5D](#pcbi-1000768-g005){ref-type="fig"}; see [Materials and Methods](#s4){ref-type="sec"}). 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 ([Fig. 6A](#pcbi-1000768-g006){ref-type="fig"}) to phase-dependent for higher frequencies ([Fig. 6B](#pcbi-1000768-g006){ref-type="fig"}). The switch occurred at a frequency of approximately 50 Hz. The average PRC of all neurons firing at low rates ([Fig. 6A](#pcbi-1000768-g006){ref-type="fig"}, 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 ([Fig. 6B](#pcbi-1000768-g006){ref-type="fig"}, thick line) exhibited a sharp peak. It is useful to compare these average PRCs ([Fig. S2](#pcbi.1000768.s002){ref-type="supplementary-material"}, thick black and red lines) with the biased ones obtained with the traditional method ([Fig. S2](#pcbi.1000768.s002){ref-type="supplementary-material"}, 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 ([Fig. 6C](#pcbi-1000768-g006){ref-type="fig"}; see [Materials and Methods](#s4){ref-type="sec"}). 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 ([Fig. 6D](#pcbi-1000768-g006){ref-type="fig"}). 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 [Fig. 6B](#pcbi-1000768-g006){ref-type="fig"} and [Fig. 6D](#pcbi-1000768-g006){ref-type="fig"}, 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.\
(A) Individual PRCs (thin dashed lines) and the average PRC (thick line) of slowly firing cells (frequency\<50 Hz, n = 10, black). The average PRC is independent of the phase of the brief current pulse within the ISI. (B) Individual PRCs (thin lines) and the average PRC (thick line) of rapidly firing cells (frequency\>50 Hz, n = 6, red). The PRC is phase-dependent with a peak at late phases. The thin solid line is the PRC of the cell with the lowest firing frequency (55 Hz) in this group (arrows in (B--D) and shows both a peak at late phases and a plateau at early phases. (C) Peak-to-baseline ratio (see [Materials and Methods](#s4){ref-type="sec"}) plotted against firing frequency: as the frequency increases the PRC switches from phase-independent to phase-dependent. Points shown in two colors correspond to PRCs from the same cell (see [Figure 7](#pcbi-1000768-g007){ref-type="fig"} and text for more details). (D) PRCs of the two groups of cells (black: low firing rate, same data as in (A); red: high firing rate, same data as in (B) plotted on the same axes with both the phase and the phase shift of the PRC in units of time, and phases aligned with respect to the second AP in the ISI. The peaks seen at high firing rates coincide at a time of approximately --3 ms from the second spike.](pcbi.1000768.g006){#pcbi-1000768-g006}

The firing rate of a Purkinje cell changes depending on modulation of its inputs [@pcbi.1000768-Thach1]--[@pcbi.1000768-Edgley1]. For example, during locomotion in cats the firing frequencies of Purkinje cells can increase from an average of about 40 Hz [@pcbi.1000768-Armstrong1] to more than 100 Hz [@pcbi.1000768-Edgley1]. 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; [Fig. 6C](#pcbi-1000768-g006){ref-type="fig"}, 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 ([Fig. 7](#pcbi-1000768-g007){ref-type="fig"}, color coding as in [Fig. 6C](#pcbi-1000768-g006){ref-type="fig"}). When Purkinje cell spiking was changed from slow (33 Hz) to fast (104 Hz), the originally square PRC ([Fig. 7A](#pcbi-1000768-g007){ref-type="fig"}), exhibited a sharp peak ([Fig. 7B](#pcbi-1000768-g007){ref-type="fig"}). 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 ([Fig. 7C](#pcbi-1000768-g007){ref-type="fig"}). Conversely, another neuron initially firing at a high rate (71 Hz) exhibited a peaked PRC ([Fig. 7D](#pcbi-1000768-g007){ref-type="fig"}), which was switched to a square shape by reducing its firing rate to 26 Hz via injection of hyperpolarizing current ([Fig. 7E](#pcbi-1000768-g007){ref-type="fig"}). When the neuron was then allowed to fire at its intrinsic firing rate (84 Hz) the sharp peak in the PRC reappeared ([Fig. 7F](#pcbi-1000768-g007){ref-type="fig"}). 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.\
(A) A slowly firing cell (33 Hz) displaying a phase-independent PRC. (B) Injection of depolarizing current into the cell in (A) increases its firing rate (104 Hz) and the PRC switches to being phase-dependent. (C) When the depolarizing current is removed, the cell relaxes back to its intrinsic spontaneous firing rate (40 Hz) and its PRC switches back to being phase-independent. (D) A rapidly firing cell (71 Hz) displaying a phase-dependent PRC. (E) When hyperpolarizing current is injected into the cell to decrease its firing rate (26 Hz), the PRC switches to being phase-independent. (F) Removal of the hyperpolarizing current returns the cell to its intrinsic spontaneous firing rate (84 Hz) and its PRC switches back to being phase-dependent.](pcbi.1000768.g007){#pcbi-1000768-g007}

Discussion {#s3}
==========

We have shown that the traditional method for calculating PRCs results in a bias, particularly in neurons exhibiting high ISI variability. We developed a corrected method for calculating PRCs which removes most of this bias. Our method can be directly applied to noisy experimental data. We used this corrected approach to measure for the first time the PRCs of Purkinje cells at various firing rates. At high firing rates, Purkinje cell PRCs were phase-dependent; however, a phase-independent PRC was observed at lower firing rates. This suggests that Purkinje cells can behave as perfect integrators at low firing rates, which has important consequences for our view of the integrative properties of these neurons.

A new approach for determining PRCs {#s3a}
-----------------------------------

We have determined Purkinje cell PRCs by injecting brief current pulses and measuring the phase change in the subsequent neuronal spiking. Since at the typical spontaneous firing rates of Purkinje cells these phase changes were small compared to the spike jitter during spontaneous spiking [@pcbi.1000768-Husser1], many trials were required. This revealed a general bias of the traditional method at late phases of the PRC in the presence of noise ([Fig. 1](#pcbi-1000768-g001){ref-type="fig"}). We characterized the effect in a model with and without noise, and showed that the bias is related to inhomogeneous phase histograms caused by interspike interval jitter ([Fig. 1](#pcbi-1000768-g001){ref-type="fig"} and [Fig. 3](#pcbi-1000768-g003){ref-type="fig"}). To correct for this, we developed a new method, which recovers periodicity in the spike jitter due to noise ([Fig. 3](#pcbi-1000768-g003){ref-type="fig"}). We showed that this method homogenizes the phase sampling in the experimental data and removes most of the bias observed in the PRCs calculated using the traditional method ([Fig. 4A](#pcbi-1000768-g004){ref-type="fig"}). Our corrected approach can be directly applied to existing experimental data in order to measure PRCs under low signal-to-noise conditions. It should be applicable to a wide range of cell types, as neuronal noise and the resulting ISI variability are not restricted to Purkinje cells [@pcbi.1000768-Schneidman1].

The use of indirect methods to obtain PRCs, for example from the spike triggered average [@pcbi.1000768-Ermentrout3] or the poststimulus time histogram (PSTH) [@pcbi.1000768-Gutkin1] are possible alternatives to the traditional method. Here we have applied a correction to the traditional method, which resulted in reliable PRC measurements in Purkinje cells. Further alternative methods for calculating PRCs exist. For example, dynamic clamp was previously used to study hippocampal spike-timing-dependent plasticity in relation to PRCs [@pcbi.1000768-Lengyel1]. In this special case, underlying subthreshold oscillations provide phase locking. Such a method is only applicable if phase information is accessible to the experimenter, independent of spiking. PRCs can also be calculated using Bayesian statistics [@pcbi.1000768-Ota1], or by injecting trains of rectangular current pulses [@pcbi.1000768-Galn3] and noisy inputs [@pcbi.1000768-Izhikevich1]. These methods result in periodic PRCs, but only because periodicity is imposed as part of the optimization (fitting) techniques employed. In conclusion, our method can be applied to noisy experimental data to calculate PRCs while avoiding possible bias or overfitting problems present in some of the currently available methods. A wide, comparative study will be required in the future to find out which methods for calculating the PRC yield the best results under different conditions.

Purkinje cell dynamics depend on firing rate {#s3b}
--------------------------------------------

Purkinje cells fire spontaneously and modulate their firing in response to synaptic input. The spontaneous firing rate of Purkinje cells varies from 10 to 180 Hz, but firing frequency can also be increased by the ∼150,000 parallel fiber synaptic inputs [@pcbi.1000768-Napper1] or decreased by molecular layer interneurons during the execution of motor tasks such as smooth-pursuit eye movements [@pcbi.1000768-Lisberger1], maintenance of posture [@pcbi.1000768-Thach2] and locomotion [@pcbi.1000768-Edgley1], [@pcbi.1000768-Udo1]. For example, the rate of Purkinje cell firing can exhibit a consistent temporal relationship with wrist movement [@pcbi.1000768-Thach1] or be monotonically related to eye velocity during smooth-pursuit eye movements [@pcbi.1000768-Lisberger1].

How is the integration of single inputs affected by the firing rate of the Purkinje cell? We have addressed this question by measuring the PRC at different firing rates. Using our new approach, we determine experimentally the PRCs of cerebellar Purkinje cells and show that their shape changes significantly depending on the firing rate (compare [Fig. 5A](#pcbi-1000768-g005){ref-type="fig"} and [Fig. 5C](#pcbi-1000768-g005){ref-type="fig"}). At high firing frequencies (\>50 Hz) Purkinje cell PRCs are monophasic ([Fig. 6B](#pcbi-1000768-g006){ref-type="fig"}). However, at low firing rates (\<50 Hz), Purkinje cell PRCs become phase-independent ([Fig. 6A](#pcbi-1000768-g006){ref-type="fig"}). To the best of our knowledge, this is the first study to report a phase-independent PRC in a mammalian neuron.

It was previously reported in a spike-frequency adaptation model of cortical neurons that an increase in firing frequency causes a shift of the PRC peak from rightward skew to the centre with a decrease in amplitude [@pcbi.1000768-Gutkin1], implying that the integrative properties of this model neuron change depending on the firing rate. Specifically, it was suggested that the model cell acts like a coincidence detector at low firing rates and more of an integrator at higher firing rates [@pcbi.1000768-Gutkin1]. Purkinje cells appear to show the opposite behaviour, acting as perfect integrators at low firing rates.

Functional implications {#s3c}
-----------------------

The shape of the PRC is thought to be linked to the type of excitability of the neuron. Neurons with type I excitability, whose f-I curves are continuous, are thought to display purely positive PRCs while neurons with type II excitability, characterized by a discontinuity in the f-I curve at the onset of firing, exhibit biphasic PRCs [@pcbi.1000768-Izhikevich1], [@pcbi.1000768-Izhikevich2], [@pcbi.1000768-Galn1]. While biphasic PRCs intuitively result in resonator behavior, neurons with purely positive PRCs act as integrators of synaptic input [@pcbi.1000768-Izhikevich1], [@pcbi.1000768-Izhikevich2], [@pcbi.1000768-Galn1], [@pcbi.1000768-Stiefel1]. Although Purkinje cells exhibit type II excitability [@pcbi.1000768-Williams1], [@pcbi.1000768-Llins1], [@pcbi.1000768-Fernandez1], their PRCs are positive at all firing rates, implying that they are integrators rather than resonators. These findings suggest that the type of excitability of a neuron is not strictly correlated with the PRC shape. Similarly, Tateno and Robinson [@pcbi.1000768-Tateno1] showed that low-threshold spiking, fast spiking and non-pyramidal regular spiking interneurons can exhibit both purely positive and biphasic PRCs which do not always strictly correspond to the type of excitability of the neuron.

The shape of the PRC has functional implications for the integration of synaptic inputs. At high firing rates, Purkinje cells are most sensitive to inputs during the last 3 ms of their firing cycle ([Fig. 6D](#pcbi-1000768-g006){ref-type="fig"}), imposing a strict relationship between the timing of the input and the timing of spike output, with direct consequences for network dynamics. It has been shown theoretically that oscillators which are described by type I PRCs and are coupled by excitatory synapses tend not to synchronize [@pcbi.1000768-Ermentrout1]. However, the opposite is true for inhibitory coupling between oscillators [@pcbi.1000768-Ermentrout1], [@pcbi.1000768-VanVreeswijk1], such as coupled Purkinje cells. Indeed, theoretical and experimental evidence indicates that Purkinje cells tend to synchronize via inhibitory inputs [@pcbi.1000768-Watt1], [@pcbi.1000768-Maex1], [@pcbi.1000768-deSolages1].

As the firing rate of Purkinje cells decreases, and the levels of synaptic and intrinsic conductances and currents are modified, the PRC switches from monophasic to phase-independent ([Fig. 6C](#pcbi-1000768-g006){ref-type="fig"}). The phase-independent PRCs at low firing rates suggest that Purkinje cells integrate their synaptic inputs independently of their timing within the interspike interval ([Fig. 6A](#pcbi-1000768-g006){ref-type="fig"}). Our results therefore support the idea that at low firing rates, Purkinje cells cannot read out the timing of their inputs, which would exclude the use of a temporal code. Instead, in this regime they are well suited for rate coding.

What are the biophysical mechanisms responsible for the switch in PRC behaviour at different firing rates? To generate an entirely flat PRC would require a neuron to effectively completely compensate for its leak conductance. This is illustrated by the example of the PRC of a simple leaky integrate-and-fire neuron in which the leak conductance was eliminated ([Fig. S1B and C](#pcbi.1000768.s001){ref-type="supplementary-material"}). However, this absence of leak is unlikely to occur in real Purkinje cells, and the biophysical implementation remains unknown. PRCs qualitatively similar to those observed in our experiments at high firing rates can be generated by the Purkinje cell model of Khaliq and colleagues [@pcbi.1000768-Khaliq1] ([Fig. S3A](#pcbi.1000768.s003){ref-type="supplementary-material"}). However, when the firing rate is lowered in the model, no qualitative switch in the shape of the PRC can be observed. A hint to the mechanisms underlying the switch in the experiment is provided by using the model of Akemann and Knöpfel [@pcbi.1000768-Akemann1] (a further development of the Khaliq et al. model): at low firing rates a 'shoulder' appears in the early phases of the PRC ([Fig. S3B](#pcbi.1000768.s003){ref-type="supplementary-material"}). However, none of these models fully capture the experimentally determined switch in Purkinje cells, perhaps reflecting the fact that both of these models represent dissociated Purkinje cells. Thus, our experimental results could aid the refinement of existing models in order to capture the full dynamic behaviour of Purkinje cells.

In conclusion, our experimental findings indicate that Purkinje cells display different dynamic behavior depending on their firing rate. At high firing rates these neurons act as coincidence detectors of synaptic inputs, with maximal sensitivity at the late phases of the interspike interval. In contrast, at low firing rates Purkinje cells are not suited for precise coincidence detection, but instead appear to perfectly integrate their inputs independently of their position within the interspike interval. Thus, at high firing rates Purkinje cells can transmit information via a temporal code whereas at low firing rates they are well-suited for rate coding.

Materials and Methods {#s4}
=====================

Ethics statement {#s4a}
----------------

All procedures were approved by the U.K. Home Office.

Acute slice preparation and electrophysiological recordings {#s4b}
-----------------------------------------------------------

Twelve- to fifteen-day-old *L7-tau-gfp* mice [@pcbi.1000768-Sekirnjak1] were anaesthetised using isoflurane, decapitated and their brains were transferred to ice-cold low Ca^2+^ artificial cerebrospinal fluid (ACSF) containing (in mM): 125 NaCl, 26 NaHCO~3~, 25 glucose, 2.5 KCl, 26 NaH~2~PO~4~, 0.5 CaCl~2~ and 3 MgCl~2~, saturated with carbogen (95% oxygen and 5% carbon dioxide gas). 230-µm-thick sagittal brain slices from the cerebellar vermis were cut on a VT1200S microtome (Leica Microsystems) and were transferred to normal ACSF containing the following (in mM): 125 NaCl, 26 NaHCO~3~, 25 glucose, 2.5 KCl, 26 NaH~2~PO~4~, 2 CaCl~2~ and 1 MgCl~2~, again bubbled with carbogen. The slices were incubated for 30--40 minutes at 37°C and were then allowed to cool to room temperature. Thick-walled, filamented, borosilicate glass electrodes (Harvard Apparatus Ltd.) were pulled to a tip resistance of 4--5 MΩ (PC-10 microelectrode puller, Narishige). Cells were visually identified with the aid of an upright infrared differential interference contrast (IR-DIC) microscope (Axioskop, Carl Zeiss) and a video camera (C2400-07, Hamamatsu). Purkinje cell somatic whole-cell patch-clamp recordings were obtained using an internal solution containing the following (in mM): 130 methanesulfonic acid, 10 HEPES, 7 KCl, 0.05 EGTA, 2 Na~2~ATP, 2 MgATP, 0.5 Na~2~GTP and 0.4% biocytin, pH-adjusted to 7.3 with KOH. All recordings were performed at 34.5±1°C in the presence of carbogen-bubbled ACSF supplemented with GABA~A~ receptor blocker SR95331 (10 µM). Recordings were made with an Axoclamp 2B amplifier (Axon Instruments) and were filtered at 3 kHz and sampled at 50 kHz using an ITC-18 DAC board (Instrutech) and Axograph 4.9 (Axon Instruments). Series resistance and pipette capacitance were carefully monitored and compensated throughout the experiment. Methanesulfonic acid was obtained from Fluka, and other chemicals from Sigma-Aldrich and BDH Chemicals.

Phase response curves {#s4c}
---------------------

Data were analyzed with MATLAB (The MathWorks). To determine how spike timing during spontaneous firing is shifted by a brief perturbation, we injected rectangular current pulses of 0.5 ms duration and 50 pA amplitude, after a baseline of 150 ms (50 ms) of spontaneous firing in subsequent trials of 350 ms (100 ms) for a slowly (rapidly) firing cell. A control PRC (cPRC) was calculated using the unperturbed part of the voltage traces and assuming a current pulse injection (0 pA amplitude) after 50 ms (25 ms) of spontaneous firing in subsequent trials of 350 ms (100 ms) for a slowly (rapidly) firing cell. The cPRC should be zero throughout all phases.

The dynamics of a neuronal oscillator can be reduced to a single variable: the phase . is calculated by dividing the time from the previous spike by the period of the oscillation; it increases linearly from 0 to 1 between two spikes. Depending on the phase of the stimulus, a change in phase, , of subsequent spiking will occur.

Traditional method: A brief current pulse is injected at a random time. The spikes before and after it are identified. is calculated by the difference between the unperturbed and the perturbed [@pcbi.1000768-Reyes1]--[@pcbi.1000768-Canavier1]. When the unperturbed is defined as the mean ISI (), a point on the PRC plot becomes:where denotes the ISI which contains the brief current pulse and is the PRC point calculated in reference to the spike just prior to the stimulus. is the time between the pulse and the preceding spike . The resulting curve is a plot of against . The curve is positive (negative) when the injected current advances (delays) the next spike. In the experimental and model (with noise) PRCs, we refer to raw data as the estimated measurements ('points') on the PRC plot. A moving average was calculated with a Gaussian kernel over the raw data.

Corrected method: A major problem with the traditional method is the loss of periodicity of the sampling reference ([Fig. 3B](#pcbi-1000768-g003){ref-type="fig"}), which results in an inhomogeneous sampling of phases in the presence of spike jittering. In order to restore periodicity, points unaffected by the stimulation pulse can be added to the ensemble of PRC points, which allows the spiking jitter to average out properly. These points can be obtained from the same data by adding PRC values when the preceding ISI is taken into account:When preceding and subsequent ISIs are taken into account as in:andperiodicity in the spiking jitter is restored, phases are sampled homogeneously and the cPRC becomes flat. In the resulting plot, the phase interval ranges from to and the PRC component affecting directly the interval corresponds to all points in the phase interval \[0,1\], termed PRC~1~. Successive PRC~2--5~, correspond to phase intervals \[−1,0\], \[−2,−1\], \[−3,−2\] and \[−4,−3\], respectively and indicate how are affected by the pulse.

*Peak-to-baseline ratio:* In order to distinguish the phase-independent PRCs from the phase-dependent ones, PRCs were classified according to the peak- to-baseline ratio. Inspired by Tateno and Robinson [@pcbi.1000768-Tateno1], local extrema at the two halves of the PRC (i.e. for and ) were calculated and were denoted as early () and late () respectively. The peak-to-baseline ratio is then defined as:

Simulations {#s4d}
-----------

Simulations were performed in NEURON [@pcbi.1000768-Carnevale1] using a model of Purkinje cells consisting of a single compartment [@pcbi.1000768-Khaliq1], [@pcbi.1000768-Akemann1]. The model includes seven voltage-gated conductances (a resurgent Na^+^ current, fast and slow K^+^ currents, P-type Ca^2+^ current, Ca^2+^-activated K^+^ current and the hyperpolarization-activated current I~h~) and one voltage-independent conductance (I~leak~), based on voltage clamp measurements from Purkinje cells [@pcbi.1000768-Khaliq1]. The membrane surface area of the neuron was modified (×13) to reproduce input resistance values close to those observed in Purkinje cells (80 MΩ). In order to mimic the noise observed in Purkinje cells, noisy current input drawn from a normal distribution with (mean) and (standard deviation) was injected at each time step of the simulation (every 25 µs) into the soma. The noise injection resulted in a coefficient of variation of ISIs of 0.05, which is comparable to the values measured in real Purkinje cells in the experiments presented here (see also [@pcbi.1000768-Husser1]). Current pulses of 0.5 ms duration and 250 pA amplitude were injected after 2500 ms, at a time at which spike jitter had randomized spiking phase. Data shown is taken from more than 15000 trials.

Additional neuron models were used in the supplementary parts of the manuscript. For [Fig. S1](#pcbi.1000768.s001){ref-type="supplementary-material"}, the Morris-Lecar model was directly implemented using parameters from [@pcbi.1000768-Oprisan2]. The adjoint was calculated using XPPAUT [@pcbi.1000768-Ermentrout4] and the PRCs with noise were integrated in MATLAB (The MathWorks). The parameters for the leaky integrate-and-fire model were: a membrane time constant of , a reset potential of , a threshold potential of , a membrane resistance of , and a steady driving current of (to result in 50 Hz firing) and was simulated at time steps of . For the non-leaky integrate-and-fire model the time constant was set to infinity and , otherwise the same parameters were used. An alternative model for Purkinje cell firing was used for [Fig. S3](#pcbi.1000768.s003){ref-type="supplementary-material"} which also includes a resurgent Na^+^-current and modified voltage-gated K^+^-conductances [@pcbi.1000768-Akemann1]. In this model, current pulses of 0.5 ms duration were injected at amplitudes of 10 pA in the low firing rate (33 Hz) case and 60 pA in the high firing rate (111 Hz) case. Simulation results were analysed in the same way as the experimental data.

Supporting Information {#s5}
======================

###### 

Validation of the corrected method. (A) Comparison between the traditional method (red line) and the corrected method (black line) to obtain PRCs and their numerical (green line) and analytical (blue dashed line) no-noise pendants (using the example of the Morris-Lecar model for which the analytical PRC can be calculated by the adjoint method). Curves have been rescaled to their maxima to aid comparison. (B) and (C) Same procedure for a non-leaky and a leaky (5 ms time constant) integrate-and-fire (I&F) model. In all cases the corrected method performs better than the traditional method. The use of the corrected method is particularly important in (B) which corresponds best to the case observed in the experimental data from Purkinje cells at low firing rates. Note that here we find phase-independent and phase-dependent PRCs in a simple I&F neuron. A square PRC can be obtained by eliminating the leak of the I&F neuron, compare (B) and (C). Purkinje cells therefore act as perfect non-leaky integrators at low frequencies (compare [Fig. 6A](#pcbi-1000768-g006){ref-type="fig"}). A phase-dependent PRC can be obtained by adding a leak to the I&F neuron (right). This suggests that Purkinje cells act like leaky integrators at high frequencies (compare [Fig. 6B](#pcbi-1000768-g006){ref-type="fig"}).

(0.52 MB TIF)

###### 

Click here for additional data file.

###### 

Comparison of the corrected and traditional methods for obtaining PRCs. (A) and (B) PRCs of individual Purkinje cells firing at low rates (A, black) and high rates (B, red) and population averages (thick lines), exactly as in [Figure 6AB](#pcbi-1000768-g006){ref-type="fig"} of the main manuscript. (C) and (D) PRCs of the same cells as in (A) and (B) but obtained using the traditional method. The population averages obtained with the traditional method (thick green lines) are qualitatively different from the population averages obtained using the corrected method (thick black and red lines). The bias is such that the conclusions obtained in this study would not have been possible without developing the new method.

(0.32 MB TIF)

###### 

Click here for additional data file.

###### 

PRCs in different model neurons. (A) PRCs obtained with the model of Khaliq et al. [@pcbi.1000768-Khaliq1] at a low firing rate (left) and a high firing rate (right) both exhibit monophasic behavior as seen in the experiment at high firing rates. (B) Frequency dependence of PRCs in an alternative Purkinje cell model from Akemann and Knöpfel [@pcbi.1000768-Akemann1]. A hint of a 'shoulder' at low firing rates indicates the trend of a frequency dependent switch in the model. However, the PRCs at low firing rates (left) are still not flat.

(0.21 MB TIF)

###### 

Click here for additional data file.

We thank Jenny Davie, Alanna Watt, Beverley Clark and Alexandre Mathy for help with experiments, and Mauricio Barahona, Boris Gutkin, Martha Nari Havenith and Simon Schultz for helpful discussions.

The authors have declared that no competing interests exist.

This work was supported by the BBSRC, EPSRC, the Humboldt Foundation, Gatsby Charitable Foundation and The Wellcome Trust. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

[^1]: Conceived and designed the experiments: EP HC AR MH. Performed the experiments: EP. Analyzed the data: EP HC. Contributed reagents/materials/analysis tools: EP HC. Wrote the paper: EP HC AR MH. Performed the modeling: EP HC AR.
