## Abstract

The propensity of a neuron to synchronize is captured by its infinitesimal phase response curve (iPRC). Determining whether an iPRC is biphasic, meaning that small depolarizing perturbations can actually delay the next spike, if delivered at appropriate phases, is a daunting experimental task because negative lobes in the iPRC (unlike positive ones) tend to be small and may be occluded by the normal discharge variability of a neuron. To circumvent this problem, iPRCs are commonly derived from numerical models of neurons. Here, we propose a novel and natural method to estimate the iPRC by direct estimation of its spectral modes. First, we show analytically that the spectral modes of the iPRC of an arbitrary oscillator are readily measured by applying weak harmonic perturbations. Next, applying this methodology to biophysical neuronal models, we show that a low-dimensional spectral reconstruction is sufficient to capture the structure of the iPRC. This structure was preserved reasonably well even with added physiological scale jitter in the neuronal models. To validate the methodology empirically, we applied it first to a low-noise electronic oscillator with a known design and then to cortical pyramidal neurons, recorded in whole cell configuration, that are known to possess a monophasic iPRC. Finally, using the methodology in conjunction with perforated-patch recordings from pallidal neurons, we show, in contrast to recent modeling studies, that these neurons have biphasic somatic iPRCs. Biphasic iPRCs would cause lateral somatically targeted pallidal inhibition to desynchronize pallidal neurons, providing a plausible explanation for their lack of synchrony in vivo.

- phase coupling
- biophysical modeling
- 555 Timer
- perforated patch
- basal ganglia
- Fourier transform
- integrate-and-fire neuron

the infinitesimal phase response curve (iPRC) is a fundamental tool in the study of neuronal synchrony (Schultheiss et al. 2012b). It describes how the phase of an oscillator is affected by brief and vanishingly weak perturbations as a function of the phase at which they were delivered. With knowledge of the synaptic coupling, the emergent network dynamics are determined by the iPRC of a neuron (Ermentrout 1996; Ermentrout and Terman 2010; Hansel et al. 1995; Kuramoto 1984; Van Vreeswijk et al. 1994). There are two qualitatively different iPRCs. Neurons for which phase can only be advanced by depolarizing inputs exhibit a type I iPRC, which is a strictly positive function of phase, whereas neurons for which phase can also be retarded exhibit a biphasic (type II) iPRC (Ermentrout 1996).

In practice, measuring the iPRC is a daunting task due to the intrinsic discharge variability of neurons (Netoff et al. 2012; Ota et al. 2011; Schultheiss et al. 2010), which introduces a jitter in the phase of the oscillator that is commensurate with the phase change caused by the perturbations. Increasing the amplitude of the perturbation runs the risk of distorting the estimate of the iPRC (Acker et al. 2003; Farries and Wilson 2012a,b; Netoff et al. 2012; Ota et al. 2011; Schultheiss et al. 2010). The problem is exacerbated with biphasic iPRCs because the negative lobe is typically considerably smaller than the positive one and therefore harder to detect reliably (Goldberg et al. 2007; Schultheiss et al. 2012b).

Two approaches have been taken to circumvent this experimental difficulty. One is to implement elaborate nonlinear algorithms to improve the estimates (Netoff et al. 2012; Ota et al. 2011; Phoka et al. 2010; Torben-Nielsen et al. 2010) or to use the effective indirect method that uses the spike-triggered average (STA) of noisy injected input to calculate the iPRC (Burton et al. 2012; Ermentrout et al. 2007; Ota et al. 2009). The other is to simulate biophysical models of neurons, where the noise level can be controlled (Farries and Wilson 2012a,b; Schultheiss et al. 2010, 2012a). Recent modeling studies have concluded that the somatic iPRC of globus pallidus (GP) neurons is monophasic (Schultheiss et al. 2010, 2012a).

Here, we propose a novel and natural method to estimate the iPRC that uses its periodicity (Galán et al. 2005; Netoff et al. 2012). We show, first analytically, then numerically with biophysical models, and finally empirically with a physical oscillator of a known design that the iPRC can be reconstructed reliably by perturbing the oscillator with extremely weak and charge-balanced sinusoidal waveforms. We apply the methodology to mouse cortical pyramidal neurons recorded in whole cell configuration and show that it correctly identifies their monophasic iPRCs. Applying it to GP neurons recorded in the perforated-patch configuration shows that the somatic iPRC of these neurons is biphasic.

## MATERIALS AND METHODS

#### Direct estimation of the spectral (Fourier) modes of the iPRC.

Consider a dynamic system (e.g., a biophysical neuronal model) represented by a set of differential equations written in vector notation as: (1)where the vector **x** is the vector of dynamic variables, and **F** is the vector notation of the right-hand side (r.h.s.) of these equations. We assume that the first component of the vector is voltage (i.e., *x*_{1} = *V*) with an equation of the form d*V*/d*t* = (...)/*C* and that there is a stable periodic (limit-cycle) solution to this system denoted **x**_{0} with period *T* [i.e., **x**_{0}(*t*) = **x**_{0}(*t* + *T*)]. It is possible to define a phase variable *ϕ* on and in the neighborhood of the limit cycle in the phase space of the system such that the phase along the limit cycle advances at a constant rate (Kuramoto 1984): (2)Expanding *Eq. 2* with partial derivatives yields: (3)The gradient of the phase that appears in *Eq. 3* is none other than the vector of the iPRCs of the various dynamic variables, denoted **Z** (· denotes a scalar product).

Now assume that a small *n*th mode sinusoidal current of amplitude Δ*I* and period *T* is added to the voltage equation in the dynamic system in *Eq. 1*: (4)where the Kronecker delta denotes that the perturbation is only to the first (i.e., voltage) equation. We can substitute *Eq. 4* into the left-hand side of *Eq. 3* and integrate it over one period of the limit cycle: (5)where *Z*_{V} is the iPRC of the voltage.

We assume that the perturbation is small enough such that the altered trajectory of the oscillator is still in the neighborhood of the unperturbed solution so that the phase is well-defined. Therefore, the first term on the r.h.s. of *Eq. 5* equals *T* by *Eq. 3*, but the second term must add a small time increment, Δ*T*_{n}^{S}, where the superscript *S* denotes that a sine wave was used, and the subscript *n* is the mode number. We therefore have: (6)where is, by definition, the *n*th sine transform of the iPRC [of course, the same analysis can be repeated for the direct current (DC) and cosine transforms].

In summary, injecting a weak current waveform shaped like a sine or cosine function with a period matching that of the unperturbed oscillator will alter the phase of the oscillation by an amount that is proportional to the corresponding sine or cosine transform of the iPRC. The iPRC can be reconstructed from these measurements as: (7)

#### Biophysical neuronal models.

We simulated two well-known biophysical models: the original Hodgkin-Huxley model (Hodgkin and Huxley 1952), shifted by −65 mV (with *I*_{app} = 11 μA/cm^{2}), and the Traub model as it appears in Ermentrout (1998; with *I*_{app} = 0.2 μA/cm^{2}). The models were solved numerically on XPPAUT (Ermentrout 2003). The iPRC was measured in three different ways: *1*) application of brief current perturbations at various phases of the cycle; *2*) using the adjoint method (Ermentrout and Terman 2010; Williams and Bowtell 1997); and *3*) spectral reconstruction. In the latter case, small amplitude current waveforms in the shape of the various cosine and sine modes were injected into the model neuron, and the effect on the timing of the ensuing spike was measured in relation to the unperturbed period and used to calculate the iPRC with *Eq. 7*. To verify the robustness of the method to neuronal jitter, we added in separate simulations a Wiener (white noise) process to the models (Ermentrout 2003) to generate a jitter in the firing of these neurons that is comparable, in the sense of having the same coefficient of variation (CV), with that observed in the GP neurons recorded in perforated-patch configuration (see below). The amplitude of the noise term needed to generate the target CV for the Hodgkin-Huxley neuron was 2.25 μA/cm^{2}, and for the Traub model it was 0.35 μA/cm^{2}. Note, however, that in the simulations, these values are multiplied by a normal random variable divided by , where d*t* is the size of the integration step. Thus when the solver advances the solution by multiplying the deterministic terms in the differential equations by d*t*, the noisy term is effectively multiplied by (Ermentrout 2003). This guarantees that regardless of step size, the total variance of the process over a given time interval is constant. These simulations were conducted in MATLAB (The MathWorks, Natick, MA). We compared the robustness to noise of the spectral estimator with that of the traditional temporal estimator by simulating the same number of trials (700) using either type of perturbation (brief pulse at random phases for the temporal estimate and extended harmonic mode for the spectral estimator). In the case of the temporal estimator, the 700 brief perturbations were distributed evenly across all phases of the mean interspike interval (ISI) of the noisy model neuron, and the scatterplot of resulting phase delays as a function of the phases at which the perturbations were delivered was calculated. The scatterplot was least-squares-fit to a model composed of a sum of 5 modes: DC, fundamental (sine and cosine), and 2nd harmonic (sine and cosine; Galán et al. 2005; Goldberg et al. 2007). In the case of the spectral estimator, the 700 trials were divided among the 1st 7 modes (DC and both cosine and sine modes of the 1st 3 harmonics) so that 100 trials were used per mode. The average phase delay generated in the noisy neuron was calculated for each mode and was used as the prefactor of that mode in the spectral reconstruction (*Eq. 7*).

#### An integrated-circuit oscillator: the astable mode of the 555 Timer.

The 555 Timer is an integrated circuit (IC) that, when configured in the so-called astable mode, can act as an oscillator for which output (pin 3) flip-flops between 0 V and the supply voltage (*V*_{CC}), which in our experiments was 3 V (http://www.ti.com/lit/ds/symlink/tlc551.pdf). The astable mode circuit was constructed with a TLC551 LinCMOS Timer (Texas Instruments, Dallas, TX) on a solderless breadboard by wiring the IC as shown below. All components of the circuits were purchased from DigiKey (Thief River Falls, MN). The voltage on the capacitor can be measured at pins 2 and 6 (which are short-circuited to each other) and nominally has the form of: (8)where *T*_{R} = ln(2)(*R*_{1} + *R*_{2})*C* is the rise time, and the period *T* equals *T*_{R} + ln(2)*R*_{2}*C*. In our experiments, the empirical period was 93.3 ms, which was close to the theoretical period of 95.0 ms calculated from the nominal values of the resistors and capacitor. The iPRC for a one-dimensional system is given by: (9)yielding a biphasic iPRC: (10)At the limit of *R*_{1} << *R*_{2} (i.e., duty cycle of 50%), which applies in our case, we have:
(11)

The voltage on the capacitor (pins 2 and 6) and the output voltage (pin 3) were recorded simultaneously, and current was delivered to pins 2 and 6, situated between the capacitor *C* and the resistor *R*_{2}, with a MultiClamp 700B amplifier (Molecular Devices, Sunnyvale, CA). Signals were digitized at 10 kHz and logged onto a personal computer (PC) with Clampex 9.2 software (Molecular Devices). A transistor-transistor logic (TTL) pulse was delivered to the reset pin (4) to control the precise timing of the onset of the oscillations. The current injected into pins 2 and 6 was in the form of either brief 5-ms, 200-nA pulses or DC and harmonic mode waveforms with the 555 Timer empirical period and with an amplitude of 20 or 50 nA. The waveforms of the harmonic modes were generated in IGOR Pro 5 (WaveMetrics, Portland, OR).

#### Whole cell and perforated-patch recordings from cortical and pallidal neurons, respectively.

Experimental procedures were approved by the Northwestern University Animal Care and Use Committee. Five-week-old C57BL/6J mice of either sex were deeply anesthetized with ketamine-xylazine and perfused transcardially with ice-cold modified artificial cerebrospinal fluid (ACSF) bubbled with 95% O_{2}-5% CO_{2} and containing (in mM): 2.5 KCl, 26 NaHCO_{3}, 1.25 Na_{2}HPO_{4}, 0.5 CaCl_{2}, 10 MgSO_{4}, 0.4 ascorbic acid, 10 glucose, and 210 sucrose. The brain was rapidly removed, blocked in the sagittal plane, and sectioned at a thickness of 275 μm in ice-cold modified ACSF. Slices containing the cortex or GP were then submerged in ACSF bubbled with 95% O_{2}-5% CO_{2} and containing, in mM, 2.5 KCl, 126 NaCl, 26 NaHCO_{3}, 1.25 Na_{2}HPO_{4}, 2 CaCl_{2}, 2 MgSO_{4}, and 10 glucose and stored at room temperature for ≥1 h before recording.

The slices were transferred to the recording chamber mounted on an Olympus BX51 upright, fixed-stage microscope and perfused with oxygenated ACSF at 32°C. A ×60, 0.9-numerical aperture water-immersion objective was used to examine the slice using standard infrared differential interference contrast video microscopy.

For the perforated-patch recordings, thin-wall borosilicate glass pipettes with no filament were front-filled with a prefiltered solution containing (in mM): 135.5 K-Me-SO_{4}, 5 KCl, 2.5 NaCl, 5 Na-phosphocreatine, 10 HEPES, 0.2 EGTA, 0.21 Na_{2}GTP, and 2 Mg_{1.5}ATP (pH 7.3 with KOH, 280–290 mosmol/kgH_{2}O). They were then back-filled with the same solution that also contained 1.5 μg/ml gramicidin B (Sigma, St. Louis, MO) with 0.1% DMSO. Pipette resistance was 2.5–3 MΩ, and >5-GΩ seals were attained on GP somata. Perforation was monitored until spikes exceeded −10 mV, after which the recordings began. For the whole cell configuration recordings, pipettes with filaments were back-filled with the same solution but lacking the gramicidin B. Recordings were obtained with a MultiClamp 700B amplifier (Molecular Devices). Signals were digitized at 10 kHz and logged onto a PC with Clampex 9.2 software (Molecular Devices).

#### Estimation of iPRCs from electrophysiological data.

The iPRCs were estimated for cortical or GP neurons in the standard method by injecting 5-ms long, 40- or 10-pA pulses, respectively, at random phases of the discharge of the neuron. The scatterplots of the phase response were fit to the sum of the first five harmonic modes (Galán et al. 2005; Goldberg et al. 2007). To estimate the spectral modes directly, a file containing the DC or one of the first two harmonic modes (sine and cosine of each) was generated using IGOR Pro 5 (WaveMetrics) and used as a current command with amplitude of 10 pA for cortical pyramidal neurons and of 1–2 pA for GP neurons. A cycle of the waveform was injected, flanked by a 0-current injection. This was repeated ∼50 times per waveform. From these, we chose for the analysis the (2–8) trials where the cell spiked within a window of 5 ms before the onset of the current waveform so that the waveform injection could be considered to overlap the ISI and averaged the resultant phase delays. The perturbed ISI resulting from this current waveform injection was normalized by the average ISI, denoted *T*, and estimated visually using the oscilloscope mode of the Clampex 9.2 software before application of the perturbations. Importantly, using visual estimation is justified because it led to an error of up to 15% compared with the average ISI measured a posteriori offline. Such an error in the estimation of *T* is inconsequential to the resulting spectral estimate of the iPRC because it only introduces a second-order correction to *Eq. 7*. This is because the contribution of this small 15% error to the spectral estimate is always multiplied by the size of the perturbation, which is also small. The following variant of *Eq. 7* was used to calculate the estimate: (7A)For cortical cell recordings, a constant 100-pA current was applied to cause them to discharge rhythmically. For some GP cells, a constant 5- to 10-pA current was applied to regularize their autonomous discharge.

#### Statistical validation of the biphasic nature of empirical iPRCs.

We have previously used the ratio of the amplitude of the fundamental mode to the amplitude of DC mode, denoted *R*_{1}, to determine whether an iPRC is biphasic (Goldberg et al. 2007). In the current notation, *R*_{1} is defined as *R*_{1} ≡ 2/Δ*T*_{0}^{C}. When *R*_{1} > 1 (i.e., the amplitude of the fundamental mode is larger than that of the DC mode), the iPRC will typically attain negative values (assuming the other spectral modes fall off faster), thus exhibiting a type II biphasic structure. However, it is possible that random fluctuations in the estimates of the amplitudes of the these modes of the iPRC (i.e., Δ*T*_{0}^{C}, Δ*T*_{1}^{C}, and Δ*T*_{1}^{S}) will cause *R*_{1} to be larger than unity spuriously. To estimate the distribution of *R*_{1} under this null hypothesis, we measured the amplitudes of the modes resulting from randomly chosen trials in which the spikes were not properly aligned with the onset of the DC or sinusoidal perturbation. We then used a two-tailed Wilcoxon rank-sum test to determine whether *R*_{1} calculated from the correctly aligned trials differed from *R*_{1} calculated according to the null hypothesis.

All offline analysis and rendering of results was done in MATLAB.

## RESULTS

#### Direct estimation of spectral modes of the iPRC.

Because the iPRC is a periodic function, it can be reconstructed as the sum of its Fourier modes. How, then, can the Fourier modes be measured directly? In materials and methods, we show that for an arbitrary oscillator, application of a weak but extended sinusoidal perturbation that spans the period of the oscillator will result in a phase change in the oscillator that is proportional to the Fourier transform of the iPRC associated with that same waveform. Let us exemplify this result with an integrate-and-fire (I&F) neuron (Fig. 1). The unperturbed voltage trajectory, *V*(*t*) = 1 − *e*^{−t/τ}, follows an exponential trajectory, with time constant *τ*, and reaches threshold *θ* at time *T*, which is the period of the I&F neuron. This threshold crossing is represented by the following equation: (12)

Injection of a current waveform in the shape of the fundamental sine mode with amplitude Δ*V* causes the perturbed voltage to reach the threshold at a later time, denoted by iPRC convention as *T*−Δ*T*. If we write the equation of the perturbed voltage and integrate it up to time *T*−Δ*T*, we can generate another equation: (13)Solving *Eqs. 12* and *13* up to first order in Δ*V* or Δ*T* yields: (14)The r.h.s. of *Eq. 14* equals half of the fundamental sine transform of the function *e*^{t/τ}. Applying *Eq. 9* to the voltage trajectory of the I&F neuron demonstrates the well-known result (Hansel et al. 1995) that *e*^{t/τ} is indeed the iPRC of the I&F neuron (Fig. 1). This iPRC is strictly positive and therefore type I, which means that brief depolarizing pulses can only advance the next spike in the I&F neuron. When applying an extended perturbation such as a sine wave, the early depolarizing portion of the sine wave will tend to speed up the approach to threshold, whereas the late hyperpolarizing portion will tend to delay it. The net effect is a delay of the spike, because the exponentially shaped iPRC is larger during the hyperpolarizing portion of the fundamental sine wave, and that portion is weighted by the larger portion of the monotonic iPRC.

Thus the I&F neuron example demonstrates that applying a perturbation shaped like the fundamental sine mode results in a phase delay that is proportional to the fundamental sine transform of the iPRC. In addition, it suggests that even though the iPRC is defined for vanishingly small and brief perturbations, on a practical (and mathematically sound) level the iPRC can be reconstructed from the response of the oscillator to a small and discrete set of extended stimuli (DC and sinusoidal waveforms). How many modes are necessary for a good reconstruction? Of course, this depends on the type of neuron. For neurons like the I&F neuron that have a discontinuous iPRC, the power in each mode should generally decay slowly (algebraically) with mode number so that a good reconstruction will require many modes. However, if the iPRC is a relatively smooth function, very few modes should capture the structure of the iPRC. To see this, we turned to biophysical neuronal models.

#### iPRCs of biophysical models can be reconstructed from the first few harmonic modes.

We used two well-known models to demonstrate the low dimensionality of the iPRCs. The first is the original Hodgkin-Huxley model (1952) shifted by −65 mV. The second is the Traub model (Ermentrout 1998). The two models differ in the nature of their iPRCs. The Traub model (Fig. 2*A*) possesses a strictly positive (type I) iPRC, which is simple to understand: depolarizing pulses bring the membrane voltage closer to threshold and therefore can only advance the next spike. The Hodgkin-Huxley model (Fig. 2*B*) gives rise to a biphasic (type II) iPRC, meaning that brief perturbations introduced early in the cycle delay the next spike. This delay occurs because the so-called h gate (which represents sodium-channel availability) has a large dynamic range in the voltage range of the oscillations. The depolarization slightly reduces the availability of sodium channels (or equivalently, inactivates them; Goldberg and Reynolds 2011), resulting in an effectively higher spike threshold, which takes longer to reach. In contrast, in the Traub model, the h gate is saturated throughout the entire subthreshold range (data not shown). Therefore, depolarizing perturbations are unable to reduce significantly the availability of sodium channels to reverse the direct effect of depolarizing, which is to approach threshold sooner.

We calculated the iPRC of these models in three different ways: *1*) with brief perturbations applied at the various phases of the cycle; *2*) numerically using the adjoint method (Ermentrout and Terman 2010; Williams and Bowtell 1997); and *3*) using the spectral reconstruction (*Eq. 7*). In both models, the DC mode plus four cosine and four sine modes were sufficient to reconstruct closely the true iPRC.

#### Empirical spectral reconstruction of an iPRC of an IC oscillator.

As its name suggests, the iPRC defines the response of a dynamic system to vanishingly small and brief perturbations. This calls its experimental validity into question for the following reason. Although numerical simulations can be studied for extremely small perturbations, “real-life” physical oscillators are subject to noise and/or have thresholds below which they are insensitive to perturbations. Hence, the predicted response to vanishingly small perturbation may have little validity in actual physical systems. Conversely, a phase response curve estimated with realistic perturbations may diverge significantly from the iPRC (Acker et al. 2003; Ota et al. 2011; Phoka et al. 2010; Schultheiss et al. 2010). It is therefore necessary to test whether the spectral estimation of the iPRC is valid in a real physical system (and not only in analytical or numerical models). Furthermore, the traditional temporal method to estimate the iPRC employs brief perturbations. In contrast, the spectral estimation requires the use of perturbations that span the entire period of the oscillator, which, one might anticipate, could give rise to an estimate the deviates from the true iPRC.

Therefore, empirically validating the spectral reconstruction should employ a physical system with a known iPRC. We therefore applied the methodology to a physical oscillator with a fully determined design.

We chose to use the 555 Timer IC (see materials and methods). When built into a circuit with two external resistors and a capacitor (*R*_{1}, *R*_{2}, and *C* in Fig. 3*A*), the output of the 555 Timer (pin 3) can be made to flip-flop, in the so-called astable mode, between 0 V and the supply voltage (3 V) with a period that depends on the resistance and capacitance of the external components (Fig. 3*B*). The voltage on the external capacitor (which alternately decays exponentially to the 2 output levels) can be measured and perturbed with external current injections at pins 2 and 6, which are short-circuited to each other (Fig. 3*B*). Because the voltage trajectory on the capacitor is known by design (*Eq. 8*), the predicted iPRC of the astable mode of the 555 Timer can be calculated analytically (*Eqs. 10* and *11*). This iPRC is biphasic and therefore type II.

The traditional method to measure the iPRC with brief current pulses revealed the type II nature of the iPRC: if the depolarizing pulse is delivered while the voltage on the capacitor is depolarizing, the upper threshold is reached sooner; conversely, if it is delivered while this voltage is hyperpolarizing, the lower threshold is reached later (Fig. 3*B*). Additionally, the iPRC calculated this way agreed well with the theoretical curve (Fig. 3*C*). Note, however, that there was a bias estimating the negative lobe of the iPRC. This resulted from the fact that the perturbations were depolarizing pulses. As expected from symmetry, when hyperpolarizing pulses were used instead, the bias shifted to the positive lobe (data not shown). This bias thus represents a deviation from the iPRC resulting from the finite size of the perturbation.

Because the iPRC of this system is discontinuous, the spectral reconstruction using the first few harmonics cannot be expected to trace the theoretical one perfectly. Nevertheless, using only the first five harmonics captures the essential biphasic nature of the theoretical curve and closely matches the sum of the first five harmonics of the theoretical curve (Fig. 3, *D–F*). Moreover, note that there is no visible bias in the spectrally reconstructed iPRC (Fig. 3*E*) compared with the traditional estimate (Fig. 3*C*). This is undoubtedly because the extended harmonic waveforms are charge-balanced (except for the DC mode), whereas the traditional pulses are not. Another advantage of the spectral reconstruction is that order-of-magnitude weaker perturbations suffice to generate it. The 555 Timer data indicate that the spectral reconstruction is a viable experimental method to estimate the essential structure of an empirical iPRC, and may even surpass the traditional temporal method, because it uses weaker perturbations and thus should distort the estimate of the iPRC less.

#### The spectral estimator is robust to physiological scale jitter in model neurons.

To determine whether the spectral estimator is valid under noisy conditions, we added to the Traub and Hodgkin-Huxley models a white-noise (Wiener) process that was tuned to give rise to a CV of 0.15 (Fig. 4, *A* and *D*, respectively). This value was chosen because it is the median CV measured from the GP neurons recorded in perforated-patch recordings (see below). To calculate the temporal estimator, perturbations were distributed uniformly over the mean ISI of the model neurons (which was approximately equal to the ISI of their corresponding noiseless neurons), and their responses were recorded as a large cloud of points. The temporal estimate of the iPRC was the periodic function that best fit this scatterplot (Galán et al. 2005). The temporal iPRCs often differed significantly from the iPRCs of the corresponding noiseless models (compare Fig. 2 with Fig. 4, *B* and *E*).

To calculate the spectral estimator, we divided the same number of trials evenly among the 1st 4 harmonic modes (7 modes when sines and cosines are counted separately) and estimated them by averaging the response to each harmonic perturbation and using that average response as the prefactor of the associated mode in the spectral reconstruction (*Eq. 7*). The 1st 3 modes are shown in Fig. 4, *C* and *G*. The most obvious difference between the 2 models is the response to the DC (constant) perturbation. For the type I (Traub model) neuron, this generated the largest response guaranteeing that the overall phase response was positive (Fig. 4*C*), whereas the response to the DC perturbation is much smaller in the type II (Hodgkin-Huxley model) neuron, giving rise to negative lobes due to the contribution of the higher modes (Fig. 4*G*).

Repeating these simulations several times demonstrated that the nonnegative nature of the type I response was quite robust and that a biphasic region was always present in the latter part of type II response (data not shown). Both of these features are shared with the corresponding noiseless model neurons. To quantify how well the two methods performed in the noisy models, we repeated the simulations several times for both the spectral and temporal estimators and calculated the correlation coefficient between the waveforms of the noisy and noiseless models from the outcome of each simulation. For both the Hodgkin-Huxley and Traub models, the median of the distribution of the correlation coefficients for the spectral method was larger than that of the temporal estimator, implying that the spectral estimator performed slightly better at capturing the underlying iPRC (Fig. 4, *D* and *H*).

#### Spectral reconstruction reveals the monophasic nature of cortical iPRCs and the biphasic nature of pallidal iPRCs.

Armed with the theoretical derivation of the spectral reconstruction and its successful application not only to noiseless and noisy numerical models, but also to a real-life physical oscillator buttresses the validity of the approach and created confidence in applying it to a physiological system such as a real neuron. To validate the method in a physiological setting, we used it to estimate the somatic iPRC of cortical pyramidal neurons recorded in whole cell configuration. Previous studies have shown that their somatic iPRC is essentially type I (Goldberg et al. 2007; Reyes and Fetz 1993). Brief 40-pA perturbations typically advanced the ensuing spike (Fig. 5*A*) and were used to construct the traditional iPRC scatterplot. The harmonic function fit to the scatterplot revealed the essentially type I response (Fig. 5*B*).

DC and sinusoidal current waveforms, for which amplitude of 10 pA was weaker than that of the brief pulses, and for which period approximated the observed average ISI, were then injected into the cell, and the ensuing phase change was recorded. Use of an extended sinusoidal current delayed the next spike, presumably because the hyperpolarizing portion had a more prominent effect on the phase than the depolarizing portion due to the type I structure of the iPRC (Fig. 5*C*). The average response to the lowest modes reveals that injection of the DC and fundamental cosine mode currents on average advanced the next spike (positive prefactors), whereas the sine-wave current delayed it (hence its negative prefactor). The similarity between the type I iPRC resulting from the spectral estimate and the iPRCs observed previously (Goldberg et al. 2007; Reyes and Fetz 1993) was striking (Fig. 5*D*).

Finally, GP neurons were recorded in the perforated-patch configuration, which best maintains the stability and reliability of the intracellular recording. Brief 10-pA pulses could sometimes delay the ensuing spike (Fig. 6*A*, *top*). Extended sinusoidal waveforms typically delayed the next spike (Fig. 6*A*, *bottom*), as is evidenced by the negative prefactor of the fundamental sine mode. The latter part of the iPRC contributes a net delay because the sine-wave current is negative (hyperpolarizing) whereas the iPRC is positive in that region as in the case of a type I response. In the GP neurons, the initial positive (depolarizing) part of the sine wave adds a net delay because it is multiplied by a negative region of the iPRC. In some cells, the traditional estimator of the iPRC resembled the spectral one (Fig. 6*B*, *top*), but in others it differed (Fig. 6*C*, *top*). However, the spectral reconstruction (*Eq. 7A*) consistently (*n* = 5/5 cells) revealed a negative lobe (following the spike) in the iPRC of GP neurons (Fig. 6, *B* and *C*, *bottom*). Using a statistic that quantifies the significance of the negative lobe (denoted *R*_{1}) demonstrated that the iPRCs of the GP neurons were truly biphasic compared with iPRC resulting from randomly shuffled unused trials (Fig. 6*C*, *inset*). As with the IC oscillator and the cortical pyramidal cell, order-of-magnitude weaker perturbations (in comparison with those used to generate the temporal estimator) sufficed to generate the spectral estimator.

## DISCUSSION

#### Comparison of the spectral method to other methods of estimating the iPRC.

We presented a method to calculate the iPRC of an oscillator by estimating its first few harmonic modes. Each harmonic mode is estimated by injecting a weak current perturbation shaped like the sine or cosine waveform (with the same period as the oscillator) that is associated with that mode and measuring the resultant phase shift. The set of phase shifts are then used as the prefactors of the Fourier series of sines and cosines to generate the spectral estimate of the iPRC.

There are several advantages to the spectral method over the traditional temporal method of applying brief perturbations at various phases of the oscillation. First, the amplitude required to perturb the oscillator is typically smaller when using the spectral estimate. We saw this to be the case in the numerical simulations of the biophysical (Hodgkin-Huxley and Traub) models (data not shown) as well as for the physical (the 555 Timer IC) and physiological (cortical and GP neurons) oscillators. The reason order-of-magnitude weaker perturbations suffice is that they are extended and span the whole period of the oscillation. This advantage is paramount because it improves the validity of the estimator that is supposed to apply for vanishingly small perturbations (which are in reality only theoretical constructs). In other words, the weakness of the perturbation safeguards against deviations of the measured iPRC from the ideal one. Moreover, these perturbations integrate the response across all phases of the oscillations, thereby introducing some averaging of the response weighted by the structure of the harmonic mode, which is expected to reduce the variability of the response.

Second, except for the DC mode, the cosine and sine waveforms are charge-balanced, so they reduce any bias in the estimator arising from using the brief unbalanced perturbations. (Incidentally, the bias can be reduced by reversing the sign of the unbalanced perturbations and splitting the difference between the response to the positive and negative perturbations.) The STA method benefits from this lack of bias as well (Burton et al. 2012; Ermentrout et al. 2007; Ota et al. 2009).

Third, for biological oscillators, obtaining a good representation of the iPRC using the temporal estimator requires collecting many (100–200) trials to sample randomly the various phases of the oscillator. In contrast, the iPRC of a continuous dynamic system (unlike the I&F neurons and the 555 Timer oscillator) is continuous. Fourier theory states that the transforms of a smooth function generally drop off more rapidly than transforms of functions that are discontinuous. Therefore, a small number of harmonics should capture the essential structure of a smooth iPRC (Galán et al. 2005; Goldberg et al. 2007). Thus the same number of trials could in principle be used more wisely by dividing them among a small set of predesigned perturbations (we used 5). In our implementation of the spectral method, we did not take advantage of this potential because we had to collect many trials and then discard the majority in which the harmonic mode waveform did not properly overlap the ISI of the neuron. However, the scheme could be improved by using real-time feedback to trigger the harmonic current injection precisely on the occurrence of the spike. Then, the efficiency of the method would be maximized because each trial could be made to count. We saw this to be the case in the simulated noisy neurons because the estimates arising from the spectral method were more similar to the iPRCs of the noiseless neurons than the estimates arising from the temporal method.

Finally, the spectral method is also appealing simply because it capitalizes on the iPRC periodicity (Galán et al. 2005; Netoff et al. 2012). The set of cosine and sine waveforms form an uncorrelated (orthogonal) set of perturbations that naturally characterize the responsiveness of the oscillator. For example, a strictly positive (type I) iPRC, such as observed in the Traub model, “resonates” strongly with the DC perturbation. Conversely, a biphasic (type II) iPRC, such as observed in the Hodgkin-Huxley model, responds strongly to the fundamental harmonic perturbation.

Recent years have seen an upsurge in published methods to estimate the iPRC (Ermentrout et al. 2007; Farries and Wilson 2012a,b; Izhikevich 2007; Netoff et al. 2005, 2012; Ota et al. 2009, 2011; Phoka et al. 2010; Torben-Nielsen et al. 2010). Many of these employ elaborate nonlinear algorithms, whereas some suggest methods to reduce the intrinsic jitter using real-time feedback such as dynamic clamp (Netoff et al. 2012). Additionally, most rely on the standard method to estimate the iPRC using brief perturbations. One interesting and robust exception is the method that derives the iPRC from the STA of continuously injected noise (Ermentrout et al. 2007; Ota et al. 2009). The STA method, which is also linear, has been shown recently to work very well (Burton et al. 2012). Another method uses white noise to generate and estimate of the iPRC directly (Izhikevich 2007; Netoff et al. 2012). Despite its intrinsic advantages, to the best of our knowledge, the spectral estimator has not been suggested before, even outside of the neuroscience literature. The method it most closely resembles is the one that employs nominally white noise to generate an estimate of the iPRC directly (Izhikevich 2007; Netoff et al. 2012). As in system identification methods, white-noise stimuli engage all the spectral modes equally and simultaneously, whereas our method focuses on the first few modes one at a time.

#### The empirically validated spectral estimate shows that GP neurons have biphasic iPRCs.

As mentioned in the Introduction, the iPRC has been criticized as a marginally valid construct because many attempts to measure it empirically lead to estimates that are distorted by the finite size of the perturbation and that do not comply with the linearity requirement (Acker et al. 2003; Ota et al. 2011; Phoka et al. 2010; Schultheiss et al. 2010). In this work, not only have we shown that the weak, charge-balanced, and averaging nature of the perturbations make the spectral estimate less susceptible to these issues, but also we have validated it empirically on a physical oscillator with a known iPRC (e.g., the 555 Timer IC) in addition to numerical and analytical models. Additionally, we validated it on somata of cortical pyramidal cells, which previous studies have shown to have a type I iPRC (Goldberg et al. 2007; Reyes and Fetz 1993). We are therefore confident that it can be implemented meaningfully in neurons.

Our analysis has shown that the somatic iPRCs of GP neurons are biphasic with a negative lobe positioned shortly after the spike. A recent modeling study concluded that somata of GP neurons should exhibit a monophasic iPRC (Schultheiss et al. 2010). Two membrane currents that commonly give rise to a type II response are the hyperpolarization-activated, cyclic nucleotide-gated (HCN) cation current and the voltage-gated sodium current (Goldberg et al. 2007; Hutcheon and Yarom 2000). Both of these currents are present in GP neurons (Chan et al. 2004; Hanson et al. 2004; Mercer et al. 2007) and were included in that study (Schultheiss et al. 2010). Further experimental and modeling work is needed to determine the potential involvement of these currents in generating a type II somatic iPRC in GP neurons. A network of neurons possessing such an iPRC would tend to be desynchronized by fast inhibition (Van Vreeswijk et al. 1994). As lateral pallidopallidal inhibitory synapses are located preferentially on perikarya (Kita and Kitai 1994; Mallet et al. 2012; Millhouse 1986; Sadek et al. 2007; Sato et al. 2000; Sims et al. 2008), having a type II response could explain why the discharge of GP neurons is normally uncorrelated in vivo (Nini et al. 1995).

## GRANTS

This work was supported by National Institute of Neurological Disorders and Stroke Grant P50-NS-047085 to D. J. Surmeier. J. A. Goldberg was supported by the IDP Foundation and CHDI.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

J.A.G. and J.F.A. conception and design of research; J.A.G. performed experiments; J.A.G. analyzed data; J.A.G. interpreted results of experiments; J.A.G. prepared figures; J.A.G. drafted manuscript; J.A.G. and D.J.S. edited and revised manuscript; J.A.G., J.F.A., and D.J.S. approved final version of manuscript.

## ACKNOWLEDGMENTS

We thank Maurizio Di Veroli, Nir Schain, and Dr. Luis Carrillo-Reid for useful discussions.

- Copyright © 2013 the American Physiological Society