## Abstract

Many neurons exhibit subthreshold membrane-potential resonances, such that the largest voltage responses occur at preferred stimulation frequencies. Because subthreshold resonances are known to influence the rhythmic activity at the network level, it is vital to understand how they affect spike generation on the single-cell level. We therefore investigated both resonant and nonresonant neurons of rat entorhinal cortex. A minimal resonate-and-fire type model based on measured physiological parameters captures fundamental properties of neuronal firing statistics surprisingly well and helps to shed light on the mechanisms that shape spike patterns: *1*) subthreshold resonance together with a spike-induced reset of subthreshold oscillations leads to spike clustering and *2*) spike-induced dynamics influence the fine structure of interspike interval (ISI) distributions and are responsible for ISI correlations appearing at higher firing rates (≥3 Hz). Both mechanisms are likely to account for the specific discharge characteristics of various cell types.

## INTRODUCTION

An ongoing challenge in neuroscience is to relate large-scale network phenomena to the properties of single neurons, including their subthreshold dynamics. It has been shown that subthreshold membrane-potential resonances influence a given brain area's ability to sustain rhythmic activity (Buzsáki and Draguhn 2004; Geisler et al. 2005), which, in turn, underlies various aspects of perception, memory, and behavior (Buzsáki 2005; Lisman 2005). Recent findings for grid cells in layer II of the entorhinal cortex (EC) even suggest that there may be a direct link between the subthreshold dynamics of single neurons and their representation of the sensory environment (Giocomo et al. 2007). Despite the increasing evidence for the interplay of single-cell and network dynamics, it is not fully understood how subthreshold frequency selectivity, as reflected in membrane-potential oscillations and subthreshold resonance, shape spike-discharge patterns in individual neurons. These discharges, however, form the link between single-cell and network activity and are therefore subject to investigation in our study.

Subthreshold resonances are usually examined by injecting an oscillatory current into a neuron, while blocking its synaptic receptors. Intracellular recordings have shown that the voltage response of most neurons, such as EC pyramidal cells, decreases with increasing frequency of the injected current. The cell membrane of these neurons, often referred to as nonresonant, therefore acts as a low-pass filter. In contrast, resonant neurons, such as EC stellate cells (Erchova et al. 2004; Haas and White 2002), have maximum response amplitudes for stimuli delivered at higher frequencies. In the last two decades many examples of resonant neurons have been found, for example, in the inferior olive (Llinás and Yarom 1986), thalamus (Hutcheon et al. 1994), hippocampus (Leung and Yu 1998), and cortex (Gutfreund et al. 1995). These resonant neurons often show subthreshold membrane potential oscillations whose frequencies are close to the resonance frequency.

Subthreshold resonance arises from the activation of delayed membrane currents (as reviewed in Hutcheon and Yarom 2000) and has been shown to influence spiking neural activity by affecting firing rate, precision of spike timing (Desmaisons et al. 1999; Haas and White 2002; Schaefer et al. 2006), and spike clustering (Chen and Shepherd 1997; Desmaisons et al. 1999; Izhikevich et al. 2003; Pedroarena et al. 1999; Wu et al. 2001). Despite these successes, the exact relationship between subthreshold resonance and spike patterns has not yet been established. The current study aims to fill this gap.

Our data underscore the close relation between subthreshold and spiking dynamics, demonstrating that fundamental characteristics of the spike-train statistics in individual neurons can be deduced from their subthreshold resonance properties. Using a resonate-and-fire model with parameters directly constrained by experiments, we quantitatively reproduce the experimentally measured distribution of interspike intervals (ISIs) and their correlations. Importantly, the model provides a straightforward mechanistic explanation for the transition from subthreshold dynamics to spike patterns in both resonant and nonresonant neurons and shows that a spike-induced reset of the phase of subthreshold oscillations accounts for the observed spike clustering in resonant neurons.

## METHODS

### Physiological data

Experiments, data acquisition, and data analysis were carried out as described in Erchova et al. (2004) and Schreiber et al. (2004a). For the present study, we recorded from 19 neurons (7 stellate, 12 pyramidal) in rat superficial entorhinal cortex. Each neuron was injected with hyper- and depolarizing current pulses lasting for 500 ms, 10 s, 20 s, and 30 s to characterize the cell type and to obtain ISI distributions. To study resonance properties, cells were injected with an oscillatory “ZAP” current of 30-s duration (1) whose oscillation frequency increased in time, *f*(*t*) = *f*_{m}*t*/2*T* (Gimbarzevsky et al. 1984; Puil et al. 1986). Here, *I*_{1} denotes the stimulus amplitude, *I*_{0} is the injected constant current, *T* is the stimulus duration, and *f*_{m} = 20 Hz is the maximum stimulus frequency.

The true membrane potential of all recorded cells was verified at the end of each recording session. The cell-voltage offset was measured after the electrode was gently pulled out of the cell and cleaned from the attached membrane. All recordings were then corrected by the value of this offset.

### Resonate-and-fire model

##### SUBTHRESHOLD REGIME.

To describe subthreshold dynamics we use an extension of an electrical circuit model that has previously been successfully used to capture the frequency-dependent subthreshold dynamics in stellate and pyramidal cells of the EC (Erchova et al. 2004; Schreiber et al. 2004a). In the circuit interpretation, the original model consists of two parallel branches (Fig. 1*A*). The first branch is characterized by a conductance *g* in parallel with a capacitance *C* and mimics the integrative properties of a leaky integrator. The second branch consists of a resistance *R*_{L} in series with an inductance *L* and captures the response properties of delayed rectifying currents.

In the present study, we model the dynamics driven by cell-intrinsic noise. Because synaptic transmission was blocked, the dominant intrinsic noise source was channel noise (Manwani and Koch 1999; White et al. 2000). Accordingly, the conductance becomes a time-dependent quantity *g*(*t*) whose dynamics we describe by an additional stochastic differential equation (Fig. 1*A*). This equation captures the dynamics of a population of ion channels that open and close independently and with constant rates. Although the stochastic conductance *g*(*t*) creates multiplicative noise, the full three-dimensional model can be reduced to a two-dimensional dynamical system with an additive noise term ξ(*t*). This simplification (derived in Verechtchaguina et al. 2007) is possible because the correlation time of conductance fluctuations (≪1 ms; see White et al. 1998) is much smaller than the dominant timescale of the subthreshold voltage dynamics (10–100 ms). The membrane potential *V*(*t*) thus evolves according to (2) Here, *I*(*t*) represents the time-dependent input current and ξ(*t*) is white Gaussian noise (〈ξ(*t*′)〉_{t′} = 0, 〈ξ(*t* + *t*′)ξ(*t*′)〉_{t′} = δ(*t*)) where δ(*t*) is Dirac's function.

In the model, each cell is characterized by its membrane capacity *C*, resting membrane potential *V*_{r}, membrane input resistance *R*_{0}, and intrinsic noise intensity *D*. Resonance properties are described by two additional parameters, γ and δ. In terms of damped harmonic oscillators, γ/*C* represents the effective damping coefficient and (1/2π) is the effective eigenfrequency. The relationship between the parameters of the circuit model and the parameters *V*_{r}, *R*_{0}, γ, and δ used here is summarized in Fig. 1*A*.

In the noise-free model (*D* = 0), stimulation with a short current pulse results in a damped subthreshold oscillation when γ < 2 (underdamped regime). The frequency of the resulting oscillation is *f*_{osc} = Ω/2π = (1/2π)/*C* and the amplitude decays exponentially with a time constant *t*_{rel} = 2*C*/γ. For γ > 2 (overdamped regime), the voltage *V* relaxes monotonically to the resting potential. Driven by intrinsic noise (*D* > 0), the model generates sustained subthreshold oscillations with fluctuating amplitude and phase in the underdamped regime (γ < 2, whereas in the overdamped regime (γ > 2) the subthreshold voltage fluctuations do not exhibit rhythmic components.

To relate the resonance properties of the model to experimental data, the theoretical frequency-dependent impedance function |*Z*|(*f*) was calculated for the noise-free case (i.e., *Eq. 2* with *D* = 0). This simplification mimics the experimental situation where impedances were calculated from voltage responses averaged over at least five trials, and results in (3) In the underdamped regime the impedance function has a maximum at the frequency *f*_{res} (4) For small γ, the resonance frequency *f*_{res} is similar to the effective eigenfrequency and the frequency *f*_{osc} of the deterministic subthreshold oscillation, with *f*_{res} slightly exceeding *f*_{osc}. Recordings from stellate cells of the entorhinal cortex (Erchova et al. 2004) confirm this prediction of our overall model framework. In the overdamped regime the impedance function |*Z*_{theory}|(*f*) decays monotonically with increasing frequency (i.e., the cell is nonresonant). Nevertheless, it may have a nonzero eigenfrequency (Erchova et al. 2004).

In the model, injecting a prolonged constant current *I*(*t*) = *I*_{0} shifts the stationary voltage to a new value *V*_{r} + *I*_{0}*R*_{0}. Since we were interested only in the stationary responses we introduced an auxiliary variable *x*(*t*) = *V*(*t*) − *V*_{r} − *I*_{0}*R*_{0}, which describes fluctuations around a zero-mean voltage (5)

Already in its noise-free form, this simple phenomenological model reproduces both resonant and nonresonant subthreshold dynamics (Erchova et al. 2004; Izhikevich 2001; Mauro et al. 1970; Richardson et al. 2003; Verechtchaguina et al. 2006b) and can be related to Hodgkin–Huxley-type conductance-based models via linearization (Mauro et al. 1970; Villacorta and Panetsos 2005). Because the model parameters depend in general on the level of membrane depolarization, experimental data for each value of the input current *I*_{0} were treated separately, similar to the approach taken by Erchova et al. (2004). The range of parameters was then estimated for each *I*_{0} value (as described in the following text).

##### SPIKING REGIME.

###### Renewal model.

For a minimal description of spike generation, the model needs to be equipped with a firing threshold and a voltage reset. Slightly extending this approach to incorporate refractoriness, we introduced three parameters describing spike generation: voltage threshold *x*_{b}, voltage reset *x*_{0}, and reset time τ_{r}. Whenever *x*(*t*) reaches *x*_{b} from below, a spike is generated, and *x*(*t*) is reset to *x*_{0}, d*x*(*t*)/d*t* to zero. After the reset period τ_{r}, the time evolution of *x* again follows *Eq. 5*. For *D* > 0, this dynamics generates a sequence of independent spikes and is referred to as the *renewal* resonate-and-fire model. Two examples of ISI distributions generated with this model are shown in Fig. 1*B*. Analyzed in detail in results, one of the model predictions is that the distance between the first and the second peaks in the ISI density is equal to the period of subthreshold oscillations. The inverse of this peak-to-peak distance is therefore approximately equal, but slightly smaller than the subthreshold resonance frequency (see explanation following *Eq. 4*). The change in peak-to-peak distance with the subthreshold resonance and oscillation frequencies is illustrated in Fig. 1*B*.

###### Model extension: nonrenewal process.

To cover ISI correlations, deterministic spike generation was replaced by a stochastic rule. This means that a spike is not necessarily generated when the membrane potential crosses the voltage threshold. Rather, spike generation fails with a probability *p*_{s}(*t*) (6) Here, *t*_{i} represents the spike times and the summation is taken over all spikes that have been generated before time *t*. If the expression under the sum exceeds 1, then *p*_{s}(*t*) = 1 and no spike can be generated.

The probability *p*_{s}(*t*) is the key new dynamic variable. Initially, *p*_{s}(*t*) is set to zero. When triggered by a single spike, *p*_{s}(*t*) instantaneously increases by an amount *p*_{0}, remains roughly at this level during the correlation time τ_{c}, and then decays back to zero on a timescale 2/β. The probability *p*_{s}(*t*) is evaluated each time the voltage variable *x*(*t*) crosses the threshold *x*_{b} from below; with probability 1 − *p*_{s}(*t*) a spike is generated. If no spike is elicited, *x*(*t*) evolves without reset according to *Eq. 5* until it again hits the threshold from below. If several spikes are generated within a time window smaller than the correlation time, the probability *p*_{s}(*t*) exhibits complex dynamics as illustrated in Fig. 6. As the example suggests, the stochastic threshold mechanism limits the number of spikes within a cluster and causes “silent overshoots” (allowing the membrane potential to briefly exceed the threshold value without producing a spike). The resulting *nonrenewal* resonate-and-fire model generates spike trains whose ISIs are negatively correlated over short timescales. The prediction that the inverse distance between the first and second peaks in the ISI density is approximately equal to the subthreshold resonance frequency also holds for the nonrenewal model.

##### PARAMETER ESTIMATION.

###### Subthreshold dynamics.

The model parameters in *Eq. 5* were estimated from impedance measurements and autocorrelation analysis. The frequency-resolved impedance was determined from the voltage response *V*(*t*) to the injected ZAP current *I*(*t*) (*Eq. 1*) as *Z*(*f*) = FFT[*V*(*t*)]/FFT[*I*(*t*)], where FFT is the fast Fourier transform algorithm. Using the analytical expression for |*Z*_{theory}|(f) (*Eq. 3*), the values of the parameters *C*, δ, γ, and *R*_{0} were obtained from least-square fits of |*Z*_{theory}|(*f*) to the experimental curve *Z*(*f*). For two examples, see Fig. 2, *A* and *B*. For larger ZAP amplitudes and at more depolarized levels of membrane potential occasional spikes were observed in the ZAP response. For impedance estimation, spikes were removed by cutting voltages at −55 mV, which did not change the impedance estimates. The estimates of resonance frequency presented in Fig. 8 are obtained as averages of resonance frequencies across different levels of membrane depolarization (based on model fits to the data). Given that resonance frequencies in stellate cells are relatively constant with depolarization (Erchova et al. 2004), this procedure is justified.

The intrinsic noise intensity *D* was estimated by comparing the experimental and theoretical autocorrelation function *r*_{x}(*t*) = 〈*x*(*t* + *t*′)*x*(*t*′)〉_{t′}. The experimental autocorrelation function was calculated from voltage responses to hyper- and depolarizing current steps, for time intervals that did not contain spikes. A Butterworth low-pass filter (cutoff frequency 40 Hz) was applied to the data before the parameter estimation. The autocorrelation function of the model can be obtained analytically. In the underdamped regime (γ < 2 it reads (7) where Ω = /*C*. Accordingly, *r*_{x}(*t*) oscillates with the period *T*_{p} = 2π/Ω and its amplitude decays with the relaxation time *t*_{rel} = 2*C*/γ (see Fig. 2*C*).

Using *Eq. 7*, the values of *r*_{x}(*t*) at *t* = 0 and at the first minimum *t*_{1}, as well as the value of *C* (from the impedance estimate), we calculated parameter values for stellate cells via the following relations (8) In the overdamped regime (when γ > 2), the autocorrelation function *r*_{x}(*t*) is given by a variant of *Eq. 7*, where the trigonometric functions are replaced by their hyperbolic counterparts. Consequently, *r*_{x}(*t*) decays monotonically (Fig. 2*D*).

All parameters were estimated separately for each steady-state depolarization. Note that parameters do not vary independently, but obey an additional constraint, *σ _{x}^{2}* =

*D*/(γδ), imposed by the model. For resonant cells γ and δ were estimated in two ways, using either impedance or autocorrelation functions. The values obtained by both methods were almost identical.

###### Spike dynamics.

The parameters of the renewal model describing spike generations—i.e., voltage threshold *V*_{b}, reset value *V*_{0}, and reset time τ_{r}—were directly estimated from the experimental data. To determine the value of the spike threshold *V*_{b}, we first estimated the SD σ_{VD} of the voltage derivative in the absence of spikes. Second, we defined a derivative threshold for spike initiation as 3σ_{VD}. Third, we defined the local voltage threshold *V*_{b} as that point where the voltage derivative crosses the derivative threshold. The procedure was performed for each individual spike. The local reset value *V*_{0} was set to the minimal voltage within 25 ms after the spike, whereas the local reset duration τ_{r} was defined as the time between threshold crossing and reaching of the reset value. All parameters were first evaluated for each spike and then averaged over all spikes generated at a given value of membrane depolarization (not less than 40 spikes). The steady-state voltage was estimated by the mean membrane potential in time intervals without spikes. The model threshold *x*_{b} and voltage reset value *x*_{0} were obtained by subtraction of the steady-state voltage from *V*_{b} and *V*_{0}, respectively.

### Spike-train statistics

##### ISI DISTRIBUTIONS.

We denoted successive spike times by *t*_{1}, *t*_{2}, … , *t*_{i}, … , *t*_{N+1}, where *i* is the position of a spike in the measured spike train. The *i*th ISI was defined as *X*_{i} = *t*_{i+1} − *t*_{i}.

A rough estimate for the ISI probability density ℱ(*T*) can be obtained from binned experimental data. The inverse distance between the first and second peaks in the ISI densities in Fig. 8 is based on distances averaged across ISI densities obtained for two to three different values of injected constant current. To adequately localize the second peak, only ISI distributions based on firing rates of ≥4 Hz were considered. This criterion was satisfied by six of seven stellate cells. The ISI densities were sampled with a bin width of 15 ms and peaks were identified at bins with the largest amplitude.

Clustering in stellate cells was reflected in multimodal ISI distributions, as shown in Fig. 2*G*. For this study, we therefore defined clustered spike patterns as those patterns where the corresponding ISI distribution showed a significant minimum between its first two peaks. Specifically, we required a confidence level of ≥95%, based on a likelihood-ratio test. Spikes separated by a distance smaller than the distance at this minimum were considered as belonging to a cluster; typical distances between spike clusters corresponded to the location of the second and potentially third peaks of the ISI density. Note that our criterion of multimodal ISI distributions is stronger than what is theoretically needed to discriminate between clustered and unclustered spike trains; however, all stellate cells with sufficient statistics satisfied this criterion. No pyramidal cells showed spike clustering, as additionally confirmed by visual inspection.

Due to the limited sample size (150–600 spikes per spike train) the exact shape of estimated ISI densities depends strongly on the bin size. For model fitting, we therefore used the cumulative ISI distribution *P*(*T*) = *∫ _{0}^{T}* ℱ(

*t*)d

*t*. Unlike ℱ(

*T*),

*P*(

*T*) can be reliably estimated from limited data as (9) where Θ(

*T*) is the Heaviside step function [Θ(

*T*) = 0 if

*T*< 0 and Θ(

*T*) = 1 if

*T*> 0].

##### ISI CORRELATIONS.

Because the same ISI distribution could potentially be generated by a renewal or a nonrenewal process, we used the following measures to test whether ISIs were independent:

*1*) The joint probability density ℱ(*T*_{1}, *T*_{2}) of two successive ISIs (i.e., the density of the ISI return map). ℱ(*T*_{1}, *T*_{2})d*T*_{1}d*T*_{2} gives the probability that successive ISIs satisfy *X*_{i} ∈ (*T*_{1}, *T*_{1} + d*T*_{1}) and *X*_{i+1} ∈ (*T*_{2}, *T*_{2} + d*T*_{2}). The joint densities ℱ(*T*_{1}, *T*_{2}) were estimated by histograms and normalized to unity.

*2*) Correlation maps *C*(*T*_{1}, *T*_{2}), characterizing the deviation of the experimental ISI sequence from the renewal ISI sequence (with independent ISIs). Correlation maps *C*(*T*_{1}, *T*_{2}) were obtained from measured joint ISI densities by subtracting the joint ISI density expected for the corresponding renewal spike train (10) Positive values of *C*(*T*_{1}, *T*_{2}) indicate that this ISI pair occurs more frequently than under the renewal assumption; negative values indicate pairs that are less frequent.

*3*) Serial correlation coefficients (SCCs) ρ_{k} describe correlations between arbitrarily distant ISIs (11) Here, averages are taken over all ISIs in a spike train and the index *k* is called a lag. For a renewal spike train ρ_{k} = 0 for all *k* > 0. To estimate the typical correlation time, we plot the SCCs as a function of time, by linearly scaling the lag axis with the mean ISI.

##### MODEL OPTIMIZATION.

To find a set of model parameters reproducing the experimental ISI distribution, we estimated the permitted parameter range from the experimental data (as described earlier). The required parameters were found through a simulated annealing (SA) technique (Press et al. 1999) described in the supplemental material.1

The membrane capacity *C* was held constant (*C* = 2.5 × 10^{−4} μF for the stellate cell and *C* = 1.7 × 10^{−4} μF for the pyramidal cell) because it typically did not depend on the membrane potential; γ, δ, and *D* were optimized within the permitted parameter range (see Fig. 3 and Table 1). Likewise, the threshold *x*_{b} and reset *x*_{0} varied within the measurement range. Since *V*_{r} depended sensitively on the threshold definition, we tested an extended range of thresholds to compensate for possible biases. To capture spike-induced effects potentially lasting beyond the time of maximum hyperpolarization, the reset time τ_{r} was allowed to vary up to the size of the minimal ISI observed in each specific spike train. The parameters for the nonrenewal threshold mechanism—τ_{c}, β, and *p*_{0}—were allowed to vary freely within the range indicated in Table 1. Note that optimal parameter values for the nonrenewal model differed from the renewal model (Table 1). The nonrenewal model admits “silent” threshold crossings without spike generation. To obtain similar firing rates, the total number of threshold crossings was therefore larger for the nonrenewal model.

To quantify the discrepancy between the measured ISI distribution *P*_{X}(*T*) and the theoretical distribution *P*_{0}(*T*), we used the Kolmogorov–Smirnov (KS) test (Ledermann 1984). The test hypothesis was that the theoretical distribution *P*_{0}(*T*) is the distribution underlying the measured ISIs. The test statistics 𝒟 in the KS test represents the largest difference between two distributions (12)

For a significance level α = 5% and a given sample size *N*, the critical value 𝒟_{K} is calculated as 𝒟_{K} = 1.36/. If 𝒟 exceeds the critical value 𝒟_{K}, the hypothesis has to be rejected. To graphically illustrate the KS test, a tube with radius 𝒟_{K} was drawn around the curve corresponding to the experimental ISI distribution *P*_{X}(*T*). If the theoretical ISI distribution *P*_{0}(*T*) laid completely inside this tube, it was regarded as a potential distribution underlying the measured ISIs.

For the renewal model, the ISI distribution *P*_{0}(*T*) was obtained analytically using the Stratonovich approximation (Verechtchaguina et al. 2006a,b) described in the supplemental material. Because the model parameters for both resonant and nonresonant cells belong to a regime of moderate damping and moderate noise intensity, the Stratonovich approximation was applicable and accurate. For the nonrenewal model, *P*_{0}(*T*) was estimated on the basis of stochastic simulations consisting of ≥3,000 spikes.

The SA procedure is described in detail in the supplemental material. For the renewal model the procedure minimizes the KS statistics 𝒟. For the nonrenewal model the procedure simultaneously minimizes the KS statistics 𝒟 and the discrepancy between experimentally measured and theoretically predicted values of the first serial correlation coefficient ρ_{1}.

Note that SA belongs to a class of heuristic algorithms that do not necessary converge to the global optimal solution. Thus the method allows us to find a model generating a spike pattern similar to that observed experimentally, but does not guarantee that this model is unique.

## RESULTS

This study aims to elucidate the relation between subthreshold and spiking dynamics of single neurons. To this end, we performed intracellular recordings of stellate and pyramidal cells, two major cell classes in the entorhinal cortex. Constant and time-varying ZAP-current injections (see methods) were used to investigate spontaneous subthreshold membrane-potential oscillations and resonances, respectively. Furthermore, constant depolarizing currents were applied to move the mean membrane potential near threshold and thus elicit spike trains driven by cell-intrinsic noise. The data were interpreted within the framework of resonate-and-fire model neurons.

### Subthreshold dynamics

Stellate cells exhibited a pronounced membrane-potential resonance in the theta range (5–15 Hz), whereas pyramidal cells showed mostly none, or at most a very small resonance at lower frequencies, between 2 and 4 Hz (Fig. 2, *A* and *B*). In all resonant cells, small constant depolarizing currents caused voltage oscillations with a rhythmic component close to the resonance frequency (Fig. 2*C*). In contrast, the voltage fluctuations of nonresonant cells did not exhibit any rhythmicity (Fig. 2*D*).

### Firing patterns

For larger depolarizing currents all investigated neurons produced action potentials. The discharge patterns, however, differed significantly between both cell types. In particular, clustered spikes were observed only for stellate cells (Fig. 2, *E* and *F*). With increasing firing rate, an increasing number of spikes belonged to a cluster (increased plateau level in the ISI distribution in Fig. 2*G*) and average intervals between clusters decreased (shortened tail of the ISI density). Remarkably, the dominant intracluster ISI of stellate cells did not vary with firing rate. In contrast to stellate cells, the ISI density in pyramidal cells was monomodal and the cumulative ISI distribution lacked any plateaus. Higher firing rates simply caused a shift of the distribution to smaller ISI values (Fig. 2*H*).

To summarize, spike patterns in stellate cells exhibited two different timescales reflecting inter- and intracluster ISIs. Increasing the firing rate decreased the intercluster ISIs but did not affect intracluster ISIs. Pyramidal cells, on the other hand, had only a single timescale that directly corresponded to their firing rate.

### Resonate-and-fire model

As the first step of our computational analysis, we reproduced the experimental ISI distributions with a minimal resonate-and-fire model (see methods). Its subthreshold dynamics are given by a linear differential equation that covers both resonant and nonresonant membrane impedances. As shown by previous studies this approach captures the different subthreshold responses of stellate and pyramidal cells surprisingly well, despite its minimal complexity. In addition, the model makes quantitative predictions that have been tested experimentally (Erchova et al. 2004; Schreiber et al. 2004a).

Here, we extend this successful approach to the spiking regime and investigate whether the model also provides a mechanistic understanding of the spike-train statistics of different types of entorhinal cortex cells. To this end, it is important to recapitulate that the parameters of this phenomenological description are allowed to change with the membrane potential so as to account for voltage-dependent cell properties. For each level of membrane depolarization, one set of parameters was therefore estimated (Fig. 3). In both stellate and pyramidal cells, the effective capacity *C* and eigenfrequency /2π depended only weakly on the membrane potential (Fig. 3, *A*–*D*). For the purpose of this study, *C* was thus considered constant. In resonant cells the damping constant γ/*C* decreased with depolarization (Fig. 3*E*), whereas it was approximately constant in nonresonant cells (Fig. 3*F*). This result agrees with the observation that the amplitude of subthreshold resonances increases with depolarization (Erchova et al. 2004). Finally, in both cell classes the intrinsic noise intensity *D* increased strongly with membrane depolarization (note the logarithmic scale in Fig. 3, *G* and *H*), in agreement with findings of Diba et al. (2004) in cultured hippocampal neurons.

To include spike generation in the model without increasing its complexity more than absolutely necessary, the subthreshold dynamics were combined with a constant voltage threshold *x*_{b}. When the membrane potential reaches *x*_{b}, an action potential is elicited, the membrane potential is reset to a fixed value *x*_{0}, and its derivative d*x*(*t*)/d*t* is reset to zero. To account for refractoriness, the subthreshold dynamics begins again after a fixed reset time τ_{r}. These three additional parameters were also estimated from the experimental traces (see Table 1).

### ISI distributions are captured by the resonate-and-fire model

The subthreshold (“resonate”) and spiking (“fire”) parameters of the resonate-and-fire model are obtained at different levels of depolarization. Since the subthreshold parameter values depend on the membrane potential, they cannot be directly inserted into the model. We rather used the ranges of parameters measured across different membrane potentials to constrain the parameter space and searched within this space for the optimal parameter set by simulated annealing (see supplemental material). The criterion of optimality was defined as the lowest KS distance (see methods) between experimental and model ISI distributions. The latter were obtained analytically as the distributions of first passage times using the Stratonovich approximation (as described in the supplemental material). Optimized model ISI distributions are consistent with the experimental data for both cell types (Fig. 4).

Structurally different ISI distributions were obtained with different parameter sets, as shown in Fig. 4*A*. Here, all model parameters were either varied independently (solid line) or constrained by a fixed upper bound for the variance of the subthreshold voltage fluctuations *σ _{x}^{2}* (dashed line). Although both optimizations captured the main ISI peak, the first parameter set reproduced the multimodal probability density better than the second. However, the better fit came at the cost of subthreshold voltage fluctuations much larger than observed in experiments (first model set: σ

_{x}= 5.6 mV; second model set: σ

_{x}= 1.6 mV; experiments σ

_{x}≈ 2.5 mV).

These results can be understood from considering the model features that shape the ISI density. In both experiment and model, spikes are more likely to occur at the maxima of subthreshold oscillations, where the distance to threshold is smallest (see also Alonso and Klink 1993; Klink and Alonso 1993). Within the resonate-and-fire model, the phase of subthreshold oscillations is reset to a fixed value after each spike. Consequently, the initial subthreshold oscillation following a spike is similar for all spikes (despite the presence of noise) and the membrane potential reaches its first maximum roughly at the same time lag *T*, for small to moderate noise levels. *T* corresponds to the intracluster ISI and is equal to the sum of the reset period τ_{r} and the time between the end of the reset period and the first maximum of the following subthreshold oscillation (occurring at approximately half an oscillation period). The second peak in the ISI density is caused by spikes occurring at the second maximum of the subthreshold oscillation. Therefore the distance between the first and the second peaks in the ISI density ℱ(*T*) reflects the period of subthreshold oscillations. However, if no spike is generated within the relaxation time window *t*_{rel} = 2*C*/γ (i.e., not within the first few periods of subthreshold oscillations), their fixed phase relation to the last spike is lost due to the noise contribution. For this reason, the large intervals between spike clusters are random and lead to the observed exponential tail of ℱ(*T*). The number of visible peaks in ℱ(*T*) therefore depends on the ratio between the relaxation *t*_{rel} of subthreshold oscillations and their period. ℱ(*T*) exhibits more peaks for large *t*_{rel} (i.e., when the effective damping coefficient γ/*C* is small).

With this insight, the different fit quality in Fig. 4*A* is easily explained: for small damping the first maximum of a subthreshold oscillation is rather large. To limit the number of action potentials triggered at this maximum, the spike threshold also needs to be high. This, however, implies that large intrinsic noise intensities are required to reach the experimentally measured firing-rate levels. Later on we show how this structural shortcoming of the current model can be overcome by a stochastic spike-generation mechanism.

For nonresonant cells, experimental ISI distributions were reproduced with ease (Fig. 4*B* and Table 1). The explanation is straightforward: for these neurons, the effective damping coefficient was relatively large, leading to a short relaxation time, which canceled the effects of any initial phase-locking; thus after each spike, the membrane voltage quickly returned to its steady-state value. Consequently, the position of the main peak in the ISI distribution was mainly determined by the noise intensity.

Firing-rate–dependent changes in spike patterns and ISI distributions can also be reproduced within the resonate-and-fire framework. In stellate cells (Fig. 2*G*), the position of the first peak of ℱ(*T*) did not vary with firing-rate increases and the tail of the density shortened. Because the subthreshold resonance frequency is fairly independent of the depolarization in stellate cells (Erchova et al. 2004; see also Fig. 3*C*), the model therefore suggests that the reset period τ_{r} also does not strongly change with depolarization. As the noise intensity grows with depolarization (Fig. 3*G*), the overall probability of spike generation increases. Consequently, the exponential tail of ℱ(*T*) becomes shorter. In pyramidal cells, an increased firing rate simply shifts the position of the maximum in the ISI density toward shorter ISIs (Fig. 2*H*): the higher noise level increases the probability of spike generation and thus reduces the length of ISIs.

### Higher-order ISI statistics

For resonant cells, a good agreement between model and experimental ISI distributions has, so far, been possible only at the cost of large voltage fluctuations. Because of the reset of all variables after each model spike, ISIs are independent of each other. The present model thus implements a renewal process. If, however, the assumption of independence of ISIs was incorrect, a modification of the model introducing ISI correlations through a nonrenewal mechanism would be required. This extension to the model may also further improve its overall performance. For this reason, we analyzed the experimental ISIs in more detail.

The probability density of the ISI return map ℱ(*T*_{1}, *T*_{2}) quantifies the joint probability of obtaining two particular successive ISIs (Fig. 5, *A* and *B*). For stellate cells, the strong peak at about (100; 100) ms and the two orthogonal arms up to ISIs of about 500 ms reflect spike clusters. For pyramidal cells, the single broad maximum at about (250; 250) ms mirrors the unimodal ISI distribution (Fig. 5*B*). Both results agree with our previous findings. Next, we analyzed correlation maps that quantify the deviation of a measured ISI sequence from the ISI sequence of a renewal spike train with the same ISI distribution (Fig. 5, *C* and *D*). As revealed by this analysis, experimental ISIs are clearly not independent. To quantify the correlations between more distant spikes, serial correlation coefficients (SCCs) were used (Fig. 5, *E* and *F*). Here, the lag *k* denotes the distance between ISIs. Interestingly, negative SCC values for *k* = 1 (neighboring ISIs) were observed at higher firing rates (3–10 Hz) in both cell types. The correlation time τ_{c} ranged from 200 to 400 ms. For lower firing rates (<2 Hz), however, SCC values for all *k* >0 did not significantly differ from zero. Finally, the number of correlated ISIs grew with increasing firing rate for both cell types (data not shown).

### ISI correlations are captured by an extended resonate-and-fire model

The experimental correlations (Fig. 5, *C*–*F*) cannot be captured by a renewal model. We therefore included a spike-history–dependent stochastic firing threshold: after threshold crossing, spikes are triggered with only a certain probability; the occurrence of several spikes within a time window shorter than the correlation time τ_{c} lowers the probability to generate another spike and, consequently, not every threshold crossing elicits a spike (Fig. 6) . The resulting nonrenewal model allows one to capture the observed ISI correlations as well as to match ISI distributions without excessive voltage fluctuations. As for the renewal model, parameters were optimized within the experimental bounds. Differences between experimental and model ISI distributions and absolute values of SCC (for *k* = 1) were now minimized simultaneously.

Although the experimental voltage fluctuations are no longer exceeded, the ISI distributions of the nonrenewal model agree better with the experimental data than the renewal-based fits (Fig. 7, *A* and *B*). This is because the nonrenewal threshold mechanism limits the number of spikes generated in close succession. It is therefore possible to produce the required number of short intervals and simultaneously ensure a small damping coefficient γ/*C* and low threshold *x*_{b}. In addition, the serial correlation coefficients suggested by the model matched the experimentally obtained results for both cell classes (Fig. 7, *C* and *E*). We also established that the model is consistent with the observation that ISI correlations vanish at lower firing rates (Fig. 7, *D* and *F*). To that end, SCCs in the optimal nonrenewal model with decreased noise (both cell types) and increased damping coefficient (only stellate cells) were analyzed. Although the nonrenewal threshold mechanisms remained unchanged, no ISI correlations were observable at lower firing rates.

It is a property of both the renewal and nonrenewal models that the distance between the first and the second peaks of the ISI density corresponds to the period of subthreshold oscillations. As previously observed in the entorhinal cortex (Erchova et al. 2004) in accordance with the subthreshold part of the model, this oscillation period is also approximately equal to the inverse of the subthreshold resonance frequency. According to the model, the distance between the first and second peaks in the ISI density should therefore be approximately equal to the inverse of the subthreshold resonance frequency. Indeed, our experimental data confirm this prediction (Fig. 8) . Across the electrophysiologically recorded set of resonant stellate cells (natural range of resonance frequency: 5–12 Hz; see Erchova et al. 2004), the inverse distance between the first two peaks in the ISI densities was found to be approximately equal to the measured subthreshold resonance frequency.

It is important to notice that the subthreshold resonance frequency in stellate cells is approximately constant with depolarization (Erchova et al. 2004). This observation implies that estimates of the resonance frequency at less depolarized levels of the membrane potential—and not only those measured close to threshold—correlate with the described firing statistics. In other cells or under different experimental conditions the subthreshold resonance frequency, however, can change with depolarization. In these cases, quantitative aspects of spike clustering may also depend on depolarization. For example, Nolan et al. (2007) reported that in stellate cells of HCN1 knockout mice, the probability of spikes occurring in clusters and the intracluster spike frequency were changed and resonance at the resting membrane potential was abolished. On the other hand, the observed intrinsic membrane-potential oscillations closely mimicked those in wild-type animals near the firing threshold. This indicates that an HCN1-independent slow potassium current is involved in spike clustering in these cells, in full agreement with our general model framework.

## DISCUSSION

The parallel analysis of experimental and modeling data shows that ISI distributions of both stellate and pyramidal cells are strongly influenced by subthreshold processes. Membrane-potential resonance—an important signature of subthreshold dynamics—causes spike clustering. A low-dimensional renewal model constrained by measured experimental parameters reproduces the observed ISI distributions and allows us to quantitatively understand the mechanisms that lead to spike clustering. In particular, we observed that a spike-induced reset of subthreshold oscillations underlies the clustered spike patterns in resonant neurons. In addition, spike-induced processes cause ISI correlations *independently* of the existence of subthreshold resonance. A spike-history–dependent modification of the threshold mechanism in the model captures these higher-order ISI statistics.

### Subthreshold dynamics shape spike-train statistics

This study was performed at membrane depolarizations in the vicinity of the spiking threshold with average firing rates ≤10 Hz. With the majority of in vivo firing rates ranging from 0 to 10 Hz (Sargolini et al. 2006), the range of firing rates explored within our study thus covers the behaviorally relevant regime. Nevertheless, the cells are able to produce fast-spike sequences, in that the maximum instantaneous rates in both cell types were up to 20 Hz (ISIs of ∼50 ms).

A relation between subthreshold resonance and clustered spike patterns has been conjectured previously for neurons in the olfactory bulb (Chen and Shepherd 1997; Desmaisons et al. 1999) and brain stem (Pedroarena et al. 1999; Wu et al. 2001). Our experimental and theoretical findings systematically show that in stellate cells subthreshold resonance is sufficient to explain spike clustering. No spike-induced activation of additional currents is necessary. In particular, the renewal model predicts that without subthreshold resonance no multimodal ISI distributions can be obtained.

Note that the proposed clustering mechanism in stellate cells is different from that of intrinsically bursting neurons. Although the latter also exhibit spike clusters at higher firing rates, the cells reported here fire regularly if sufficiently depolarized (Richardson et al. 2003; Schreiber et al. 2004b). Within the low-firing range, however, stellate cells in entorhinal cortex can potentially profit from the same computational advantages discussed for intrinsically bursting neurons. These include coding with intra- and intercluster time intervals (Kepecs and Lisman 2003), increase of reliability of synaptic transmission (Lisman 1997), and frequency-specific routing of information if ISIs within a cluster of presynaptic neurons match the membrane resonance of the postsynaptic cell (Izhikevich et al. 2003).

### Implication of negative interspike-interval correlations

Negative ISI correlations have been reported for many different cell types, starting as early as Kuffler et al. (1957). Negative ISI correlations are observed even under constant stimulation (Ratnam and Nelson 2000) and it has recently been suggested that they enhance neuronal coding capabilities (Chacron et al. 2001). The correlations between ISIs observed here on a timescale of 200–400 ms are likely to be caused by slowly deactivating currents that are activated by the preceding spike. Such patterns, e.g., slow potassium currents, are known to mediate firing rate adaptation (Brown and Adams 1980; Madison and Nicoll 1984). Interestingly, these currents can also contribute to subthreshold resonance (Hutcheon and Yarom 2000). Our study, however, clearly shows that both phenomena can exist independently of each other: pyramidal cells exhibited negative correlations but did not show any subthreshold resonance.

The mechanistic insight gained from the nonrenewal model explains why the ISI correlations are visible only at higher firing rates in both cell types: spike-induced currents responsible for the temporal relations between spikes are likely to be active across a wide range of depolarizations. Nevertheless, these currents cannot influence spike patterns if the mean ISI is larger than the deactivation timescale. This is true at low firing rates. In contrast, for higher firing rates the mean ISI is sufficiently small for spike-induced currents to interact with spike generation. Along these lines, (*i*) mimicking a decreased firing rate in the optimal nonrenewal model made the negative correlations disappear (Fig. 7) and (*ii*) fits of experimental data with low firing rates were equally good with the renewal and nonrenewal models (data not shown).

Interestingly, for higher firing rates the nonrenewal spike mechanism also improved the representation of ISI distributions within the resonate-and-fire framework. This indicates that correlations between spikes are involved in shaping the fine structure of ISI distributions, whereas the fundamental features, such as spike clustering, are determined by the subthreshold frequency preference.

### Model scope, limitations, and outlook

Why do we use a simple and transparent mathematical model for describing activity in cells of the entorhinal cortex? First, on a practical level, the model is analytically tractable and, at the same time, permits an accurate phenomenological description of both resonant and nonresonant neurons. Due to its generality it can be applied to different cell types, independently of a particular ionic realization of subthreshold resonance properties. Noticeably, the model parameters used to describe a given neuron are tightly constrained by experimental data and obtained from the very same cell. This is not possible for biophysically realistic models whose many parameters cannot be estimated with measurements from a single neuron.

Second, on a conceptual level, the current model provides a quantitative framework to investigate the relation between subthreshold and spiking activity. Here, neglecting microscopic details allows one to reveal general neural mechanisms that go far beyond any specific biophysical processes. This is important because the combinations of conductances may be vastly different between cells, but still lead to the same input–output transformation (Marder and Prinz 2002). In particular, the model clearly demonstrates that the well-known phenomenon of subthreshold resonance can be sufficient to fully explain clustered spike patterns and multimodal ISI distributions. Furthermore, it enables us to disentangle the contributions of subthreshold and spike-induced dynamics to the observed spike patterns. Although it is possible that specific ionic currents contribute to both spike-induced and resonance-related dynamics, our model shows that separating these dynamics on a phenomenological level still captures firing statistics. In particular, describing spike-induced effects by a reset time and a spike-history–dependent spike probability results in a surprisingly good approximation of experimental data.

Interestingly, the mathematical model also suggests that the clustered spike patterns of EC stellate cells are obtained because each spike induces a phase reset of the intrinsic subthreshold oscillation so that the subsequent spike is likely to be generated after a fixed time lag, roughly equal to the reset period plus half an oscillation period. This picture differs from an alternative situation where spikes, such as those in the CA1 region, are preferentially generated without phase reset at the maxima of an ongoing subthreshold oscillation superimposed with a slowly varying excitatory input (see, e.g., Buzsáki 2002; Leibold et al. 2008). If that description also applied for spikes caused by the combination of intrinsic subthreshold oscillations and constant current injections, as investigated in the present study, the main peak of the intracluster ISI distribution should coincide with the subthreshold oscillation period. Our data show, however, that the dominant intracluster ISI is shorter than the oscillation period (Fig. 2). In contrast, the time between the first and second peaks is similar to the oscillation period, in agreement with the model.

We also point out that, according to the proposed model, the existence of persistent subthreshold oscillations relies on the presence of damped oscillations of the corresponding noiseless system (*D* = 0). For the observed oscillations to be visible, however, deviations from the rest state are required, which for synaptically blocked cells are triggered when (*i*) the cell is depolarized and (*ii*) intrinsic noise allows large deviations from the steady state that effectively translate into persistent oscillations. The difference to models like that of Fransén et al. (2004) and Rotstein et al. (2006) lies in the fact that persistent subthreshold oscillations in such a model can exist independently of noise. In agreement with our view that noise is needed to trigger the damped oscillations, Dorval and White (2005) have shown that channel noise from persistent sodium channels is essential for subthreshold membrane potential oscillations.

Within the current model framework, spike-induced changes in the voltage trajectory are described by a constant reset period and a fixed reset voltage. Spike-induced currents, such as those underlying afterhyperpolarization, may, however, be active longer than the model reset period and thus also contribute to the generation of later spikes. Other nonrenewal models either explicitly include slow spike-induced ionic currents (Benda and Herz 2003; Izhikevich 2004) or modify the threshold values dynamically (Chacron et al. 2004). These models, however, are not consistent with two observations in our experiments: First, the firing threshold was not higher for shorter than that for longer ISIs, as would be expected with a dynamic threshold. Second, the voltage sometimes crossed the threshold without producing spikes, resulting in “silent overshoots,” as shown in Fig. 6. The stochastic threshold used in this study accounts for both observations and directly reflects channel noise at the spike-initiation zone. This noise is particularly high near the firing threshold and known to severely influence spike generation (Schneidman et al. 1998).

As demonstrated by our results, spike patterns of single EC stellate and pyramidal cells are shaped by a combination of subthreshold resonances and spike-induced ionic currents. The generality of the model suggests that similar mechanisms may also be relevant for other neural systems. It remains an open challenge—to be tackled in future studies—whether the signatures of resonance and nonresonance, respectively, are behaviorally relevant. Were this indeed the case, and recent results on grid cells point in that direction (Giocomo et al. 2007), information processing in the entorhinal cortex might be understandable on the basis of a simple mathematical framework. If, on the other hand, behaviorally relevant time-varying inputs generate responses that cannot be understood in terms of a resonate-and-fire mechanism, the interesting question would remain why stellate and pyramidal cells exhibit such strikingly different dynamics, both below and above firing threshold.

Findings from the present study support the first of these two alternatives: all our data suggest that both stellate and pyramidal cells respond only rather weakly to high-frequency inputs. In particular, we did not observe any resonance or strong intrinsic oscillation of the membrane potential above 20 Hz. Rapid synaptic input fluctuations, which might occur under in vivo conditions, would thus only marginally influence the cells's firing activity. On the other hand, slower changes of the mean membrane depolarization are likely to be encoded in the cell-type–specific output patterns of stellate and pyramidal neurons. In such a scenario, the current model framework offers a powerful tool for studying neural dynamics and signal processing under artificial as well as behaviorally relevant conditions. The simplicity of the model will also make it possible to readily extend those studies from the single neurons to the network level.

## GRANTS

This work was supported by Deutsche Forschungsgemeinschaft Grants SFB 555 and SFB 618, Bernstein Center for Computational Neuroscience Berlin, Biotechnology and Biological Sciences Research Council, Engineering and Physical Sciences Research Council, Medical Research Council/Doctoral Training Center at the University of Edinburgh, Lloyd Trustee Savings Bank Group plc, and a personal fellowship to I. Erchova from the Royal Society of Edinburgh.

## Acknowledgments

We thank U. Heinemann and I. M. Sokolov for many stimulating discussions and M. Stemmler for comments on a previous version of the manuscript.

Present address of A. Herz: Department of Biology, Ludwig-Maximilians-Universität Munich and Bernstein Center for Computational Neuroscience, Munich, Germany.

## Footnotes

↵* These authors contributed equally to this work.

↵1 The online version of this article contains supplemental data.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society