A Comparison of Methods to Determine Neuronal Phase-Response Curves

The phase-response curve (PRC) is an important tool to determine the excitability type of single neurons which reveals consequences for their synchronizing properties. We review five methods to compute the PRC from both model data and experimental data and compare the numerically obtained results from each method. The main difference between the methods lies in the reliability which is influenced by the fluctuations in the spiking data and the number of spikes available for analysis. We discuss the significance of our results and provide guidelines to choose the best method based on the available data.


INTRODUCTION
The phase-response curve (PRC) of a regularly fi ring neuron quantifi es the shift in the next spike time as a function of the timing of a small perturbation delivered to that neuron. The PRC is an important measure for several reasons. First, the ability of neurons to synchronize in excitatory coupled pairs, chains or networks can be predicted from the PRC: type-I PRCs (purely positive, all excitatory perturbations lead to an acceleration of spiking) do not allow synchronization while type-II PRCs (biphasic, acceleration or delay of spiking depending on the phase of the perturbation) allow synchronization with excitatory connections and short delays (Hansel et al., 1995). Furthermore, the PRC is informative about the type of bifurcation leading from rest to spiking (Izhikevich, 2007), thus constraining quantitative models of the neuron under investigation. Also, the PRC is correlated with the type of excitability of a neuron (Hodgkin and Huxley, 1952;Marella and Ermentrout, 2008).
More precisely, a regular fi ring neuron can be seen as a stable oscillator with period T and only the phase φ to describe its state. T results from the characteristic angular velocity ω of the oscillator, thus dφ/dt = ω. In the absence of inputs a regularly fi ring neuron fi res exactly when φ = kT (with k an integer and T corresponding the average ISI, ISI ). Now suppose an input to that neuron with a small amplitude at phase φ, Π(φ). Then, the infl uence of this perturbation on the next spike time is described by dφ/dt = ω + Π(φ)Z(φ) where Z(φ) is the PRC. In other words, the time required to reach the next spike deviates from T according to the perturbation and the PRC. Since the exact spike times, perturbations Π(φ) and ω (i.e.,2π/T) are known, we can estimate the PRC Z(φ) from these data.
Several methods have been proposed to compute the PRC from experimental or modeled data. For basic neuronal models, the PRC can be directly computed from the underlying differential equations by the adjoint method (Ermentrout, 1996), but for all To alleviate this problem, novel methods have been proposed that use predictions of how spike times will be altered by incoming pulses, and, methods that use continuous fl uctuation signals to obtain a more robust PRC measurement based on less spikes. The fi ve reviewed methods consist of one variation of the direct method (Galan's method), two methods that use spike-time predictions to reconstruct the PRC [the modifi ed-Izhikevich method and the standardized error prediction (STEP) method], and two methods that derive the PRC from the incoming continuous fl uctuating signal [the spike-triggered average (STA) and weighted spike-triggered average (WSTA) method]. The different methods are outlined below and illustrated in Figure 1.

Galan's method
Galan's method (Galán et al., 2005) uses pulses as perturbations (see the top panel of Figure 1), and fits the PRC to the spike time shifts as a function of the phase of the perturbation. This is one of the methods which is an extension of the direct method for noisy data (Reyes and Fetz, 1993;Tsubo et al., 2007;Stiefel et al., 2009). The PRC Z(φ) is substituted by a truncated Fourier series, i.e., Z n n cos and the parameters describing the curve (a and b) are optimized using the Euclidean distance between the data points and the curve as an error signal. The resulting curve is the best approximation of the PRC (only constrained by the length of the expansion of the Fourier series).

Modifi ed-Izhikevich method
Izhikevich proposed an inverse solution to compute the PRC in Izhikevich (2007). The method relies on predicting the next spike time and minimizing the error between the predicted spike time and the true next spike time. The prediction is based on the sum of the phase shifts that a small perturbation (part of a continuous fl uctuation) would cause. However, the proposed method does not converge to a correct solution after a reasonable number of fi tting rounds (e.g., in about 2-3 h of computation while other methods converge within minutes). Therefore, we modifi ed the method and used it with less complex perturbation data. In the modifi ed-Izhikevich method, one pulse x(t) is injected per phase and the next spike is predicted according to a candidate PRC, i.e., s x z t c is the candidate PRC. Then, the candidate PRC is optimized to match the spiking data by computing an error signal proportional to the difference between the predicted next spike time and the actual next spike time, e.g., Err φ = − ( ) + + s n n s 1 1 .

Spike-triggered average method
In the STA method (Ermentrout et al., 2007) the STA is computed from continuous low amplitude current fl uctuations (see Figure 1). Then, this STA is numerically integrated to produce the PRC. The connection between the integral of the STA and the PRC is proven for regularly fi ring neurons and small perturbations in Ermentrout et al. (2007) holds. In contrast to both Galan's method and the modifi ed-Izhikevich method, only a single pass over the complete noise signal and the voltage trace is required because there is no optimization step. From this single pass, the STA is computed and subsequently integrated.

Weighted spike-triggered average method
The WSTA method devised in Ota et al. (2009a) is an extension of the STA method and also integrates the continuous low amplitude current fl uctuations to derive the PRC. However, in the WSTA method, the fl uctuations in between the different spike times are normalized to the average ISI ISI ( ) of all spikes in s ( s , ISI τ with instantaneous ISI τ i ). Then, as Ota et al. (2009a) prove, with an appropriate weighting function the weighted sum of these normalized stretches of current fl uctuations constitutes the PRC (for regularly fi ring neurons). The weighing function is α τ τ = − ( )

Standardized error prediction method
The STEP method (Torben-Nielsen et al., 2010) is an extension of the modifi ed-Izhikevich method to work with continuous fl uctuation data (in a different way than originally proposed by Izhikevich). In this method, instead of using a prediction error averaged over all ISIs and all phases, the temporal information in the error is preserved by binning the errors of all ISIs independently per phase. The fl uctuations are binned equidistantly on ISI (i.e., normalized) and are treated as if they were independent, i.e., for each phase bin of each ISI, the predicted next spike time is computed and the mismatch with the true next spike time is used to optimize the parameters of a curve representing the PRC. More precisely, all inter-spike intervals are normalized to [0,2π] and discretized into N bins. For each bin, an independent prediction is made about the next spike time: , where φ corresponds to phase of bin j . Subsequently, one obtains a two-dimensional array with the dimensions given by N bins and M spikes. Finally, a leastsquares fi tting algorithm is used to minimize the 2-D prediction error array and obtain a PRC that predicts the recorded spike time shifts (Torben-Nielsen et al., 2010). Table 1 summarizes the fi ve implemented methods (plus the direct method) and how they relate to each other. In addition to the reviewed methods, there are several other published methods which we omitted because they proved impractical (with respect to experimental demands), e.g., the MAP-estimation algorithm (Ota et al., 2009b) and the post-stimulus time histogram method (Gutkin et al., 2005).

DATA SETS
We tested the fi ve methods (and the direct method) with three different data sets whenever possible. The fi rst two data sets contain model data while the third set contains experimental data. We used the single-compartmental model as developed by Golomb and Amitai (1997) and modifi ed by Stiefel et al. (2009) to generate the data. This model uses a Hodgkin-Huxley-type formalism to model neural spiking behavior:   Table 1.
where V is the membrane potential, g x the maximum conductance for ion x and E x the reversal potential for ion x. The parameter values can be found in Golomb and Amitai (1997) and Stiefel et al. (2009). By turning the adaptation current on or off, this model switches between type-II or type-I excitability (Stiefel et al., 2009), respectively.
The fi rst data set contains noise-free model data in which a single compartmental model neuron (see below) is perturbed at different phases. This set contains 128 pulses evenly spaced over [0,2π]. The second data set contains modeled data from the same single-compartmental model but with an additionally injected fl uctuating current. The fl uctuations are generated through a stationary Ornstein-Uhlenbeck process around a given mean value and parameterized by the reversion rate (g = 0.1) and four different volatility levels (D = 1e −4 ,5e −4 ,1e −5 ,5e −5 ). The advantage of the noisy modeled data is that the excitability type is known with certainty because small perturbations do not change the PRC type (Izhikevich, 2007). The injected fl uctuating current and the resulting spike trains are illustrated in Figure 2. The different noise levels result in four groups of data each containing approximately 950 spikes. The noise and the resulting spike trains are illustrated in Figure 2. The last data set contains experimental data recorded from a layer 3/4 pyramidal cell of the mouse visual cortex with the whole cell patch-clamp technique in vitro. Standard patch-clamp techniques as in Stiefel et al. (2008) were used. Membrane potential voltage data and the injected fl uctuations were digitized at 40 kHz, and two levels of current (µ = 50 and 100 pA) were injected as fl uctuations. The fl uctuations consisted of white-noise low-pass fi ltered at 200 Hz. In both the model data and experimental data, the fl uctuations are on top of a step current (I s ) which is required to get the model/cell into a regime of regular fi ring.

QUANTITATIVE COMPARISON
To investigate which method produces the most reliable result, we compared the different methods with the calibrated PRC resulting from the direct method. The difference resulting from the fi ve implemented methods with the directly observed PRC can be quantifi ed by examining the Euclidean distance (by taking the mean-squared error, MSE) and correlation (with the Pearson correlation) between the obtained PRCs.
It is important to note that Galan's method and the modifi ed-Izhikevich method cannot be used with continuous fl uctuating data because they are designed to work solely with pulsed perturbation data. However, the methods intended for continuous fl uctuation data (e.g., STA, WSTA and STEP) can be used to compute the PRC from both perturbation data and fl uctuation data because the former is a simplifi ed case of the latter: instead of a continuous stream of fl uctuations only a single fl uctuation per period is injected. Hence, we can compare the PRCs resulting from the fi ve implemented method with the directly observed PRC on the perturbation data, but only the STA, WSTA and STEP on the more complex continuous fl uctuation data (i.e., noisy model data set and experimental data set).  To compare the methods with each other, we normalize the PRCs in a post-processing step. All but the STA method are defi ned (along the x-axis) over [0,2π]. We scale the time of the STA produced PRC to the same interval. Along the y-axis we normalize the PRC so that the (positive) peaks are equal to 1. Since most researchers are primarily interested in the type of the PRC (type-I vs. type-II), and because the normalization does not affect qualitative features of the PRC such as the slope and the positive and negative areas, the normalization of the results is valid.

IMPLEMENTATION DETAILS
We implemented the different methods in Python in combination with the Numpy/Scipy and Matplotlib 1 . In three methods (Galan's method, modifi ed-Izhikevich and STEP) a curve is optimized to fi t the data. Any smooth curve such as a polynomial or a Fourier series can be used for this purpose. To be consistent with the implementation in previously published studies (Galán et al., 2005;Izhikevich, 2007) we use the third expansion of the Fourier series (n = 3) in the remainder of this manuscript. Larger expansion would provide better fi ts in some cases (when there is a steep slope in the PRC) but it would also be more prone to overfi tting and hence the third expansion seems suitable. Moreover, Galan's method and (the original) Izhikevich method do not prescribe a particular optimization algorithm although Galan uses least-squares optimization. We follow his work and also use least-squares optimization in Galan's method, the modifi ed-Izhikevich method and STEP method.
For the WSTA method, the authors suggest to fi t a polynomial to the raw outcome of their algorithm because this raw output is noisy with a smaller number of spikes. For the sake of clarity we show the raw outcome to illustrate the true capabilities of this method.
All methods require confi guration of the estimated inter-spike interval ISI 2 . In our implementation of the different methods, the ISI can be given as an argument to the algorithm or automatically computed.
The automatic computation straightforwardly takes the mean and, therefore, works only for highly regular fi ring neurons. In addition we exclude ISIs that do not satisfy 0.1 ISI ISI ISI × ≤ ≤ × 2 because larger spread of ISIs generally causes the methods to fail (remember that the PRC is a characteristic of regularly fi ring neurons). Figure 3 illustrates the PRCs as computed by all implemented methods for noise-free model data with both type-I and type-II parameters. Table 2 quantifi es the difference between each PRC and the directly observed PRC by means of the MSE and Pearson correlation. For brevity, the modifi ed-Izhikevich method is labeled as 'IzhiLQ' and Galan's method is referred to as 'GalanLQ'; in both cases the LQ suffi x indicates the use of least-squares optimization. The PRCs resulting from the direct method is plotted as a dashed line and serves to calibrate the results. It can easily be verifi ed from Figure 3 that the fi ve methods compare qualitatively; it can be verifi ed from Table 2 that they are also quantitatively similar.

PERFORMANCE ON NOISE-FREE MODEL PERTURBATION DATA
The best results using this type of data are obtained by the methods designed to work with pulse-perturbation data only, namely Galan's method and the modifi ed-Izhikevich method. The results of these two methods are better in terms of the MSE and the Pearson correlation compared to the results of the methods intended for continuous fl uctuating data. Moreover, the quantitative analysis shows equal results of the modifi ed-Izhikevich and Galan's method on type-I data. This result is due to the simplicity of the curve to fi t and the low dimension of the search space (i.e., three values for a and b of the Fourier series). On type-II data, the modifi ed-Izhikevich method performs best. From methods intended for continuous fl uctuating data, the STEP method has the best performance on both the type-I and type-II model data set as the MSE and Pearson correlation indicate closest resemblance to the directly obtained PRC.

PERFORMANCE ON NOISY MODEL DATA
The noisy model data are obtained by continuous fl uctuating current injection. Here we compare the three methods designed to analyze such data. We used four different noise levels for the comparison. Each noise level has the same mean and only differs in the variance  around the mean 3 . Figure 4 illustrates the PRCs produced by the STA, WSTA and STEP method at the four levels of fl uctuations. From this fi gure it is clear that for low noise levels all methods agree qualitatively on the PRC type as all methods correctly indicate type-II PRC in the model. However, for the two higher noise levels, the WSTA method results in a (mostly) nonnegative, type-I, PRC while the two other methods correctly assess the neuron as type-II. Hence, we can say that the STA and STEP method cope better with high amplitude fl uctuations. Table 3 quantifi es the difference between the computed PRCs and the directly observed PRC. In contrast to the qualitative observation that STA and STEP perform better with fl uctuation, the WSTA method has the highest resemblance to the directly observed PRC in terms of MSE and Pearson correlation (except at the highest fl uctuation level where the STEP method obtains the best Pearson correlation). We explain this observation as follows: with higher amplitude fl uctuations, the regularity of the spikes decreases. As a result, the PRC as computed from this less regularly fi ring data is different than the PRC from the noise-free data (Tateno and Robinson, 2007;Tsubo et al., 2007). However, the PRC is a quantitative measurement of regular fi ring neurons and hence, the true PRC of a neuron is the PRC measured from spikes in the regular fi ring regime. Therefore, we compare the PRCs always to the PRC directly observed in the noise-free case although the shape of the PRC at higher fl uctuation levels might look different. The PRC produced by the WSTA method has more resemblance to the directly observed PRC although is does provide a wrong categorization of the PRC type at higher amplitudes. For low-amplitude fl uctuations we conclude that the WSTA method produces the most reliable PRC, but for higher amplitude-fl uctuations (NL = 3,4) the STEP and especially the STA method produce more reliable results (Figure 4). The reliability of the STA method has to be inferred from a visual inspection of the produced PRC and knowledge of the noise-free, directly observed PRC because the quantitative analysis always indicates low similarity between the PRC produced by the STA method and the noise-free, directly observed PRC. This phenomenon is due to the fact that the STA method is computed on the interval [0,ISI max ] and only afterwards scaled to 0,ISI ⎡ ⎣ ⎤ ⎦ . As a consequence, the STA is always relatively fl at at the beginning of the normalized interval and the PRC is shifted to later phases. Even though the interpretation of the PRC is beyond the scope of this review, the shift to higher phases does not matter for the reliability as the shift is consistent and proportional to increasing irregularity of the spike times.

PERFORMANCE ON EXPERIMENTAL DATA
The most interesting test case is the comparison of different PRC methods applied to real, experimental data. We compare the STA, WSTA and STEP method on data from layer 2/3 neurons. Figure 5 illustrates the results with two different noise levels. With 650 spikes and considerable spread of ISIs, the WSTA method produces a noisy outcome while the STA and STEP result in smooth PRCs (as they are optimized Fourier series of small expansion). The two noise levels produce PRCs that very similar; all three methods indicate a type-II PRC although at both noise levels the WSTA produced PRC is rather noisy. In the case of experimental data, there is no calibrated data to test against and hence we only look at the Pearson correlations between the different methods to assess the level of agreement between the different methods ( Table 4). For the lower noise level, the Pearson correlation between the three methods is always higher than 0.85, indicating good correspondence between the three methods. Also, for the higher noise level, the correlations stay above 0.73 which still indicates good agreement between the different methods. The STEP method has the highest Pearson correlation with the other two outcomes and can therefore be seen as a sort of 'average' of the other two methods. Additionally, the high resemblance between the three methods corroborates that these PRCs are reliable.

RELIABILITY AND PARAMETER SENSITIVITY
Here we address the reliability of the various methods for estimating the PRC. Moreover, we investigate how the reliability is affected by parameter sensitivity in the algorithm.
The reliability of the estimated PRC depends strongly on the number of spikes available for analysis. The data used to obtain the PRCs in Figures 3, 4 and 5 contain a reasonable number of spikes 4 , and the results indicate that all three methods are -to a certain extent -capable of producing reliable PRCs even under the presence of higher-amplitude fl uctuations and more diverse ISIs. However, the number of spikes available for analysis can be limited in experimental data because the experimental protocol is often not exclusively used to gather data to determine a PRC but for other scientifi c goals. Therefore, we examined the performance of the STA, WSTA and STEP method with less spikes, namely 50, 100 and 500 spikes. Figure 6 illustrates these results. The columns illustrate the PRC with (from left to right) 50, 100 and 500 spikes (the spikes are the fi rst 50, 100 and 500 of the available spikes in the data set). The PRCs from the top row use spikes from the continuous fl uctuating data set at the lowest fl uctuation level. The PRCs from the bottom row use experimental data at the second fl uctuation level.
For the model data, the WSTA and STEP methods give a fair result when using as little as 50 spikes while the STA method produces no usable PRC 5 . With 100 spikes and more, all three methods produce a reliable PRC on this set of model data with little variance in the ISIs. For the experimental data, a different view emerges. For both 50 and 100 spikes, the WSTA output is useless because of the high noise. Moreover, the STEP method wrongly classifi es the PRC as type-I with only 50 spikes; a negative part in the STEP-produced PRC using 100 or 500 spikes correctly indicates type-II PRC. With 500 spikes, all three method provide a reliable PRC. In contrast to the PRCs from the modeled data, the STA method produces fair PRC with as little as 50 spikes. This result demonstrates that the produced PRCs are infl uenced by regularity of the data and the number of spikes, rather than claiming superiority of any method of the other. With a different combination of spike times (i.e., not the initial 50, 100, or 500) the results look slightly different (not shown).
The results presented in this paper demonstrate that all methods are capable of producing reliable PRC on pulsed-perturbation data. Moreover, provided a suffi cient number of spikes are available, the STA, WSTA and STEP methods also work well with continuous fl uctuating data. This result might be, however, misleading since the reliability of the tested methods is sensitive to algorithm parameters such as the a priori estimated ISI ISI ( ) and the estimated injected Shown are the Pearson correlation between the PRCs generated by three different methods. All show strong correlation indicating similar trends in the three curves.
FIGURE 6 | The infl uence of the number of spikes on the reliability of the produced PRCs. The top and bottom row represent noisy modeled data and experimental data, respectively. We tested 50, 100 and 500 spikes and observe that the PRC type is correctly assessed by the STEP and STA method for data sets containing more than 100 spikes. The reliability increases for higher number of spikes but with 500 spikes the STA and STEP methods seem to converge on both modeled data and experimental data while the WSTA method is still very noisy on the experimental data set with 500 spikes.
step current I s (on top of which the fl uctuations are modulated). These parameters can be set manually or computed by the algorithm itself. Below we describe the effects of these two parameters on the different methods to estimate a PRC. We observed that the STA method is highly sensitive to a correct estimation of the current step, i.e., I s the injected current to make the neuron fi re regularly. Figure 7 (top left) illustrates this effect. With a very small deviation of the estimated mean from the real injected DC the outcome becomes unstable. In most cases, this effect will be evident to the researcher: when the amplitude of the STA (in the STA method) is close to 0, minor offsets in the estimated DC may shift the STA (to either all positive or all negative), and, in turn, shift the PRC resulting in a wrong indication of the PRC type. The top left panel in Figure 6 also clearly demonstrates this effect: the PRC drops steeply (which is accentuated by the normalizing of the positive peak to 1). However, when the STA is further away from 0, or when the STA crosses 0, a faulty PRC is hard to observe because the resulting PRC will resemble a PRC but will not be representative of the data. Hence, a few different settings for the estimated mean (for instance, the calculated mean from the injected signal in simulations, or, straightforwardly the injected current from the experimental setup) should be used and the resulting PRCs should be compared to PRCs produced by the other methods. In addition, we observed that the PRC produced by the STA is diffi cult to interpret for two reasons. First the PRC is always shifted on the x-axis towards later phases because it is computed on the interval [0,ISI max ] and later normalized to 0,ISI ⎡ ⎣ ⎤ ⎦ . Second, the y-axis provides the integral of the STA; it is not obvious how this relates to the exact delay or advance in spike times. Unfortunately, this second diffi culty also arises for the WSTA method where the y-axis is defi ned by the non-symmetric (see below), weighted sum of the input fl uctuations. The STEP method has a straightforward interpretation of the y-axis as it stands for the phase shift.
We observed that the WSTA method is highly sensitive to the estimated ISI (see Figure 7, top right). In effect, the ISI shifts the resulting PRC along the y-axis. This observation can be explained as follows: spikes with relatively short ISIs compared to the ISI (i) are stretched to the ISI and become straight lines with little variation, and (ii) receive a high weight in the weighted sum (Figure 7, bottom left). These two effects combined lead to shorter ISIs pulling the PRC upwards. On the other hand, relatively long ISIs compared to the ISI make the PRC resulting from the WSTA method noisy because they are compressed to fi t on the normalized interval; during this compression smooth fl uctuations become steep and FIGURE 7 | Infl uence of confi guration settings in different PRC methods. All PRCs are generated from noisy modeled data with NL = 2. Top left: the STA method is highly sensitive to errors in the estimated current step (I s ). Top right: the WSTA method is sensitive to estimates of the ISI . Bottom left: the WSTA weighing function is not symmetric and non-symmetric distributions of ISIs will lead to upward and downward drifts of the PRC. Bottom right: different outcomes of the WSTA method depend on a combination of the estimated ISI and the step current which can be either manually specifi ed as the predetermined mean (man) or computed from the input fl uctuations (set).
are added up to the PRC. For highly regularly fi ring neurons, the estimated ISI ISI ( ) is simply the mean of the ISIs ISI ( ) . However, in the more realistic cases with more variation and a skewed ISI histogram, it becomes less clear what the 'best' a priori estimated ISI should be: mean, median, or mode of the ISIs? In type-II cases this estimation might cause unreliable results because the crossing point (i.e., the point where the curve changes sign) can be shifted along the x-axis and the ratio of positive to negative surface will clearly be modifi ed, even up to the point that the negative part may became insignifi cant. In addition, Figure 7 (bottom right) illustrates the accumulated effect of different settings for the ISI and the step current I s . One problem is that all of the PRCs produced by the WSTA method shown in that panel resemble what a researcher might anticipate: some PRC estimates are of type-II while others are of type-I. The difference is caused by changing two confi guration options in the algorithm.
In general, the fact that one can configure two settings ISI and I s ( ) opens a possibility for a bias in the resulting PRC: as we just illustrated a PRC can be easily 'tuned' under the presence of considerable fluctuations to a particular PRC type by changing ISI and I s ( ) . Therefore, one might tune the settings in a way as to prove a particular PRC type. Moreover, even without predispositions about the PRC type, all of the PRC methods (and especially the methods that optimize a smooth curve) can produce PRCs that appear realistic without being representative of the data.

DISCUSSION
We reviewed fi ve different methods to determine the PRC from experimental or modeled data. Two methods are variations of the direct method and require the perturbation-stimulation protocol to gather data. The three other methods use continuous fl uctuation data, which requires the continuous injection of (a step current and) a fl uctuation into a regularly fi ring neuron. We found that on noise-free modeled perturbation data, all methods worked well and showed little difference. Moreover, the two methods requiring pulse perturbation data produced the best results, but this comes at the cost of a specialized stimulation protocol and the requirement of a large number of spikes in order to cover all the phases. In contrast, the methods that can use the continuous fl uctuation data use all available data effi ciently as random fl uctuations are by defi nition delivered at every phase and thus require less spikes. For instance, one study (Galán et al., 2005) uses 7000 (highly regular) spikes while we show that the continuous fl uctuation methods provide reliable results after a few hundreds spikes (e.g., 500). Hence, for experimental situations where little time is available it is best to use the continuous fl uctuation protocol. When an experiment is done solely for the purpose of determining a PRC, the perturbation protocol should be used.
However, we also demonstrated that the estimated PRC not only depends on the method employed, but also on the settings of the method I s and ISI

(
) and the regularity of the spikes (as altered by the amplitude of the fl uctuations). Moreover, we demonstrated that the different techniques might generate PRC curves that appear plausible but are not representative for the data. A 'panel of experts' strategy can be applied to enhance the reliability: one can run all the different methods (that are appropriate for the given data set) with slightly different confi gurations. A stable PRC for the widest range of settings can be considered the most correct. And, pitfalls such as upward or downward shifts in the PRC can be detected by trying several settings. W also suggest researchers dealing with PRC to inspect carefully the spiking data and obtain good estimates for the ISI and the DC step current before running the analysis and using any of the PRC estimation methods. Interpretation of the PRCs is beyond the scope of this review. Briefl y, however, different criteria are used to classify PRC curves into type-I curves and type-II curves. For instance, the ratio between the negative amplitude and the positive amplitude (Tateno and Robinson, 2007) or the ratio between the positive and negative surface (Tsubo et al., 2007) have been proposed as PRC categorization criteria. Moreover, some reports suggest that the exact shape of the PRC, skewness, zero-crossings and other features contain information about the underlying system (e.g., Gutkin et al., 2005;Tateno and Robinson, 2007). These features can only be reliably interpreted after obtaining a PRC in a reliable manner following guidelines and avoiding potential pitfalls as outlined above.
We hope that this review of methods for determining the PRC will motivate researchers to determine this measure of spiking behavior for the neural cell types they investigate. The PRC is a measure with considerable predictive power for the behavior of a single neuron in a network (Ermentrout, 1996). Given knowledge of a neuron's PRC and its synaptic connections (excitatory/inhibitory, waveform duration, and delay), it is possible to analytically determine its participation in population activities such as synchronization, asynchrony and beating. Specifi cally, pairs, chains or networks of neurons with a type-I PRC will synchronize when coupled by fast inhibitory synapses. For type-II neurons, synchrony ensues with excitatory synapses. Strictly speaking, these predictions only hold for homogeneous networks of regularly spiking neurons perturbed by small synaptic potentials. Nevertheless, they allow for insights about network activity emerging from single neuron properties in many activity regimes. Thus, reliable knowledge of PRCs from more neural cell types will aid the understanding of how single neurons contribute to brain function.