## Abstract

Experimental evidence indicates that the response of subthalamic neurons to excitatory postsynaptic potentials (EPSPs) is well described by their infinitesimal phase response curves (iPRC). However, the factors controlling the shape of that iPRC, and hence controlling the way subthalamic neurons respond to synaptic input, are unclear. We developed a biophysical model of subthalamic neurons to aid in the understanding of their iPRCs; this model exhibited an iPRC type common to many subthalamic cells. We devised a method for deriving its iPRC from its biophysical properties that clarifies how these different properties interact to shape the iPRC. This method revealed why the response of subthalamic neurons is well approximated by their iPRCs and how that approximation becomes less accurate under strong fluctuating input currents. It also connected iPRC structure to aspects of cellular physiology that could be estimated in simple current-clamp experiments. This allowed us to directly compare the iPRC predicted by our theory with the iPRC estimated from the response to EPSPs or current pulses in individual cells. We found that theoretically predicted iPRCs agreed well with estimates derived from synaptic stimuli, but not with those estimated from the response to somatic current injection. The difference between synaptic currents and those applied experimentally at the soma may arise from differences in the dynamics of charge redistribution on the dendrites and axon. Ultimately, our approach allowed us to identify novel ways in which voltage-dependent conductances interact with AHP conductances to influence synaptic integration that will apply to a wide range of cell types.

- phase response curve
- subthalamic nucleus
- basal ganglia

infinitesimal phase response curves (iPRCs) offer a simple representation of an autonomously oscillating neuron's response to synaptic input (Galán 2009; Gutkin et al. 2005; Rinzel and Ermentrout 1998; Winfree 2001), giving the rate at which an input current changes the neuron's oscillation phase as a function of the phase at which the input is delivered. Despite their simplicity, iPRCs can be used to predict some aspects of the collective behavior of neuronal populations, including their propensity for generating synchronized activity (Ermentrout 1996; Hansel et al. 1995; van Vreeswijk et al. 1994). This is particularly valuable for subthalamic neurons, which are embedded in a network of oscillating neurons (Bevan et al. 2002b) in which the degree of synchronous activity correlates with Parkinsonian pathology (Rivlin-Etzion et al. 2010). iPRCs can be measured experimentally or, given a realistic biophysical model, can be derived mathematically. Experimental measurement has the advantage of not requiring detailed knowledge of the cell's biophysical mechanisms but does not provide an understanding of the iPRC structure and may not even provide an accurate estimate of the iPRC, depending on the size and nature of the stimulus used. Mathematical derivation typically relies on analyzing the impact of inputs that are so weak that their effects are manifest only when repeated over many oscillation cycles, using singular perturbative techniques (Ermentrout 1981; Hoppensteadt and Izhikevich 1997; Neu 1979a, 1979b; Schwemmer and Lewis 2012). Although this kind of perturbation analysis requires a biophysical model, it does not directly reveal how biophysical parameters control the shape of the iPRC, nor does it disclose the range of stimuli over which the iPRC remains a reasonably accurate representation of a cell's response to input.

Subthalamic nucleus (STN) neurons are autonomously active, and their response to excitatory synaptic input appears to be well described by experimentally measured iPRCs (in the companion paper to this article: Farries and Wilson 2012). We developed a biophysical model of STN neurons to help us understand the structure of their iPRCs and to see why iPRCs appear to work so well as a predictor of their response to synaptic input. We sought to develop a method for predicting the iPRC of STN neurons that makes manifest the relationship between the cell's complement of ionic conductances and the structure of its iPRC; we compared the results of this effort with the iPRC derived from perturbation methods. We explored the range of stimuli under which the iPRC can accurately predict the model's response to input and examined sources of error in this prediction. We devised an experimental approach that allowed us to test theoretically predicted iPRCs against iPRCs estimated from the response to excitatory postsynaptic potentials (EPSPs) or injected current. We used this approach to address a number of questions raised by our data on experimentally estimated iPRCs, including the divergence between estimates made with synaptic and current pulse stimuli and the great diversity iPRC structure we encountered (Farries and Wilson 2012). This work has applications beyond the study of autonomously oscillating neurons and their response to weak perturbations; it provides new insights into how intrinsic properties shape a neuron's response to synaptic input.

## MATERIALS AND METHODS

All experimental procedures followed National Institutes of Health guidelines and were approved by the Institutional Animal Care and Use Committee of the University of Texas San Antonio. The methods for slice preparation, perforated-patch recording, and data analysis are described in the companion article (Farries and Wilson 2012).

#### Subthalamic neuron model.

We constructed a single-compartment model of a subthalamic neuron that qualitatively reproduced the behavior of STN cells recorded in brain slices. Aside from the spike-generating conductances (inactivating Na^{+} channels and repolarizing K^{+} channels) and leak conductances, five major classes of conductances have been identified in STN cells. Persistent Na^{+} conductances, which are available at subthreshold potentials (Beurrier et al. 2000; Bevan and Wilson 1999; Do and Bean 2003) and inactivate very slowly (Barraza et al. 2009; Do and Bean 2003), provide the major drive for autonomous pacemaking activity in these cells (Beurrier et al. 2000; Bevan and Wilson 1999; Do and Bean 2003). The small-conductance Ca^{2+}-dependent K^{+} (SK) current is a major source for the spike afterhyperpolarization (AHP) in subthalamic neurons (Bevan and Wilson 1999; Hallworth et al. 2003). STN cells express high-voltage-activated (HVA) Ca^{2+} channels (Bevan and Wilson 1999; Hallworth et al. 2003; Song et al. 2000) that are evidently the primary source for the Ca^{2+} that activates SK currents (Bevan and Wilson 1999; Hallworth et al. 2003). Finally, subthalamic neurons express hyperpolarization-activated, cyclic nucleotide-gated (HCN) channels and T-type Ca^{2+} channels (Atherton et al. 2010; Beurrier et al. 1999; Bevan and Wilson 1999; Nakanishi et al. 1987). We omit HCN and T-type Ca^{2+} currents from our model, as they are largely deactivated or inactivated, respectively, over the voltage range encountered during autonomous oscillations or during excitatory synaptic input (Atherton et al. 2010; Bevan et al. 2002a). This model differs from our past subthalamic model (Terman et al. 2002) in several ways, incorporating several new experimental findings on the nature of the subthalamic AHP (Hallworth et al. 2003; Teagarden et al. 2008), the overall current-voltage (*I-V*) curve (Farries et al. 2010), and the limited role of HCN and T currents in ongoing spontaneous activity (Atherton et al. 2010; Bevan et al. 2002a).

The model's current balance equation includes Na^{+}, K^{+}, Ca^{2+}, and leak conductances:
*I*_{inj} is an extrinsically applied current. The equilibrium potentials for the ion species are as follows: *E*_{Na} = +50 mV, *E*_{K} = −90 mV, *E*_{Ca} = +100 mV, and *E*_{leak} = −60 mV. The capacitance *C* was 100 pF, comparable to the effective cell capacitances estimated from the membrane potential change caused in STN cells by small current pulses used in our companion article (Farries and Wilson 2012; 109 ± 27 pF, range 57–191 pF, *n* = 66 cells; see below). Both inactivating (spiking) and persistent Na^{+} conductances were modeled with instantaneous activation; that is, their activation states were direct functions of voltage equal to their steady-state activation:
*h* is the Na^{+} inactivation gating variable, *m*_{Na}^{∞}(*V*) and *m*_{NaP}^{∞}(*V*) are the steady-state activation functions for the inactivating and persistent Na^{+} conductances, respectively, and *ḡ*_{Na} and *ḡ*_{NaP} are their maximal conductances. The K^{+} conductance was the sum of a delayed rectifier conductance (KDR) and SK:
*n* is the delayed rectifier activation variable, *m*_{SK} is the SK activation variable, and *ḡ*_{KDR} and *ḡ*_{SK} are their maximal conductances. Our model's Ca^{2+} conductance was noninactivating and active only at very depolarized potentials (i.e., HVA); like *g*_{Na}, its activation was treated as an instantaneous function of voltage:

The steady-state values for voltage-dependent gating variables were all sigmoidal functions of voltage,
*V*_{H} is the half-activation (or half-inactivation) voltage and *V*_{S} is the slope factor. The values of these parameters for each gating variable are given in Table 1. All dynamical variables except the membrane potential (Na^{+} inactivation *h*, KDR activation *n*, SK activation *m*_{SK}, and a Ca^{2+} sensor variable *s*, discussed below) obeyed differential equations of the form
*x*^{∞} is the steady-state value of the variable *x* and τ_{x} is its equilibration time constant. For voltage-dependent gating variables, the equilibration time constant had the form
_{base}, *V*_{H1}, *V*_{H2}, *V*_{S1}, and *V*_{S2} are given in Table 2.

Because we include the Ca^{2+}-dependent SK current, it is necessary to model Ca^{2+} dynamics in some way. If SK currents in STN cells tracked intracellular Ca^{2+} as measured, for example, by fluorescence of a Ca^{2+} indicator such as fura 2, it would be straightforward to design a simple model whose Ca^{2+} dynamics reproduced the behavior of the SK current in STN cells. However, STN SK currents do not behave in this way: although dependent on internal [Ca^{2+}], SK currents decay faster than measured Ca^{2+} transients, and SK currents do not accumulate with high rates of repetitive firing, whereas measured Ca^{2+} does accumulate (Teagarden et al. 2008). Indeed, SK channels in subthalamic neurons behave as if the Ca^{2+} driving their activation disappears very quickly so that SK currents cannot accumulate even at very high firing rates, and their decay is governed by the channel's own deactivation kinetics rather than by the decay of available Ca^{2+} (Teagarden et al. 2008). Such a situation can arise if Ca^{2+} entering the cell during an action potential cascades through a series of Ca^{2+}-binding molecules of varying Ca^{2+} affinity and kinetics: Ca^{2+} could initially interact with faster but lower affinity binding sites before being handed off to slower but higher affinity sites (Goldberg et al. 2009). This suggests that a model that correctly accounts for Ca^{2+} dynamics and SK activation in STN cells would be quite complicated. Rather than attempt to construct such a model, which would be very underdetermined by the data anyway, we have devised a simpler model linking Ca^{2+} currents to SK activation that reproduces three basic properties of SK currents in STN cells: *1*) they peak 10–15 ms after the spike; *2*) they decay with a time constant of about 30 ms irrespective of any Ca^{2+} dynamics; and *3*) they do not accumulate at high firing rates (Hallworth et al. 2003; Teagarden et al. 2008). To do this, we add one dynamical variable that represents the binding of Ca^{2+} to the SK Ca^{2+} sensor. The fraction of that sensor bound, which we call *s*, drives activation of the SK channel. It decays rapidly, representing the loss of Ca^{2+} to higher affinity sites (and preventing accumulation), but the rate of SK activation is saturated for most values of *s*, allowing SK activation to peak relatively late even though *s* decays rapidly (in real STN cells, SK activation is probably delayed by the action of lower affinity Ca^{2+} binding sites that are faster than the SK-linked sensor). The dynamics of *s* are governed by an on-rate, α_{s}, that depends on the Ca^{2+} current, *I*_{Ca}, and a fixed off-rate, β_{s}:
*s*-dependent on-rate, α_{SK}, and a fixed off-rate, β_{SK}:
^{2+} sensor binding *s* and SK activation *m*_{SK} are given by the usual relations:

Figure 1 illustrates the behavior of our SK model. The Ca^{2+} current flowing during a spike is sufficient to cause the level of sensor binding *s* to saturate; it reaches a plateau and then begins to decay after the action potential repolarizes (Fig. 1*A*, blue curve). The on-rate for the SK conductance (α_{SK}) is basically saturated for *s* > 0.15, and thus *g*_{SK} continues to rise even as *s* decays, achieving a maximum about 12 ms after spike initiation (Fig. 1*A*, red curve). Our modeled SK conductance shows negligible accumulation at high firing rates (Fig. 1*B*; compare to Fig. 8*B* of Teagarden et al. 2008).

#### Numerical estimation of iPRCs.

The iPRC of the STN model described above was computed using two techniques. The first technique employed the perturbative method to derive the adjoint equation, which was then solved numerically using XPPAUT, a freely available program for analyzing dynamical systems (Ermentrout 2002). This is often called the “adjoint method” for finding the iPRC. Our second method, known as the “direct method” (Achuthan et al. 2011), simulates the response of the model to small, brief current injections and measures the resulting phase shift. The iPRC at a given phase is the phase shift divided by the charge injected, as long as the injected current is small and brief enough to produce phase shifts that scale linearly with injected charge (Achuthan et al. 2011). We used Gaussian current waveforms with a peak value of 1 pA and widths (measured by standard deviation) of 0.01–0.03 ms (25–75 fC of total charge); we confirmed that phase shift scaled linearly over this range. The input phase was sampled at intervals no greater than 0.01 phase units, and the rapidly changing portions of the model's iPRC (φ < 0.05 and φ > 0.95) were sampled at much shorter intervals (0.0001–0.001). We define the phase of 1 (and 0) to be the point during the rising phase of the action potential when the rate of rise reaches its peak. The simulations were performed with Mathematica's (Wolfram Research) numerical differential equation solver using the LSODA method, with relative and absolute tolerances of 10^{−13} and 10^{−20}, respectively. These two methods for computing iPRCs produce virtually identical results once expressed in equivalent units. Our method for analyzing the biophysical basis of iPRC structure focuses on a segment of the periodic orbit that excludes the action potential itself (the spike spans φ = 0.998 to φ = 0.015, constituting <2% of the total period), beginning shortly after spike repolarization and ending at spike threshold. We define spike threshold in the model exactly as we do in our recordings: it is the voltage at which the second derivative of the voltage trajectory approaching a spike reaches 50 mV/ms^{2}.

#### Assessing the performance of phase models.

We evaluated the ability of iPRCs to predict the response of the full STN model by comparing their responses to current waveforms of two basic types. First, we used synaptic current waveforms to simulate the effect of EPSPs of varying amplitude. The functions specifying the synaptic current had the form
*t* = 0. τ_{rise} and τ_{decay} control the EPSC's onset and decay kinetics, respectively. We set τ_{rise} to 0.25 τ_{decay} unless otherwise noted. EPSP amplitude was defined as the integral of *I*_{syn} divided by the model capacitance (100 pF); *A* was adjusted to achieve the desired EPSP size. These synaptic current waveforms were used to generate normalized PRCs by simulating the response of the models to this input at 60 evenly spaced input phases (sampling φ = 0.03 to φ = 0.98 at intervals of 0.016); the resulting phase shifts were divided by the EPSP amplitude.

We also used noisy current waveforms to simulate the effect of synaptic bombardment by many presynaptic neurons. Noisy current waveforms were generated by convolving Gaussian white noise (generated at intervals of 0.4 ms) with a Gaussian kernel 2 ms in standard deviation; we used interpolation to sample the smoothed waveform at any interval desired. This waveform was scaled and shifted to obtain the desired mean and standard deviation. To compare iPRC predictions to the full model, we first drove the full model with noisy current until 1,000 spikes had been generated. We then took the current waveform during each interspike interval (ISI) and applied it to our phase models (see results for details of these models) to generate predicted ISIs. If a phase model failed to predict a spike before the end of the full model's ISI, a common occurrence, we continued the current injection using the waveform that spilled into the next ISI, until a spike was predicted. However, when considering the response of a phase model during the next ISI, we began with the current waveform as it was at the beginning of the full model's next ISI, rather than using a direct continuation of the current that completed the phase model's previous ISI. This method allows us to measure the prediction error that accumulated during individual ISIs without worrying about error accumulation across ISIs. Note that we do not force the full model to reset its state after each spike, because the full model is driven by continuous current waveforms of long duration without any explicit return to initial conditions. Thus the phase model prediction error includes error attributable to the (possibly erroneous) assumption that the state of the cell is the same after each spike in addition to the error that accumulates during the ISI.

#### Estimation of IV curves and AHP conductances.

We have developed a method for deriving a cell's iPRC from its instantaneous *I-V* relationship and its AHP conductances (treated as a function of time since the last spike). The instantaneous *I-V* curve describes the cell during rapid changes in membrane potential. Kinetically fast gating variables (e.g., Na^{+} activation) are treated as direct (instantaneous) functions of voltage, whereas slow gating variables (e.g., HCN activation) are treated as parameters (Borisyuk and Rinzel 2005). In our model, most non-AHP gating variables have very fast kinetics, with the exception of the Na^{+} inactivation variable. We nevertheless incorporated the voltage dependence of Na^{+} inactivation into the *I-V* curve (treating it as instantaneous), but with a modified voltage dependence (see results). Thus we could use a single instantaneous *I-V* curve to represent the cell's state throughout the ISI. To test this method experimentally, we estimated the instantaneous *I-V* curve and AHP conductances from current-clamp data. We did this by applying hyperpolarizing current pulses (−20 to −100 pA, 50–100 ms in duration) delivered late in the cell's autonomous ISI to reset the membrane potential after the AHP is largely complete; the subsequent voltage trajectory can be used to infer the *I-V* curve (see results). In some cases, this relatively long-duration current pulse was preceded by a shorter but stronger pulse (2 ms, −200 to −250 pA) to help abort spike generation when the pulse was initiated very late in the ISI. Current pulses were initiated at a fixed time during a trial, not at some particular time relative to the last spike, so current pulses initiated late in the ISI were obtained by selecting the subset of trials in which this had happened. Specifically, we selected trials in which *1*) the end of the current pulse occurred at least 80 ms after the time of the preceding spike (we used a longer criterion in most cases, typically 120–150 ms), and *2*) the membrane potential 5 ms after pulse termination was 0.1–10 mV below the minimum voltage achieved during the average spontaneous ISI. We aligned the traces from these selected trials to the time of the first spike following the current pulse and averaged them to get the voltage trajectory used to estimate the *I-V* curve.

The average spontaneous ISI voltage trajectory was computed by rescaling the time of individual ISI trajectories to make each one match the mean ISI duration. Each ISI trajectory began and ended with the point where the rate of rise during the action potential was maximal; this point corresponds to our definition of φ = 0/1, i.e., it defines the beginning and end of the oscillation cycle. The individual ISI voltage trajectories, all having identical durations after rescaling, were then averaged to produce the average spontaneous ISI voltage trajectory. The time derivative of this average trajectory was used to calculate the total current flowing during the ISI. The average voltage trajectory applied to the estimated *I-V* curve was used to calculate the current flowing through non-AHP conductances; the difference between the total current and the current predicted by the *I-V* curve yielded an estimate of the AHP currents. When the AHP current *I*_{AHP} was outward (hyperpolarizing), the AHP conductances could then be computed from
*V*(*t*) is the average voltage trajectory during the ISI and *E*_{K} is taken to be −90 mV.

The initial estimate of *g*_{AHP} always decayed to a nonzero plateau before diverging sharply upward or downward near the end of the ISI. We attribute both effects (the nonzero plateau and the late divergence) to errors in the initial estimate of the *I-V* curve (discussed in more detail in results). To obtain a better estimate of *g*_{AHP}, we identify a decaying phase for *g*_{AHP}, beginning 5–10 ms after the *g*_{AHP} peak and ending at a time well before the late divergence (typical decay phase duration was ∼100 ms). We fit this decay phase with an exponential, *a*·exp[−(*t* − *t*_{0})/τ] + *b*, where *t*_{0} is start time for the decay phase and time is measured relative to the beginning of the ISI. The fitted value of *b* gives the conductance level of the plateau; this is subtracted from the original estimate of *g*_{AHP}. In addition, *g*_{AHP} is set to *a*·exp[−(*t* − *t*_{0})/τ], the fitted function without the offset component *b*, for the region following the decay phase, thereby getting rid of the divergence late in the ISI. The transition from *g*_{AHP} determined by data to *g*_{AHP} defined by a decaying exponential function was smoothed by a weighted average of the two versions over a 10-ms window centered on the end of the decay phase. During this window, the weight assigned to the exponential fit rose linearly from 0 to 1 (and declined from 1 to 0 for the other component). Even after these adjustments to *g*_{AHP}, there were a few cases in which *g*_{AHP} was briefly negative early in the ISI. We attributed this to the resurgent Na^{+} current (Do and Bean 2003; Raman and Bean 2001) and accordingly assigned this portion of the AHP to a Na^{+} conductance; during this phase we assume that the *g*_{K} component of the AHP is zero, and outside of it we assume that the *g*_{Na} component is zero.

Once the final estimation of *g*_{AHP} is complete, the sum of the current predicted by this conductance and the current predicted from the original estimate of the *I-V* curve will no longer equal the total current flowing during the ISI, i.e., there will be a component of the total current that is not accounted for by either the AHP or the other intrinsic conductances. We attribute this current to errors in our original estimate of the *I-V* curve. We obtain a corrected *I-V* curve by subtracting the newly calculated AHP current from the total current and associating the value of that difference current at each time with the corresponding membrane potential drawn from the average ISI voltage trajectory, much as the original *I-V* curve was derived from the average voltage trajectory following hyperpolarizing current pulses and the derivative of that trajectory. However, the spontaneous ISI trajectory, unlike the trajectory following hyperpolarizing current pulses, often does not rise monotonically to spike threshold because of the action of the SK current. This means that there can be a range of voltages that are sampled twice, with slightly different associated current values. We solve this problem by smoothing the voltage-current data points with a Gaussian filter 0.1 mV wide.

#### Comparison of theoretical and empirical iPRCs.

Final estimates of the *I-V* curve and AHP conductance are used to obtain a theoretical prediction of the iPRC, which is then compared with the iPRC measured from the response to current injection and synaptic input. This comparison is complicated by the fact that the empirically measured iPRCs are not truly infinitesimal; they are finite-stimulus PRCs normalized by the size of the voltage change caused by the stimulus, and they approximate iPRCs only to the extent that the iPRC remains constant during the phase change occurring during the stimulus. For that reason, we often refer to these empirical iPRC estimates as “normalized PRCs.” To compare theoretical iPRCs with empirical normalized PRCs, we use the theoretical iPRC to predict the response to the stimulus current and divide the predicted phase shifts by the experimentally measured membrane potential change. This gives the theoretically predicted normalized PRC, in units of phase shift per millivolt, which can be directly compared with empirical normalized PRCs. In the case of synaptic input, the stimulus current waveform must be inferred from the membrane potential trajectory. The EPSP waveform is isolated by subtracting the average spontaneous voltage trajectory from the trajectory with EPSPs, as described in our companion article (Farries and Wilson 2012). The time derivative of this isolated EPSP voltage trajectory provides an estimate of the underlying synaptic current. Methods for measuring normalized PRCs with current pulses and synaptic input are described in our companion article, as are the methods for statistical analysis of PRCs.

Converting the time derivative of a voltage trajectory into a current, as we have done repeatedly, requires an estimate of the cell's capacitance. To obtain an estimate of the effective capacitance of the perisomatic region accessed by our recording electrode (which clearly includes more than just the somatic membrane but much less than the membrane of the entire cell), we used the response of the cell to the small, brief current pulses used to measure the PRC. To get the effective capacitance, we divided the integral of these current pulses (the total charge delivered) by average membrane potential change they caused (measured 5 ms after pulse termination). For the 10 cells that were tested with brief current pulses (the relatively long hyperpolarizing pulses used for *I-V* curve estimation were not suitable for this task), the effective capacitance was 101 ± 19 pF (range 76–128 pF). For the one remaining cell that did not receive suitably brief current pulses, we assumed a capacitance of 100 pF.

#### Simulation of a model with neurites.

We briefly examined the effect of charge redistribution in a spatially extended model, comparing the membrane potential in the axon initial segment following somatic current injection or dendritic synaptic conductance activation. This model was implemented in NEURON, a freely available neural simulation environment (Carnevale and Hines 2006). The dendritic tree was represented by a single cylinder 2.5 μm in diameter and 3 mm long, the soma was 15 μm in diameter and 15 μm long, and the axon was 1 μm in diameter and 1 mm long. The specific membrane capacitance was 1 μF/cm^{2}; the cytoplasmic resistivity was 100 Ω·cm. We were only interested in using this model to qualitatively study charge redistribution over time, not intrinsic oscillations or spike generation, so no intrinsic membrane conductances were included. In any case, zero membrane conductance is not an unreasonable approximation for subthalamic neurons during much of their autonomous oscillation, at least as measured from the soma (Farries et al. 2010). Somatic current injection was modeled as a 2-ms square pulse of +50 pA. The synaptic conductance was modeled by an alpha function (α = 1 ms) with a reversal potential of 0 mV and a peak conductance of 2 nS; the synapse was placed on the dendrite, 1.5 mm from the soma.

## RESULTS

#### Phase response curve of the STN cell model.

Our subthalamic neuron model reproduced the qualitative behavior we observed in subthalamic neurons recorded in brain slices. It fired spontaneously at 6.5 Hz and achieved a minimum membrane potential 27 ms after the preceding spike (Fig. 2*A*). When stimulated with a synaptic current waveform, it responded with a nondecaying potential change characteristic of STN EPSPs (Fig. 2*B*). The model's iPRC (Fig. 2*C*) was assessed with both small injected current waveforms and the adjoint method (Ermentrout 2002), producing identical results. Over most input phases, the model's iPRC exhibited a structure that was common in our sample of experimentally measured iPRCs: a type I PRC (i.e., one in which excitation always advances the time of the next spike) that is relatively insensitive to stimuli at early phases (but far from zero sensitivity for φ > 0.01), growing to be more than three times as sensitive at middle phases before declining again at later phases. However, the model's iPRCs diverged qualitatively from our experimentally measured iPRCs at very early (φ < 0.01) and very late (φ > 0.97) phases. At very late phases, the iPRC declined sharply to zero, a pattern we rarely saw in our experimentally measured iPRCs. At very early phases, the difference was more profound: the model iPRC exhibits a narrow type II region, that is, a range of input phases over which excitatory input would actually delay the next spike (Fig. 2*D*). Figure 2*D* plots the membrane potential trajectory (gray trace) with the iPRC on the same time scale, demonstrating that this type II region is confined to part of the repolarizing phase of the action potential, a property common in other models that are otherwise type I (Achuthan et al. 2011; Oprisan and Canavier 2002; Rinzel and Ermentrout 1998).

If iPRCs of real subthalamic neurons exhibit rapid changes at very early and very late input phases, like the ones exhibited by the model, such changes would be difficult to observe experimentally. This is primarily because the stimuli used to estimate iPRCs, especially EPSPs, are extended in time and thus will tend produce iPRC estimates that blur sharp changes in the true iPRC (noise in the measured iPRC resulting from ISI variability will also play a role). Of course, what really matters is the impact synaptic potentials have on spike times, not necessarily the fine structure of the true iPRC. Specifically, it is crucial to determine whether EPSPs could ever delay the next spike, a factor that has important implications for the establishment of synchronous activity in populations of interconnected neurons (Hansel et al. 1995). In STN cells recorded in brain slices, EPSP-associated currents arriving in the soma (inferred from the time derivative of the membrane potential) decayed with very short time constants (1.3 ± 0.8 ms, *n* = 89); the fastest current we measured decayed with a time constant of just 0.3 ms. This probably reflects an interaction with dendritic conductances (e.g., Wilson 1995) and does not directly represent the kinetics of a synaptic conductance. However, even the fastest synaptic current we ever measured was too slow to cause a spike delay in our model, and a current that was twice that fast produced a minuscule delay (Fig. 2*E*). Despite the presence of a nominally type II region in its iPRC, our model's response to any realistic EPSP is type I. Overall, the iPRC exhibited by our model resembled iPRCs commonly seen in experiments on real subthalamic neurons.

Rapid depolarizations in STN cells are followed by a transient drop in spike threshold (Farries et al. 2010). Because of this phenomenon, EPSPs delivered at late input phases should trigger spikes at a lowered threshold, an effect we observed experimentally. Our model also exhibits this phenomenon: when stimulated with an excitatory synaptic current waveform (same amplitude and time course as the one used in Fig. 2*B*), spikes following late phase inputs were triggered at a reduced threshold (Fig. 2*F*). As in real STN cells, this effect was strongest for spikes triggered during the EPSP (Fig. 2*F*, gray circles). However, this dynamic threshold phenomenon was quantitatively much weaker in the model than it was in real STN cells. This relative lack of threshold dynamics in our model probably explains its essentially zero second-order iPRC (not shown); experimentally observed second-order iPRC effects appeared to be due to persistent threshold changes (see companion article, Farries and Wilson 2012). In the model, spike threshold accommodation is caused by the relatively slow kinetics of the Na^{+} inactivation variable (*h*) so that more rapid depolarizations encounter a greater fraction of available (noninactivated) Na^{+} conductance at a given membrane potential. In cortical pyramidal neurons, on the other hand, it is due in large part to Kv1 K^{+} channels expressed in the axon (Higgs and Spain 2011). This may be the case in STN cells, as well, and it may allow much larger threshold shifts than Na^{+} inactivation dynamics alone. Such an effect cannot be represented directly in our single-compartment model, and this is one aspect of the STN's response to excitatory input that is not captured by our model.

#### Reduction of the STN model.

Most of the voltage-dependent currents expressed in STN cells that are active in the voltage range covered by its autonomous oscillation exhibit rapid kinetics and can be approximated by direct functions of the membrane potential. Indeed, this aspect is already partially incorporated into our model, because we treat the activation variables of Na^{+} and Ca^{2+} currents as instantaneous functions of voltage. Of course, a model whose gating variables are all direct functions of voltage cannot produce an action potential, but if we skip over the spike, treating it as a process that is initiated at a fixed threshold potential and taking a fixed (but small) amount of time before the cell is returned to a fixed resetting potential, we can reduce all voltage-dependent conductances to a single dynamical variable and use this reduced model to study the structure of the iPRC. We did this for our model, using the spike threshold measured during the autonomous oscillation (−36.2 mV), achieved 0.27 ms before the end of the cycle (φ = 0.998). We used the membrane potential at φ = 0.015 (2.29 ms after the start of the cycle) to define the resetting potential (−65.6 mV); the early part of the oscillation that is skipped by this choice is shown in Fig. 2*D*. We end up with a model that covers more than 98% of the oscillation, allotting 2.56 ms to the spike, which is ignored (Fig. 3*A*). We refer to the time from threshold crossing to the end of the cycle (0.27 ms) as *t*_{threshold} and the time from the beginning of the cycle to the voltage reset (2.29 ms) as *t*_{reset}.

Although many conductances can now be treated as direct functions of voltage, AHP conductances cannot. Those conductances consist of the K^{+} delayed rectifier, which is still in a process of rapid deactivation at our stating phase of 0.015, and the Ca^{2+}-dependent SK conductance, which is early in its activation at that point; Fig. 3*B* shows the values of these conductances during the ISI. Although the AHP conductances cannot be regarded as functions of voltage, they can be viewed purely as functions of time since the last spike. At early phases, the delayed rectifier is decaying toward its steady-state activation value, which is effectively zero for all *V* < −40 mV. Although the time constant of decay of KDR is voltage dependent, it is very fast (<1 ms) near the resetting potential, causing it to disappear quickly in a stereotyped fashion (i.e., it would do so even if the membrane potential were perturbed by an extrinsic input at early phases, unless that input was extraordinarily powerful and depolarizing). The SK conductance is governed by a Ca^{2+} sensor whose activity is essentially decoupled from the membrane potential except during the spike itself; it too acts as an autonomous function of time between action potentials. The SK conductance does decay slowly, raising the possibility that it could accumulate across ISIs if the cell was firing fast enough; that would mean that it could not be treated as a simple function of time since the immediately preceding spike. However, SK currents show very little accumulation in STN cells even at firing rates far higher than those used here (Teagarden et al. 2008), and we built that feature into our model (Fig. 1*B*). Thus AHP conductances in our model can be represented as a simple function of time (Fig. 3*B*).

Voltage-dependent gating variables are normally treated as direct functions of voltage by equating them to their steady-state values. However, this will not work for the Na^{+} inactivation (or availability) gating variable *h*. Figure 3*C* compares the steady-state value of *h* with its actual trajectory during the ISI. The largest deviation occurs early in the ISI, as the cell recovers from inactivation, but the error that is functionally significant comes late in the ISI, as the cell becomes depolarized enough to get significant Na^{+} activation; the slow kinetics of *h* cause a growing gap between *h*(*t*) and the steady-state value of *h*. We correct this error by equating *h* to a function of *V* that matches *h*(*t*) to *V*(*t*) over the region where *h*(*t*) is greater than its steady-state value (*t* > 52 ms, *V* > −66 mV) and sets *h* to its steady-state value otherwise; the resulting function *h**(*V*) is shown in the *inset* of Fig. 3*C*. Using this function for Na^{+} inactivation and equating all other gating variables to their steady-state values, we can express all non-AHP currents as direct functions of voltage (Fig. 3*D*). The resulting reduced STN model is described by a one-dimensional (but nonautonomous) differential equation in which membrane potential is the only remaining dynamical variable:
*f*(*V*) is the total membrane current due to voltage-dependent conductances (Fig. 3*D*, thick black curve). The ISI trajectory generated by the reduced model is almost identical to that of the full model (Fig. 3*E*), and its iPRC is very close to the full model iPRC over the phase range it covers (Fig. 3*F*).

#### Infinitesimal phase response curve of the reduced model.

To analyze the iPRC of the reduced model, we first consider the model without AHP conductances; the model then falls within the generalized integrate-and-fire paradigm considered by van Vreeswijk et al. (1994) and as a special case by Oprisan and Canavier (2002). Figure 4*A* shows the voltage trajectory generated by the model without AHP conductances; the period of oscillation is reduced from 152.8 to 84.2 ms. Because the only conductances left are direct functions of voltage, the membrane potential trajectory after a stimulus, and hence the duration of the remaining ISI, is completely determined by the voltage reached by the end of the stimulus. This also means that there is a one-to-one relationship between the model's membrane potential and its phase. If we consider the instantaneous delivery of charge at a certain point in the cycle, causing a sudden voltage jump Δ*V*, the resulting change in phase can be directly read off of the unperturbed trajectory (Fig. 4*A*). The magnitude of the phase shift is inversely proportional to the slope of the trajectory (Fig. 4*A*), and hence inversely proportional to *f*(*V*). The exact iPRC for this model can be derived as follows. Let *V*(*t*) be the unperturbed trajectory of this model (i.e., the solution to *Eq. 1* for *g*_{AHP} = 0), defined over *t*_{reset} ≤ *t* ≤ *T*_{1} − *t*_{threshold}, where *t* = 0 is the time of the spike starting the ISI and *T*_{1} is the period of oscillation without AHP currents. We denote the phase of the model without AHP as θ, defined as θ ≡ *t*/*T*_{1}. We thereby obtain the relationship between voltage and phase for this model, *V*_{θ}(θ) ≡ V(θ*T*_{1}) (this is the trajectory shown in Fig. 4*A*). The relationship between voltage and phase is invertible so long as *f*(*V*) < 0 for all *V* ∈ [*V*_{reset}, *V*_{threshold}], and *f*(*V*) must remain negative over that voltage range if the model is to generate autonomous oscillations. We denote the inverse of *V*_{θ}(θ) as θ_{V}(*V*), plotted in the *inset* of Fig. 4*B*.

The iPRC, *z*(θ), gives the rate at which an applied current *I*_{stim} changes the phase θ. Now the rate at which *I*_{stim} directly changes *V* is just *I*_{stim}/*C*. This can be related to the rate of change in phase by using the mapping from *V* to θ, θ_{V}(*V*):
_{V}(*V*) is derived from the voltage trajectory generated by the intrinsically expressed currents, the derivative of θ_{V} with respect to voltage is given by
*Eqs. 2* and *3*, and using the relationship between *V* and θ, we obtain an expression for the iPRC of the model without AHP currents:
*B*. This iPRC is valid for arbitrary stimulus currents; no assumptions about the magnitude of the current have been made. Thus AHP currents can be treated as an extrinsically applied current acting through this iPRC without any loss of accuracy.

Introducing the AHP into our phase model immediately brings a complication: the AHP current pulls the membrane potential below the minimum voltage achieved during the ISI without AHP (Fig. 5*A*), and thus *V* leaves the range over which the relationship between θ and *V* is defined. In the phase model, the AHP current, acting through the iPRC, would drive θ below zero. Because of the circular topology of the θ state space, that would push the cell through the spiking state (θ = 0/1) in reverse and leave it at a relatively late phase, ready to spike again at short latency. This problem is solved by changing the state-space topology of the phase model (Oh and Matveev 2009): we add a θ < 0 region that is accessed when θ = 0 is crossed from above (i.e., with θ̇ < 0; Fig. 5*B*, *inset*). The relationship between θ and *V* is defined for θ < 0 by integrating *Eq. 1* backward in time for *g*_{AHP} = 0. This is tantamount to computing the trajectory of the model without AHP currents for starting voltages below *V*_{reset}, but always defining the time at which *V*_{reset} is reached to be *t*_{reset}; the relationship to phase is obtained by using *t* = θ*T*_{1} (Fig. 5*B*). The iPRC is still given by *Eq. 4* (Fig. 5*C*).

We are now free to treat the AHP as an applied current acting through the iPRC. The phase of our reduced STN model, driven by a combination of *I*_{AHP} and an extrinsic current, *I*_{stim}, obeys
*k*_{1} ≡ 1/*T*_{1}.

For reasons that will quickly become clear, we have changed notation slightly and labeled the iPRC given by *Eq. 4* as *z*_{1}(θ); this iPRC is multiplied by −*I*_{AHP} because *z*_{1}(θ) assumes the sign convention used for injected currents. *Equation 5* does not represent the fact that *I*_{AHP} is voltage dependent via its dependence on the K^{+} driving force, *V* − *E*_{K}. However, because there is a one-to-one mapping between θ and *V*, we can express the driving force as a function of θ and use it to define an iPRC for the AHP conductance *g*_{AHP}: *z*_{AHP}(θ) ≡ −[*V*(θ) − *E*_{K}]*z*_{1}(θ). The phase model with AHP thus becomes
*I*_{stim} = 0 is longer than *T*_{1}; we denote this lengthened period *T*_{2} (it is the period of the original model, 152.8 ms). This is a perfectly good phase model of the reduced biophysical model, but it does not use the cell's natural phase variable, φ ≡ *t*/*T*_{2}, and the cell's iPRC in those natural coordinates is not given by *z*_{1}. From the point of view of an experimenter studying this cell, θ is a hidden variable and *z*_{1}(θ) cannot be found using conventional PRC measurement techniques. We need to derive the model's iPRC in natural phase coordinates, denoted *z*_{2}(φ), from *z*_{1}(θ) and *g*_{AHP}(*t*). The relationship between φ and θ can be obtained by solving *Eq. 6*: let θ(*t*) be the θ trajectory when *I*_{stim} = 0; then θ_{0}(φ) ≡ θ(φ*T*_{2}) gives θ as a function of φ in the absence of input (Fig. 5*D*).

One way in which the AHP alters the iPRC is simply a direct result of changing the oscillation period. Late in the oscillation period, when *I*_{AHP} ≈ 0, that is the only effect. If *I*_{AHP} is small enough to be neglected, then θ simply grows linearly with φ (Fig. 5*D*). In that case, a stimulus-induced phase shift, δθ, can be converted to a shift in the natural phase coordinate φ by using the original definitions of these variables in terms of time: dθ = *dt*/*T*_{1} and *d*φ = *dt*/*T*_{2}, so δφ = (*T*_{1}/*T*_{2})δθ. Under these conditions, *z*_{1}[θ_{0}(φ)], which gives the change in θ per unit charge delivered as a function of φ (Fig. 5*E*), can be converted into a change in φ simply via multiplication by the factor *T*_{1}/*T*_{2}:
*I*_{AHP} ≈ 0 for the remainder of the cycle. This is not true if *I*_{AHP} ≠ 0 for any part of the remaining cycle, because *z*_{1}[θ_{0}(φ)] only gives the immediate stimulus-driven change in θ; that initial change δθ will be further altered by the action of the AHP as it interacts with the iPRC. Specifically, the initial change δθ shifts θ relative to *g*_{AHP}(*t*) so that *g*_{AHP}(*t*) now acts on a slightly different part of the iPRC, *z*_{AHP}(θ). In this way, the direct change caused by the stimulus is supplemented with indirect changes that accumulate as a phase-shifted *g*_{AHP}(*t*) acts through *z*_{AHP}(θ). If we know the net change in θ caused by the stimulus, combining direct and indirect effects, it can be converted to a change δφ with the factor *T*_{1}/*T*_{2}. In other words, if we know how δθ(*t*) evolves over time following stimulation, we can obtain δφ by multiplying δθ evaluated at the end of the cycle [≈δθ(*T*_{2}) if the stimulus is small] by *T*_{1}/*T*_{2}.

To obtain δθ(*t*), we consider the impact of a very brief stimulus current that is effectively a delta function centered on *t* = *t*_{stim}, delivering a small amount of charge *q*_{stim}. The result, derived in the appendix, is
*t*_{stim}/*T*_{2} is the stimulus input phase. Note that the solution is a product of the direct stimulus-driven change in θ, *q*_{stim}*z*_{1}[θ_{0}(φ)], and a time-dependent exponential factor that gives the additional effects arising from the shift in *g*_{AHP}(*t*) relative to *z*_{AHP}(θ). If this factor is evaluated at the end of the oscillation period (*t* = *T*_{2}), then it becomes a function of the stimulus input phase φ, which we now denote *g*(φ):
*t*_{stim} = φ*T*_{2}. *g*(φ) is shown in the *inset* of Fig. 5*E*. Multiplying *g*(φ), *z*_{1}[θ_{0}(φ)], and the factor *T*_{1}/*T*_{2} gives us an explicit expression for the iPRC of the model with AHP (*z*_{2}) in terms of the iPRC without AHP (*z*_{1}) and the AHP conductance *g*_{AHP}(*t*):

The result is shown in Fig. 5*F* (green curve); it agrees quite well with the iPRC of the full STN model (Fig. 5*F*, thick black curve). Note that the method we used to transform the iPRC for injected current, *z*_{1}(θ), into an iPRC for AHP conductances could just as easily be used to derive iPRCs for synaptic conductances (Achuthan et al. 2011). Although this may be advantageous for describing the response to synapses located near the soma, the data in our companion article (Farries and Wilson 2012) show very little dependence of EPSP amplitude on somatic membrane potential, presumably due to the relatively distal dendritic location of the activated synapses (Rall 1967), and thus we continue to treat synaptic input as an injected current.

Although *Eq. 8* looks formidably complex, the effects of the AHP on the iPRC fall into four intuitively intelligible categories. First, the AHP lowers the membrane resistance, reducing the voltage change caused by a given amount of current. Second, by lengthening the ISI, the AHP reduces the phase shift at all input phases because a given change in spike time Δ*t* [predicted, for example, from the instantaneous *I-V* curve *f*(*V*)] constitutes a smaller fraction of the overall oscillation period; this effect is captured by the factor *T*_{1}/*T*_{2}. Third, AHP controls what part of the *I-V* curve the cell occupies at a given phase, and hence which part of the *I-V* curve-defined iPRC *z*_{1}(θ) is presented when the stimulus arrives; this effect is expressed by *z*_{1}[θ_{0}(φ)], plotted in Fig. 5*E*. Finally, there is the altered interaction between the AHP current and *z*_{1}(θ) that occurs when a stimulus arrives, expressed by the function *g*(φ). The exact effect of this interaction is given by *Eq. 7*, but the role of the derivative of iPRC in this expression can be understood qualitatively, treating the AHP as if it were an injected current for the sake of clarity (replacing *z*_{AHP} with −*z*_{1} and *g*_{AHP} with *I*_{AHP} in *Eq. 7*). The part of the iPRC *z*_{1} that the AHP current acts on is determined by the phase trajectory θ(*t*). When a stimulus is delivered, the entire subsequent trajectory θ(*t*) is shifted relative to θ_{0}(*t*), which causes the time-dependent AHP current to act on a slightly different part of the iPRC. If the iPRC is constant over phase, i.e., d*z*_{1}/dθ = 0, then this does not matter; the value of *z*_{1} multiplied by the remaining AHP current is unchanged [*g*(φ) = 1]. If, on the other hand, the iPRC rises with increasing θ (d*z*_{1}/dθ > 0), then positive phase shifts (for example) cause the remaining AHP current to act on a more sensitive region of *z*_{1} and thus lengthen the ISI more than they would otherwise, counteracting the phase advance [*g*(φ) < 1]. If d*z*_{1}/dθ < 0, then phase advances push *I*_{AHP}(*t*) into a region of reduced sensitivity, lessening the effect of the AHP and thereby amplifying the effect of the initial phase advance [*g*(φ) > 1].

#### Accuracy of the phase model representation.

Unlike the derivation of *z*_{1}(θ), which is valid for arbitrarily large stimulus currents, we invoked the smallness of the stimulus current at several points in our derivation of *z*_{2}(φ) (see appendix). It is therefore important to determine how well *z*_{2}(φ), and the other iPRCs, predict the response of our full STN model to finite stimuli. Figure 6*A* shows normalized PRCs of the full model generated with simulated EPSPs of varying amplitude. The shape of the normalized PRC does vary with the size of the stimulus used to measure it, but most of this does not reflect deviation from the predictions made by the phase model representation; a similar variation with EPSP size is evident when the phase model itself (with the iPRC shown in Fig. 2*C*) is used to make normalized PRCs (Fig. 6*A*, *inset*). Instead, this variation reflects the fact that larger EPSPs act on extended regions of the iPRC: phase advances caused by early parts of the synaptic current cause later parts of the current to act on later sections of the iPRC. That is not to say that large EPSPs do not cause deviations from the phase model prediction. Figure 6*B* compares the normalized PRCs generated by the full model and by the full model's iPRC for small (0.5 mV) and large (8 mV) EPSPs; the phase model prediction differs slightly from that of the full model, and this difference is greater for the larger EPSP.

We evaluate the error in normalized PRCs for three phase models. The first simply uses the iPRC of the full model (Fig. 2*C*); this is the model used to generate the “iPRC prediction” of Fig. 6, *A* and *B*. The second uses the iPRC of the reduced model, that is, the derived iPRC *z*_{2}(φ) given by *Eq. 8*. The third uses the *I-V* curve-based iPRC *z*_{1}(θ) with the AHP treated as an extrinsically applied K^{+} conductance (“iPRC with AHP”), i.e., the model defined by *Eq. 6*. The “iPRC error” for each of these models is defined as the normalized phase shift they predict minus the normalized shift generated in the full model, using synaptic current waveforms of varying size. Figure 6*C* plots this error as a function of input phase for 8-mV EPSPs. All phase models underestimate the phase shift generated in the full model, albeit slightly. The performance of the reduced model iPRC is comparable to the full model iPRC, but the model with an applied AHP conductance performed markedly better at all but late input phases (φ < 0.8). When the error magnitude is averaged over the entire normalized PRC and plotted as a function of EPSP size (Fig. 6*D*), we see that the average error grows more rapidly for the reduced and full model PRCs than it does for the applied AHP model, although the error magnitude is comparable across models for small EPSPs (<1 mV). Curiously, the error associated with the applied AHP model actually drops with increasing EPSP size up to EPSPs of ∼3 mV before rising very slowly with EPSP size for larger EPSPs; we do not know why.

If phase models are to be useful in understanding neuronal population dynamics in vivo, they should produce at least qualitatively correct behavior when driven by noisy current waveforms that approximate the effect of synaptic bombardment a neuron might experience in an intact brain. We tested our models with noisy current waveforms of varying mean and standard deviation. As with the normalized PRCs, the performance of three phase models was evaluated by comparison with the behavior of the full STN model. In this case, the spike train generated by the full model (Fig. 6*E*, *top* trace) was compared with the spike trains predicted by the phase models (Fig. 6*E*, spike rasters) when driven by the same current waveform (Fig. 6*E*, *bottom* trace). Histograms of the spike time prediction error of the three phase models relative to the full model spike times are shown in Fig. 6*F*. Each column shows the results of using current waveforms of increasing standard deviation, producing ISI distributions with an increasing coefficient of variation (CV). As CV increases, the distribution of prediction errors both broadens and becomes more highly skewed toward extreme values for the reduced and full iPRC models. The median error magnitude grows roughly linearly with CV for these models (Fig. 6*G*), and although it remains modest even as the CV approaches 1 (∼8 ms at a CV of 0.94), the 90th percentile error grows much more rapidly (from ∼1.6 ms at a CV of 0.05 to ∼70 ms at a CV of 0.94). In contrast, the applied AHP model continues to perform well as the stimulus current becomes more variable (Fig. 6*G*, blue circles), and even the 90th percentile error remains under 7 ms for the current waveforms we tested. Altering the firing rate by adjusting the mean current level had much less effect on spike prediction error for all models (Fig. 6*H*).

Comparing the performance of these models allows us to identify distinct sources of error that arise from the use of phase models and to determine their relative magnitudes. Under most circumstances, the phase model with an applied AHP was more accurate than models that incorporate the effect of the AHP into their iPRC. Because the applied AHP model is isomorphic to the reduced STN model described by *Eq. 1*, its error can be attributed to the assumptions of that reduction, namely, that the gating variables are direct functions of voltage, that the spike is triggered at a fixed voltage threshold, and that the AHP conductance is just a function of time since the last spike. The magnitude of this error remains quite small over a wide range of input regimes (Fig. 6, *D–H*) and is largely encountered at late input phases (Fig. 6*C*), presumably because that is where the neglect of Na^{+} availability (*h*) dynamics takes its toll. The additional error associated with the reduced model iPRC stems from neglect of the changing interaction between the AHP and the cell's voltage-dependent conductances (embodied in the *I-V* curve-derived iPRC *z*_{1}). This error grows with increasing input variability as the assumption of small phase shifts used to derive *Eq. 8* becomes increasingly invalid, and it is by far the largest source of error when the model is driven by noisy current waveforms (Fig. 6*G*). However, this error is not very large when the model is driven by individual EPSPs (Fig. 6, *B–D*). Because the blurring effect of finite-sized EPSPs on normalized PRCs is relatively minor for EPSPs less than 4 mV in amplitude (Fig. 6*A*), this suggests that iPRCs estimates derived from the response to individual EPSPs, normalized by EPSP size, could be close to the true iPRC.

#### Experimental estimates of the instantaneous I-V curve and AHP conductance.

We have devised a single compartment biophysical model that behaves qualitatively like subthalamic neurons recorded in brain slices (Fig. 2), and we have presented a theoretical approach that allows us to understand how the conductances present in that model control the structure of its iPRC (Figs. 4 and 5). We would like to use *Eq. 8* to predict the iPRCs of real subthalamic neurons and test it against their empirical iPRCs estimated from the response to synaptic input or somatic current injection. Use of *Eq. 8* requires knowledge of the cell's *I-V* curve and its AHP conductances. If we hope to explain the structure of a neuron's specific iPRC, which can vary considerably from cell to cell, we need the *I-V* curve and AHP conductances specifically of that neuron, rather than their averages derived from separate experiments on different cells. We can estimate these things from a simple current-clamp experiment, and we did so in 11 STN cells recorded using the perforated-patch technique.

Figure 7 illustrates our technique for estimating the functions used by *Eq. 8* from experimental data in an example cell. We seek to model the cell's behavior during the ISI, starting from some postspike resetting point we select (2 ms into the cycle in this case, Fig. 7*A*, *inset*) to the point where spike threshold is crossed. If the cell's intrinsic conductances in this region are truly just a combination of purely time-dependent AHP conductances that eventually decay to zero and non-AHP conductances that are well approximated by direct functions of voltage, then the *I-V* curve can be derived from the time derivative of the voltage trajectory after the AHP is over. A hyperpolarizing current pulse initiated near the end of an ISI can bring the membrane potential below the minimum achieved during the unperturbed ISI; once the pulse is terminated, the cell passes through the full voltage range experienced during spontaneous ISIs, but in the absence of AHP currents (Fig. 7*A*). We take a segment of this trajectory, starting well after the end of the current pulse but still below the minimum potential achieved during unperturbed ISIs (Fig. 7*A*, green circle) and ending at spike threshold (red circle), and use its time derivative as a function of voltage to form an initial estimate of the instantaneous *I-V* curve (this also requires an estimate of the cell's capacitance).

This initial estimate of the *I-V* curve (Fig. 7*B*, main curve) was used to estimate the current flowing through intrinsic voltage-dependent conductances (but not AHP conductances) during spontaneous ISIs (Fig. 7*B*, *inset*, red curve). Because the total current flowing into the perisomatic region accessed by our recording electrode can be determined from the time derivative of the spontaneous ISI (Fig. 7*B*, *inset*, black curve), the *I-V* curve-predicted current can be subtracted from the total current to obtain an estimate of the AHP current (Fig. 7*B*, *inset*, blue curve). If we assume that this difference current is generated by a K^{+} conductance with an equilibrium potential of −90 mV, we can compute the time-dependent conductance that generated it. The resulting conductance estimate (Fig. 7*C*, gray curve) closely resembles what we would expect from a combination of spike-repolarizing K^{+} conductances and the SK conductance (cf. Fig. 3*B*), but it appears to decay to a nonzero plateau before becoming sharply negative just before the spike. Both of these effects are probably due to conductance changes caused by the hyperpolarizing pulse. For example, the *I-V* curve-predicted current likely includes a small amount of inward HCN current that is not present during the spontaneous ISI, which must be compensated by a small outward component in the difference current. We correct this (and the late negative conductance) by assuming that AHP conductances must decay to zero asymptotically (even if they do not necessarily reach zero within the unperturbed ISI). This correction entails discarding a constant offset component of the conductance estimated from the difference current as well as the rapid changes occurring at the end of the ISI; what is left is our estimate of *g*_{AHP}(*t*) (Fig. 7*C*, black curve). The estimated *g*_{AHP}(*t*) always decayed with a time constant consistent with the behavior of SK conductances in STN neurons (32.9 ± 13.9 ms, range 15.9–58.8 ms).

Our correction of the AHP conductance leaves a part of the total current unaccounted for by either the *I-V* curve of Fig. 7*B* or *g*_{AHP}(*t*). We assume that this remaining current reflects the difference between the initial estimate of the *I-V* curve derived from the voltage trajectory following hyperpolarizing current pulses and the real *I-V* curve operating during unperturbed ISIs. We use this current, now treated as a function of voltage rather than of time, to correct our initial estimate of the *I-V* curve and produce our final estimate of membrane current as a function of voltage, *f*(*V*) (Fig. 7*D*). At hyperpolarized potentials, adjustment of the original estimate entails removal of a small amount of inward current (Fig. 7*D*, *inset*) corresponding to the HCN conductance activated by the hyperpolarizing current pulse; this was true of all cells subjected to this analysis (*n* = 11). At depolarized potentials, the corrected *I-V* curve is shifted to the left, predicting a larger inward current than the original estimate at a given voltage. In this respect, the example cell shown in Fig. 7 is somewhat unusual; in many of our cells, the corrected *I-V* curve was shifted to the right, possibly corresponding to compensation for T-type Ca^{2+} conductance deinactivatived by hyperpolarizing current. With the corrected *I-V* curve, we can produce a model of the cell without AHP currents (Fig. 7*E*). Starting from our selected resetting time and potential (2 ms, −61.5 mV; Fig. 7*A*, *inset*), we integrate *f*(*V*) forward in time until spike threshold (−39.4 mV) is reached and backward in time to the minimum voltage achieved during the unperturbed ISI (−64.4 mV), thereby obtaining the shortened AHP-free oscillation period *T*_{1} (70.4 ms) and voltage as a function of the AHP-free phase variable θ (Fig. 7*E*, *inset*). We can then apply *Eq. 4* to get *z*_{1}(θ) (Fig. 7*E*) and, finally, use *Eq. 8* to get the theoretically predicted iPRC *z*_{2}(φ) (Fig. 7*F*). We also obtained iPRC estimates in the same cell using synaptic and current pulse stimuli; these normalized PRCs were converted to the same units as *z*_{2}(φ) via an estimate of the effective cell capacitance (Fig. 7*F*).

The procedure outlined above was applied to our entire sample of 11 cells largely as described, but with one important exception. In two of these cells, the difference current used to estimate the AHP conductance exhibited a brief period of inward current early in the ISI, wedged between an initial phase dominated by a rapidly decaying outward current and a late phase that begins with the rise of the outward SK current. In both cases this inward current region lasted ∼4 ms and, after subtraction of the outward current plateau required to make *g*_{AHP}(*t*) approach an asymptotic value of zero, achieved peak values of −5.3 and −9.5 pA. Ascribing the entire AHP to a K^{+} conductance would require *g*_{AHP}(*t*) to be negative during this inward current phase; clearly, the AHP conductance of these cells has additional components. In fact, such a component has already been identified in subthalamic neurons: the resurgent Na^{+} current (Do and Bean 2003). The resurgent Na^{+} current is thought to result from a distinct inactivation state of Na^{+} channels, called the “open blocked” state, that can be exited only through the open channel state; this causes Na^{+} channels in the open blocked state to admit current briefly after an action potential as they recover from inactivation and then deactivate (Raman and Bean 2001). We therefore attributed the brief inward component AHP current, when it appears, to a Na^{+} conductance; *z*_{2}(φ) is still calculated as before, except that *g*_{AHP} and *z*_{AHP} are broken into distinct K^{+} and Na^{+} components. It is likely that the early part of the AHP comprises a mixture of K^{+} and Na^{+} conductances regardless of whether the net AHP current is inward or outward. However, in the absence of a method for determining what this mixture is, we simply assume that the AHP conductance is purely K^{+} when its current is outward and is purely Na^{+} when inward.

#### Comparison of predicted and measured iPRCs.

One way to compare normalized synaptic or current pulse PRCs with the theoretically predicted iPRC is to multiply their phase shift per millivolt by an estimate of the cell's effective capacitance to obtain the phase shift per unit charge delivered, as we did for the example shown in Fig. 7. However, this does not account for the way normalized PRCs generated by finite stimuli depend on stimulus amplitude even in pure phase models (Fig. 6*A*, *inset*). Although this blurring effect by finite stimuli is small, we chose to avoid it by using the theoretical iPRC and an estimate of the stimulus current to generate the theoretically predicted phase shift for a given stimulus and then normalize by stimulus amplitude. Figure 8*A* shows the average theoretically predicted normalized PRC generated by synaptic current in our sample of 11 cells, together with the actual normalized synaptic PRCs. The average synaptic and theoretical PRCs agree quite well over most input phases but diverge at late input phases. Specifically, the difference between theoretical and synaptic PRCs is not statistically significant at the *P* = 0.05 level for any φ < 0.68 but remains significant at that level for all φ ≥ 0.68 (Fig. 8*B*, green circles). The difference between synaptic and theoretically predicted PRCs is never significant at the Bonferroni-corrected level of *P* = 0.001, but if differences averaged over early (φ < 0.7) and late phases (φ > 0.7) are considered, we find that the phase-averaged difference is statistically significant at late phases (*t*-test, *P* = 0.007) but not at early phases (*P* = 0.19). In 9 of 11 cells, we also were able to measure PRCs using injected current pulses. Normalized current pulse PRCs exhibited markedly less stimulus sensitivity than both theoretical and synaptic PRCs over most input phases (Fig. 8*A*). When normalized current pulse PRCs are compared with theoretical predictions generated by square current waveforms of matching amplitude and duration, we find that the difference between them is significant at the *P* = 0.05 level for all φ < 0.72 (Fig. 8*B*, black circles) and meets the Bonferroni-corrected criterion (*P* < 0.001) over most of that range (0.22–0.60).

Although theoretically predicted and synaptically measured iPRCs largely agree on average, the iPRCs predicted by *Eq. 8* may not necessarily account for the cell-by-cell variation in synaptic iPRC structure we reported in our companion study (Farries and Wilson 2012). In Fig. 8, *C* and *D*, we examine two cells whose synaptic iPRCs differ qualitatively: one is most sensitive to inputs delivered at early phases with a gradual decline in sensitivity as phase increases (Fig. 8*C*), whereas the other exhibits an increase in sensitivity with increasing phase until a peak is reached at φ ≈ 0.5 (Fig. 8*D*). In both cases, the theoretically predicted normalized PRCs exhibit the same qualitative pattern as the normalized synaptic PRCs. This suggests that our theoretical PRCs can do more than just predict the average iPRC and gives us an opportunity to identify the specific biophysical differences that allow two cells of the same physiological type to exhibit such distinctive patterns of responsiveness to synaptic input. The crucial difference between the cells shown in Fig. 8, *C* and *D*, is found in their *I-V* curves. The *I-V* curve of the cell shown in Fig. 8*C* has a negative slope over the entire subthreshold voltage range encountered during the ISI, whereas the *I-V* curve of the other cell (Fig. 8*D*) has a weakly positive slope at relatively hyperpolarized potentials (Fig. 8*E*) before becoming negative at more depolarized potentials. The slope of the *I-V* curve *f*(*V*) controls the derivative of the AHP-free iPRC *z*_{1}(θ), and from *Eq. 4* (and the fact that *V* is a monotonically increasing function of θ), one can see that a voltage range over which d*f*/d*V* > 0 implies a θ range over which d*z*_{1}/dθ > 0 (Fig. 8*E*, *inset*). The AHP magnifies this difference because its ability to reduce the cell's sensitivity to input is amplified when d*z*_{1}/dθ > 0 and curtailed when d*z*_{1}/dθ < 0. Thus a comparatively minor change in the balance of negative-slope (e.g., persistent Na^{+}) and positive-slope (e.g., leak or HCN) conductances can have a profound effect on the cell's sensitivity to synaptic input at different phases.

#### Why do iPRCs estimated from synaptic and current pulse stimuli differ?

In our companion study (Farries and Wilson 2012), we found that synaptic input is much more effective than current injection at changing phase for all but late input phases, even after normalization by stimulus size (illustrated in Fig. 8*A*). Our finding that predicted iPRCs more closely match synaptic iPRCs than they do current pulse iPRCs suggests that somatically injected current may introduce some kind of artifact that distorts the measured iPRC. The problem may lie in the way somatically injected charge redistributes itself within the cell over time. The initial membrane potential change Δ*V* caused by current injection partially decayed over the first few milliseconds following the end of the pulse (Fig. 6*B* of Farries and Wilson 2012), whereas the Δ*V* caused by EPSPs usually did not decay at all. We attribute the partial decay of Δ*V* following current injection to the gradual redistribution of charge from the site of injection to more distal parts of the neuron. This decay appeared to be largely complete 5 ms after the end of the pulse, so we measured Δ*V* at that point to normalize current pulse iPRCs. However, in dendrite-bearing neurons, the process of charge equalization does not take place with a single time constant (Rall 1969), and given STN neurons' unusual membrane properties, namely, a slope conductance close to zero (Farries et al. 2010), some of these time constants could be quite large. Thus the redistribution of somatically injected charge could continue throughout most or all of the ISI.

If the issue of charge redistribution was just a question of how to measure Δ*V* to correctly normalize the current pulse PRC, then such PRCs would have the same shape as synaptic PRCs and the difference between them would be proportionately the same at all phases. Yet, this is not the case. In subthalamic neurons, spikes are initiated in the axon initial segment (Atherton et al. 2008), and thus stimuli are effective at advancing spike time only to the extent that they depolarize that part of the cell. Because of gradual charge redistribution into the dendrites, current injected into the soma can initially have a large immediate impact on the membrane potential in the axon initial segment that decays over time. We illustrate this in a model with a soma, dendrite, and axon (Fig. 8*F*, *inset*) but zero membrane conductance. Current injected into the soma initially causes a large depolarization in the axon, but that change decays as charge flows into the dendrite; even after 20 ms, a gradual loss of axonal Δ*V* continues (Fig. 8*F*, black trace) as charge trickles into yet more distal regions of the dendrite. A current pulse can cause relatively more phase shift at late input phases because spike threshold will be reached before as much charge has escaped into the dendrites. The same is not true of synaptic input, because it charges up the dendritic membrane before depolarizing the axosomatic region (Fig. 8*F*, green trace). The additional charge used to depolarize the dendrites does not directly contribute to a phase shift (and is never seen by our somatic recording electrode), but it does prevent the charge that does make it to the axosomatic region from flowing back into the dendrites, allowing even early phase EPSPs to be fully effective at altering spike timing. We emphasize that this is just a plausible hypothesis explaining the differences between synaptic and current pulse PRCs; a more certain explanation would require further research.

## DISCUSSION

The iPRC is generally regarded as a valid description of a neuron's input-output relationship only in the limit of weak inputs, and its mathematical derivation typically invokes perturbations so feeble that their effects are manifest only after summation over many oscillation cycles. Our data on the response of subthalamic neurons to EPSPs (see our companion article, Farries and Wilson 2012) suggest that an iPRC might adequately represent the response even to relatively strong synaptic inputs, at least under our recording conditions. This suggestion is corroborated by the behavior of our STN model, whose response to powerful fluctuating current injection was reasonably well predicted by its iPRC. It appears that the behavior of subthalamic neurons can be represented by a very simple one-dimensional model over a significant input range. This indicates that phase models are likely to be useful tools for understanding the STN, and perhaps for understanding the basal ganglia as a whole.

Part of the reason for the success of the iPRC description is the fast kinetics of many of the STN's voltage-dependent conductances. This allows the current flowing through these conductances to be approximated by a direct function of voltage. A model governed by such an *I-V* curve and otherwise behaving like an integrate-and-fire neuron (fixed spike threshold crossing followed by a return to a fixed resetting voltage) is a “pure phase model” (van Vreeswijk et al. 1994); i.e., it is isomorphic to a phase model without any assumptions about the strength of the input. Subthalamic neurons are not integrate-and-fire neurons, but their biophysical mechanisms can be split into one component that does act like a generalized integrate-and-fire neuron and a second component, the AHP, whose contribution can be understood from its action through the phase model defined by the first component. This approach allows us to explain why the iPRC description works as well as it does (many conductances have fast kinetics) and to identify the major source of deviation from the iPRC description (input-induced shifts in the interaction between AHP currents and voltage-dependent conductances). It also allows us to understand, in a transparent way, how the neuron's biophysical properties shape its iPRC.

The correspondence, in real STN neurons, between iPRCs predicted by our theory and iPRCs estimated from the response to EPSPs demonstrates that our approach can account for much of the response of subthalamic neurons to synaptic input. On the other hand, the large and consistent divergence between iPRCs estimated using current pulses and both synaptic and theoretical iPRCs suggests that the current injection method of iPRC measurement is flawed. We hypothesize that this flaw arises from the differing patterns of charge distribution generated by the two kinds of stimuli. In this view, somatically injected charge gradually spreads into the distal dendrites and away from the spike trigger zone, whereas synaptically delivered charge largely remains once it arrives in the perisomatic region, because the synaptic conductance will have already charged the dendritic membrane. The nature and size of the error introduced by this flaw must depend on the spatial distribution of active conductances and membrane potential during the oscillation, and it may be smaller when charge loss across the perisomatic membrane is large compared with loss into the distal dendrites; perhaps this explains why this error was not observed in cortical pyramidal neurons (Reyes and Fetz 1993). In any case, the error in the STN is evidently quite large, and any PRC measurement technique that employs somatic current injection, including sophisticated methods like dynamic clamp, should probably be validated in the cell type being studied by comparison with the response to the synaptic input of interest.

In addition to correctly predicting the average iPRC, our theory is able to account for the variability in iPRC structure we observed experimentally. Without knowing why iPRCs might be so variable, it was difficult to be certain how much of the experimentally observed variability was mere measurement error. Our analysis demonstrates how subtle biophysical changes can have a profound effect on the iPRC, thereby bolstering confidence in the experimental findings. There is, however, one aspect of experimentally measured iPRCs that is not accounted for by our theory: the continued sensitivity to inputs at late phases. The theory predicts that the iPRC should always decline toward zero at late input phases, due to the growing negative slope in the *I-V* curve as spike threshold is approached. It is not clear why many iPRCs estimated using EPSPs or current pulses fail to show this expected decline. Because this part of the iPRC is estimated from trials in which a spike was often triggered before the stimulus was complete (see companion paper, Farries and Wilson 2012), the measured iPRCs at late phases are probably less reliable than at earlier phases and may include some systematic error that explains the late phase divergence from the predicted iPRC. Alternatively, this divergence could be caused or exacerbated by the dynamic threshold phenomenon, which is expected to increase sensitivity to depolarizing inputs at late input phases; our theory assumes a static spike threshold.

Although our theoretical analysis was developed to explain the structure of subthalamic iPRCs, it has considerably greater generality than might at first appear. Our analysis can explain the iPRC structure of any neuron whose intrinsic properties are well approximated by an instantaneous *I-V* curve and time-dependent AHP conductances. Not all cells have AHPs that so closely resemble functions of time since the last spike, but it is a reasonable approximation in many cases. Furthermore, the theory can accommodate voltage-dependent AHP conductances, as long as that voltage-dependence is kinetically fast. As for the *I-V* curve, many voltage-dependent conductances are too slow to be treated as direct functions of voltage. However, some of those are slow enough to be treated as approximately constant over the ISI; the dynamics of such conductances would be expressed by changes in the iPRC over time. Conductances too fast to be dealt with in that way might be approximated by modified functions of voltage derived from their average behavior during ISIs, as we did in the STN for the Na^{+} availability variable *h*. Application of our theory does not even require that the neuron be autonomously active. There is no fundamental distinction between intrinsic and extrinsic currents in our approach, so the mean level of the extrinsic conductance (or current) that drives activity in the neuron can be incorporated into its *I-V* curve; the resulting iPRC would predict how input fluctuations about that mean alter the neuron's input-driven spiking pattern.

There will inevitably be situations in which a one-dimensional phase model fails to capture important aspects of a neuron's response to its inputs. Indeed, this is likely to occur in the STN. For example, prolonged hyperpolarization of subthalamic neurons will engage HCN and T-type Ca^{2+} currents (Nakanishi et al. 1987); HCN activation might still be incorporated into a dynamic iPRC, but T channel inactivation might not. Even within the inhibition-free regime we have considered, the iPRC description probably fails for sharp increases in excitatory drive that trigger large drops in spike threshold; spike threshold dynamics can have a potent effect on the response of subthalamic neurons to excitation (Farries et al. 2010). Our STN model exhibits much less spike threshold accommodation than real STN cells, and so it probably overstates the ability of phase models to predict spike times under strong fluctuating input currents. It may be possible to deal with some of these issues without entirely abandoning the phase description by adding additional dynamical variables that are coupled to phase.

The reduction of the full complement of a neuron's biophysical mechanisms to a phase model is not just a specialized technique for analyzing autonomously oscillating neurons coupled by weak synapses, nor is it just a method of model reduction. Phase is an explicit representation of when the next spike is expected to happen, and thus the iPRC directly manifests the relationship between synaptic input and spiking output. If one can connect biophysics to iPRC structure in a transparent way, as we have done, then phase models become a tool for understanding how a cell's intrinsic physiological properties shape the connection between input and output. In this regard, phase reductions are likely to provide valuable insights into neural function even when they cannot predict spike times from input current with high accuracy. One such insight is the rather counterintuitive revelation that depolarization-activated depolarizing currents, such as the persistent Na^{+} current, do not boost the response to EPSPs as measured by changes in spike time. On the contrary, they suppress responsiveness by reducing the amount of phase shift elicited by a given amount of synaptically delivered charge. Another achievement that is specific to our approach is the identification of novel mechanisms whereby the AHP can influence synaptic integration. This, too, offers some unexpected results. Although a hyperpolarizing AHP reduces input sensitivity overall, locally it can either boost or suppress sensitivity depending on the mix of other conductances it interacts with. As these examples illustrate, analysis of phase models can deepen our understanding of the way a neuron's intrinsic physiological properties shape its response to synaptic input.

## GRANTS

This work was supported by National Institute of Neurological Disorders and Stroke Grant NS-47085.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## AUTHOR CONTRIBUTIONS

M.A.F. and C.J.W. conception and design of research; M.A.F. performed experiments; M.A.F. analyzed data; M.A.F. interpreted results of experiments; M.A.F. prepared figures; M.A.F. drafted manuscript; M.A.F. and C.J.W. edited and revised manuscript; M.A.F. and C.J.W. approved final version of manuscript.

## APPENDIX

To obtain δθ(*t*), we consider the impact of a very brief stimulus current that is effectively a delta function centered on *t* = *t*_{stim}, with *I*_{stim}(*t*) = 0 for all *t* ≠ *t*_{stim}, delivering a small amount of charge, *q*_{stim}. Let θ_{p}(*t*) be the θ trajectory when perturbed by *I*_{stim}(*t*) and θ_{0}(*t*) be the trajectory in the absence of any stimulus; δθ ≡ θ_{p} − θ_{0} is the difference between these trajectories. Both θ_{p}(*t*) and θ_{0}(*t*) must obey *Eq. 5*:
*t*):
*t*) between the perturbed and unperturbed trajectories is zero for *t* < *t*_{stim}; we are only interested in *t* > *t*_{stim}. Because *I*_{stim}(*t*) = 0 for *t* > *t*_{stim}, *I*_{stim}(*t*) merely sets the initial condition of θ_{p} at *t* = *t*_{stim}. Using the fact that *q*_{stim} is small, we get
_{p} − θ_{0} is also small, which means that the difference in *z*_{AHP} evaluated along the two trajectories can be approximated by δθ multiplied by the derivative of *z*_{AHP} with respect to θ, evaluated along the unperturbed trajectory:
*Eq. A1*, we obtain a relatively simple differential equation for the difference δθ between the perturbed and unperturbed trajectories for *t* > *t*_{stim}:
*z*_{1}(θ) is treated as a function of time because it is evaluated along the trajectory θ_{0}(*t*). δθ(*t*) is just the net change in θ caused by the stimulus, from both direct and indirect effects, that has accumulated up to time *t*; the total change in θ can be found by evaluating this function at the end of the cycle: δθ(*T*_{2}).

- Copyright © 2012 the American Physiological Society