Computational Reconstruction of Pacemaking and Intrinsic Electroresponsiveness in Cerebellar Golgi Cells

The Golgi cells have been recently shown to beat regularly in vitro (Forti et al., 2006. J. Physiol. 574, 711–729). Four main currents were shown to be involved, namely a persistent sodium current (I Na-p), an h current (I h), an SK-type calcium-dependent potassium current (I K-AHP), and a slow M-like potassium current (I K-slow). These ionic currents could take part, together with others, also to different aspects of neuronal excitability like responses to depolarizing and hyperpolarizing current injection. However, the ionic mechanisms and their interactions remained largely hypothetical. In this work, we have investigated the mechanisms of Golgi cell excitability by developing a computational model. The model predicts that pacemaking is sustained by subthreshold oscillations tightly coupled to spikes. I Na-p and I K-slow emerged as the critical determinants of oscillations. I h also played a role by setting the oscillatory mechanism into the appropriate membrane potential range. I K-AHP, though taking part to the oscillation, appeared primarily involved in regulating the ISI following spikes. The combination with other currents, in particular a resurgent sodium current (I Na-r) and an A-current (I K-A), allowed a precise regulation of response frequency and delay. These results provide a coherent reconstruction of the ionic mechanisms determining Golgi cell intrinsic electroresponsiveness and suggests important implications for cerebellar signal processing, which will be fully developed in a companion paper (Solinas et al., 2008. Front. Neurosci. 2:4).


INTRODUCTION
Oscillations are one of the most prominent features of brain activity (Buzsaki, 2006). Numerous neurons in the brain show pacemaker activity exploiting a variety of ionic mechanisms (Hutcheon and Yarom, 2000;Koch, 1999;Lliná s, 1988;Yamada et al., 1998). In a well-known case represented by relay thalamic neurons (RTN), large low-voltage oscillations are generated by the interplay of h-current with lowthreshold activated (LVA) calcium channels McCormick and Huguenard, 1992;McCormick and Pape, 1990). The h-current determines the depolarizing ramp leading to LVA channel activation and to the calcium spike, which triggers a spike bursts before inactivating and concluding the cycle. These pacemakers cannot be easily perturbed and drive fast sodium spike bursts. Some other, probably more diffused pacemaker neurons, show smaller oscillations in the subthres-hold membrane potential region generating a regular low-frequency beating, usually in the theta-frequency band. These include several neuronal types in the basal ganglia and brain cortex (Bennett et al., 2000;Dickson et al., 2000a;Surmeier et al., 2005;Wilson and Callaway, 2000). The mechanisms are thought to employ primarily persistent sodium currents in combination with some slow repolarizing mechanism, although oscillations driven by high-voltage activated calcium channels have also been reported. The variety of pacemaker mechanism may be related to the additional behaviors the neuron is able to express. It is therefore important to determine how a specific neuron generates pacemaking and how this combines with other aspects of electroresponsiveness.
The Golgi cell (Golgi, 1883) is an inhibitory interneuron regulating information flow in the granular layer circuit of cerebellum (Eccles et al., 1967;Marr, 1969), which has been recently shown to beat regularly in vitro , see also Dieudonné , 1998). Electrophysiological and pharmacological analysis  showed that four main currents were involved, namely a persistent sodium current (I Na-p ), an h current (I h ), an SK-type calcium-dependent potassium current (I K-AHP ), and a slow M-like potassium current (I K-slow ). However, the precise role of these currents was unclear. In particular, both I Na-p and I h may be able to drive the oscillation and the repolarizing action of I K-AHP and I K-slow may also be redundant.
As a first step for investigating the mechanisms of Golgi cell excitability, we have developed a computational model incorporating state-of-the-art knowledge on Golgi cell ionic channels and electro-responsiveness. The model predicted that pacemaking was sustained by subthreshold oscillations generated by I Na-p and I K-slow . I h was less important but maintained the model into the appropriate membrane potential range for oscillation. The oscillations were tightly coupled to spikes thereby causing a prominent activation of I K-AHP and a stabilization of the pacemaker cycle. The combination with other currents, in particular a resurgent sodium current (I Na-r ) and an A-current (I K-A ), allowed a precise regulation of response frequency and delay. These results provide a coherent reconstruction of the ionic mechanisms determining Golgi cell intrinsic electroresponsiveness. The implications of this complex ionic complement will be further clarified in a companion paper (Solinas et al., 2007), in which I K-AHP and I K-slow will assume the additional and more specific roles of determining phase-reset and theta-frequency resonance.

MATERIALS AND METHODS
This work presents a combined experimental and modeling analysis of Golgi cell electroresponsiveness. In order to develop a comprehensive hypothesis on how the Golgi cell generates its excitable response, an extensive re-analysis of patch-clamp recordings and pharmacological data obtained in a previous paper  has been complemented with a computational model. The model allowed us to investigate neuronal behaviors hard to resolve experimentally, like the apparent close coupling of spikes with intrinsic oscillatory mechanisms and the effect of persistent Na þ currents. Moreover, while conductance measurements and dynamic-clamp have been successfully applied to electrotonically compact neurons (e.g., in isolated Purkinje cells: Raman and Bean, 1997), the large clamp escape in Golgi cells in slices (data not shown) prevents the application of these methods.

Golgi cell modeling
A computational model of rat Golgi cell electroresponsiveness was constructed using the NEURON simulator (version 5.9; Hines and Carnevale, 2001). The model consisted of five compartments allowing for a minimal description of the Golgi cell electrotonic structure. The soma was connected to an ''electrode'' to reproduce realistic current-clamp conditions. All voltage-and Ca 2þ -dependent mechanisms (see below) were placed in the somatic compartment. With this approach, the model reproduced satisfactorily basic aspects of Golgi cell electroresponsiveness elicited by somatic current injection. It should also be noted that, although differences have been reported at the histochemical level (Geurts et al., 2001;Simat et al., 2007), the basic electroresponsive properties were homogeneous in a large majority of Golgi cells , see also figures in this paper). Thus, we have reconstructed a ''canonical'' Golgi cell model which simulates the electrophysiological behavior for this major Golgi cell subgroup and lays within the scatter of physiological parameter values. The model was calibrated on the cellular response corrected for the liquid-junction potential.
The model included 12 voltage-dependent conductances placed into the soma (Table 1). Gating kinetics were corrected using a Q 10 ¼ 5 for I Ca-LVA activation (Destexhe et al., 1994) and a Q 10 ¼ 3 for all the other currents according to the relation Q ðT sim ÀTexpÞ=10 10 (Gutfreund et al., 1995) to account for temperature differences between experimental recordings and the model (238C). Nernst equilibrium potentials were calculated from ionic concentrations used in current-clamp recordings. Only the Ca 2þ Nernst potential was updated during simulation. The maximum density voltage-dependent conductances were regulated to fit the Golgi cell responses to various stimulations (Achard and De Schutter, 2006;Traub and Lliná s, 1979;Traub et al., 1991;Vanier and Bower, 1999).
Ionic currents were modeled following the Hodgkin and Huxley formulation of ion channel gating using mathematical methods reported previously (D'Angelo et al., 2001;Nieus et al., 2006). Voltage was obtained as the time integral of the equation (Yamada et al., 1998) where V is membrane potential, C m membrane capacitance, g i are ionic conductances, and V i reversal potentials (the subscript i indicates different channels), and i inj is the injected current. I AHP was simulated using a Markov gating scheme reproducing the dynamics reported by Hirschberg et al., (1998), while all other membrane conductances were represented using Hodgkin-Huxley-like models (Hodgkin and Huxley, 1952) of the type where G max i is the maximum ionic conductance, x i and y i are state variables (probabilities ranging from 0 to 1) for a gating particle, and z i is the number of such gating particles in ionic channel i. x and y (the suffix i omitted) were related to the first order rate constants a and b by the equations where a and b are functions of voltage. The equations used to parameterize a and b and the state variables x 1 , t x , y 1 , t y for different ionic channels are shown in Àb Ca ð½Ca 2þ À½Ca 2þ 0 Þ where d is the depth of a shell adjacent to the cell surface of area A, b Ca determines the loss of Ca 2þ ions from the shell approximating the effect of fluxes, ionic pumps, diffusion, and buffers (De Schutter and Smolen, 1998;McCormick and Huguenard, 1992;Traub and Lliná s, 1979), and [Ca 2þ ] 0 is the resting [Ca 2þ ].

Model mechanisms
Passive properties and model structure. In cerebellar Golgi cells, charge displacement during voltage steps followed a tri-phasic trajectory (see Figure 1). According to Dieudonné (1998), the fastest transient was related to the soma, the intermediate to the dendrites, and the slowest to the axon. According to proportions reported by Dieudonné (1998) the capacitances of soma, dendrites, and axon were set to 23, 32, and 90 pF, respectively, summing up to a total cell capacitance of 145 pF . The specific axial resistance of axon and dendrites was set to 100 V Á cm (Roth and Hä usser, 2001), the specific membrane resistance was set to 47.6 KV Á cm 2 (corresponding to the experimental determination of leakage reported below) and specific membrane capacitance was set at 1 mF/cm 2 . Compartments equivalent to the complex structure of dendrites and axons (Koch, 1999;Yamada et al., 1998) were generated to comprise a spherical soma (diameter ¼ 27 mm), a single axonal compartment (1200 mm long and 2.4 mm thick), and three dendritic compartments (each 113 mm long and 3 mm thick). With these parameters, input capacitance and input resistance in the model equaled those of Golgi cells and the model could appropriately reproduce current transients recorded in vitro in the subthreshold region ( Figure 1). The micropipette was simulated as a passive cable with sealed end, null capacitance, infinite wall resistance and an axial resistance of 25 MV equivalent to ''access resistance'' (Neher and Sakmann, 1995).
Voltage-dependent mechanisms. These were based on experimental observations including current-clamp  and preliminary voltage-clamp recordings. A first current set was needed to generate action potentials. This consisted in a standard description of I Na-t and I K-V and was completed through the introduction of I Ca-HVA and I K-C , which improved the fast phase (1-2 ms following the upstroke) of spike afterhyperpolarization (AHP) preventing Na þ channel inactivation and discharge blockage at high frequency. I K-A regulated the first spike delay and I Na-r intensified high-frequency discharge. I Ca-LVA was needed to enhance rebound depolarization. A second current set was needed to generate pacemaking and, according to previous experimental analysis , included I Na-p , I h , I K-slow , and I K-AHP . The procedures  The table reports the equations used to calculate a x , b x , a y , and b y (see Equations (3) and (4)) for the membrane conductances used in the model (in a few cases x 1 and y 1 are reported). Only, the SK channel was built using Markov multistate transition model in order to use the data available in literature (Santoro et al., 2000). The number of gating particles (n), maximum conductance (G max ), and ionic reversal potential (* resting value for Ca 2þ ) used to calculate ionic currents are also shown. The opening and closing rates are corrected for temperature (23 8C).
The presence of I h in Golgi cells was reported by Dieudonné (1998) and confirmed by Forti et al., (2006). I h in Golgi cells is manifest during the application of hyperpolarizing current steps causing sagging inward rectification with both a rapid and a slow phase. Slow and fast response components were reproduced by two voltage-dependent mechanisms corresponding to HCN1 and HCN2 currents, respectively (Santoro et al., 2000, see also Chan et al., 2004;Dickson et al., 2000b;Luthi and McCormick, 1998;McCormick and Pape, 1990;McCormick and Huguenard, 1992). After regulating current densities to match the sag and the steady-state level of hyperpolarizing responses, the model generated a robust rebound depolarization, as observed experimentally (a reinforcement of the early rebound phase was obtained with I Ca-LVA , see below).
I Na-p . I Na-p was matched to the I-V curve reported by Forti et al., (2006). I Na-p kinetics were derived from the model previously developed for cerebellar granule cells (D'Angelo et al., 2001;Magistretti et al., 2006;Nieus et al., 2006) and pyramidal cells (Gutfreund et al., 1995).
The presence of a slow repolarizing retigabine-sensitive M-like current (I K-slow ) in Golgi cells was suggested by pharmacological observations . The I K-slow model was obtained from previous models (D 'Angelo et al., 2001;Gutfreund et al., 1995) and maximum conductance was set in proportion to I Na-p (Nieus et al., 2006). I K-slow contributed to shape the ISI and spike frequency adaptation.
I K-AHP . The presence of a slow Ca 2þ -dependent repolarizing apaminsensitive K þ current (encoded by SK-type channels; for a review see Bond et al., 2005;Stocker, 2004) sustaining the slower phase (peaking 30-40 ms after spike peak) of the spike AHP in Golgi cells was indicated by pharmacological observations . Although, in rat Golgi cells, mRNAs encoding the SK3 subtype are particularly abundant (Stocker and Pedarzani, 2000), the corresponding current in the model, I K-AHP , was simulated using the available SK2 model reported by Hirschberg et al., (1998). I K-AHP was activated by Ca 2þ entering through HVA Ca 2þ channels (Marrion and Tavalin, 1998;Stocker, 2004). After regulating intracellular Ca 2þ dynamics, I K-AHP was tuned in order to fit the slow AHP phase and spike frequency adaptation.
I Ca-LVA . Preliminary observations suggest that Golgi cells express LVA calcium channels reinforcing rebound excitation on return from hyperpolarization (Forti, Cesana and D'Angelo, unpublished observations). Indeed, the model could not properly reproduce the initial reinforcement of rebound firing with I h alone. The I Ca-LVA model was taken from Destexhe et al., (1994) and current density adjusted to appropriately reinforce rebound excitation. Due to its voltage-dependent inactivation, I Ca-LVA did not contribute to any other excitable process of the Golgi cell.
Spike-related currents. The complement of spike-related currents was reconstructed with mechanisms taken from a previous granule cell model (D'Angelo et al., 2001;Nieus et al., 2006). The density of the transient Na þ current (I Na-t ) was adjusted to reproduce the action potential threshold and overshoot. A preliminary investigation revealed that, as well as granule cells, Golgi cells also express a proportionate resurgent Na þ current (I Na-r ) component (Magistretti, Castelli and D'Angelo, unpublished observations). Its insertion into the model was required to enhance the initial doublet in response to depolarization raising the maximum firing rate into the physiological range measured in vitro (up to $200 Hz, Forti et al., 2006). The delayed rectifier K þ current, I K-V , was needed to generate spike repolarization and balanced with I K-C to match the spike shape. The requirement of a Ca 2þ -dependent voltage-dependent K þ current, I K-C , emerged from the inability of I K-V to fully account for the fast phase of spike AHP and to guarantee Na þ channel re-priming at high discharge frequency. I K-A was needed to regulate spike delay in response to depolarizing current injection and during rebound excitation.
Ca 2þ dynamics. Ca 2þ dynamics were adapted to yield Ca 2þ transients of about 1 mM following spike-dependent opening of HVA channels (Gabbiani et al., 1994). Parameters used in Equation (6) were similar to those used in other neurons (e.g., De Schutter and Smolen, 1998;D'Angelo et al., 2001;Traub and Lliná s, 1979;Traub, 1982). With this asset, there was calcium accumulation during high-frequency discharge but not when the model fired at its intrinsic pacemaker rhythm. The computed intracellular Ca 2þ concentration [Ca 2þ ] i activated I K-C just after the spike upstroke and generated the fast AHP. I K-AHP was activated using 1/3 of the computed [Ca 2þ ] i value on the assumption that local concentration is smaller on SK than on channels supporting I K-C (see Stocker (2004) for a discussion on effective local [Ca 2þ ] i concentrations). This seems reasonable in view of the higher [Ca 2þ ] i affinity of SK channel subtypes with respect to channels supporting the fast spike AHP (reviewed in Stocker, 2004). The slow dynamics of I K-AHP , which contributed to the slower phase of the spike AHP, was accounted for by slow channel closing kinetics rather than by the permanence of elevated Ca 2þ levels. I K-AHP activation was enhanced during high-frequency repetitive discharge due to the associated Ca 2þ build-up, thereby determining spike-frequency adaptation (see Figure 4). It should be noted that, in the model, I Ca-LVA was not allowed to contribute to the Ca 2þ pool controlling I K-C and I K-AHP , since such a coupling would unrealistically depress rebound excitation (data not shown). The recording (black trace) shows a current transient elicited in a Golgi cell by a 50-ms voltage step from À70 to À80 mV. The model (gray trace) faithfully matches the response. The plot shows the relationship between input resistance (R in ) and input capacitance (C in ) in 39 Golgi cells (black circles, the larger one corresponding to the trace shown on top). The model prediction (large gray dot) falls in the middle of the experimental data distribution (p > 0.8, t-test).

Solinas et al.
Membrane noise in vitro. In order to reproduce the ISI irregularity observed during Golgi cell recordings in vitro (CV ISI ¼ 0.13, median value from loose cell-attached recordings, Forti et al., 2006; CV ISI ¼ 0.16, median value from WCR, this study), Gaussian noise was added to the model in selected simulations aimed at reproducing the effect of apamin application (cf. Figures 2E and 5C). The NEURON random number generator was used to produce a sequence of numbers cast from a Gaussian distribution with 0 mean and unitary variance generating a sequence lacking temporal correlation. The random numbers were scaled by a 32 pA coefficient yielding a noise current injected into the soma. In the canonical model, this produced a CV ISI ¼ 0.16. The Golgi cell model NEURON code is available from ModelDB (senselab.med.yale.edu/ModelDb/).

Data analysis
The processing of experimental and simulated traces was automated using dedicated scripts written in MATLAB (MathWorks, Natick, MA, USA). As in Forti et al., (2006), spike threshold was detected at the point where the depolarization rate reached 5 mV/ms. Spike maximum and minimum (AHP trough) were computed and the AHP rise-time was estimated as the time needed to reach the AHP trough from threshold crossing during the down phase of the spike. The sag amplitude was calculated as the difference between the minimum and steady-state value of membrane potential during responses to hyperpolarizing currents. Firing rate adaptation in response to depolarizing currents was calculated as the ratio of the steady-state and initial discharge frequency (i.e., the inverse of the last and first ISI).
The match of model to experimental data was evaluated by testing the null hypothesis that their difference was a random sample from a normal distribution (Lilliefors test) with mean equal to 0 and unknown variance (Student's t-test). All data are reported as mean AE SD.

RESULTS
The Golgi cell was recently proposed to express a characteristic set of ionic currents (comprising I Na-p , I h , I K-AHP , and I K-slow ; Forti et al., 2006). Each one of these is endowed with unique gating and kinetics properties (see for a comprehensive description Koch, 1999;Yamada et al., 1998): I Na-p activates rapidly with depolarization in the subthreshold region, I h slowly deactivates with depolarization, I K-AHP activates briskly upon Ca 2þ entry and remains activated for tens of milliseconds, and I K-slow activates slowly and progressively in the subthreshold region during sustained depolarization.
In order to elaborate a hypothesis on how these ionic mechanisms interact, the electrophysiological and pharmacological analysis of intrinsic excitability was assisted by a computational model incorporating prototypical ionic current properties (see Methods, Figure 1 and Table 1).

Modeling the basic properties of Golgi cell electroresponsiveness
The basic properties of Golgi cell electroresponsiveness in acute brain slices (henceforth called in vitro) were inferred from the previous paper of Forti et al., (2006) and from novel WCR data, and could be summarized as follows. (i) Golgi cells paced rhythmically at a frequency between 1 and 8 Hz at room temperature ( Figure 2A). In the pacemaking regime, spikes were followed by a protracted AHP ( Figure 2B). The average ISI distribution obtained from different cells in WCR showed a prominent peak around 200 ms ( Figure 2C). There was no clear correlation between the Relationships among ISI parameters. Note that the model data points (large gray dots) fall within the distribution of spike amplitude versus spike threshold and within the distribution of AHP trough versus AHP rise-time measured experimentally (black circles, the larger ones corresponding to the trace in A-B). The model did not significantly differ from the data (p > 0.12, t-test). (E) The experimental relationship between CV ISI and firing rate (n ¼ 47 WCR) showed a negative correlation. The model was injected with a noisy current (32 pA SD, see Methods) and the spontaneous frequency was varied among the values observed during WCRs by changing the leakage reversal potential from À67 to À55 mV. Simulations (gray line) show that the model could appropriately fit the experimental measurements.
rise-time and depth of AHP or between spike size and threshold ( Figure 2D). The variability of the ISI in each cell was inversely correlated with the average firing frequency ( Figure 2E), and was small (CV ISI ¼ 0.16, median value from WCR). (ii) When injected with hyperpolarizing currents, Golgi cells showed sagging inward rectification and, at the end of hyperpolarization, generated a rebound excitation transiently intensifying spike generation ( Figure 3A). Rebound excitation consisted of an acceleration of spike frequency, which progressively decayed back to baseline levels. There was no clear correlation between the amplitude of the hyperpolarizing sag and the first spike delay or the duration of the first ISI ( Figure 3B). (iii) When injected with a depolarizing current, Golgi cells generated repetitive spike discharge showing fast adaptation at highcurrent intensity. At the end of discharge, Golgi cells generated an AHP ( Figure 4A). The first spike delay decreased progressively with current injection ( Figure 4B, left) while the firing rate increased ( Figure 4B, middle). The comparison of instantaneous frequencies at the beginning and end of discharges showed that adaptation could reduce firing frequency by 78 AE 4%. Inspection of Figure 4A also showed that the discharge initiated with a high-frequency doublet and that most adaptation occurred after the first few spikes resulting in an adaptation that was inversely dependent on the initial firing rate ( Figure 4B, right).
Since the inspection of data plots in Figures 1-4 did not allow to discern subclasses of Golgi cells on the basis of pacemaking frequency, ISI shape, rebound excitation or discharge frequency, a ''canonical'' model was generated to simulate the average electrophysiological behavior. The model was able to quantitatively reproduce all the properties summarized above, as demonstrated in Figures 1-4. The average pacemaker frequency coincided with the mode of the experimental data distribution, and model predictions about spike shape, inward rectification, rebound excitation, firing delay, and frequency adaptation during depolarization fell within the experimental data distributions with p > 0.12 (t-test). It should be noted that the AHP at the end of a depolarizationinduced discharge deviated from the shape observed in neurons ( Figure 4A) suggesting that the model does not completely represent the Ca 2þ dynamics developing over longer time-scales. This is, however, a minor problem in the present context, since the physiological implications of AHP are considered in response to short depolarization and the model could indeed precisely predict the subsequent silent pause (see Solinas et al., 2007).

Testing the model versus pharmacological challenges
The adequacy of the model in terms of intrinsic mechanisms was further tested by comparison with pharmacological observations using specific ionic channel blockers (cf. Forti et al., 2006). A 50% I h reduction, simulating sub-maximal block by ZD7288, reduced the pacemaker frequency (73% of control) and slowed down the ISI trajectory ( Figure 5A), similar to in vitro observations. As expected from their gating properties reported in vitro, 50% I Na-p reduction arrested pacemaking ( Figure 5A). The action of retigabine, an agonist of M-type currents which reduced firing frequency and altered the ISI during current injection in vitro, was matched by a À6.5 mV shift in I K-slow activation (Tatulian et al., 2001) ( Figure 5B). Finally, since the most evident effect of the SK channel blocker apamin was to reduce pacemaker regularity, the role of I K-AHP in the model was investigated after adding noise, to reproduce the standard deviation of membrane potential during the ISI ( Figure 5C; see Methods for details). The simulated pacemaking thus acquired a slight irregularity (CV ISI ¼ 0.16), which was markedly increased by I K-AHP blockage (CV ISI ¼ 0.36), as observed in the experiments in vitro. It is thus probable that the AHP, by temporarily displacing membrane potential below threshold, prevents the perturbing action of stochastic channel gating and spontaneous synaptic currents in the subthreshold regime. In summary, simulations show that the model predictions for the effects of I Na-p , I k-slow , I k-AHP , and I h on pacemaking agree with pharmacological analysis in vitro.

Model currents and intracellular calcium
The model currents and intracellular calcium changes during Golgi cell activity are shown in Figure 6. The currents dominating the ISI during pacemaking were small ( 5%) compared to those generating the spikes. During the ISI, I Na-p showed a progressive increase, while I h remained almost constant. I k-AHP was driven by the spike, arose rapidly and then decreased toward the middle of the ISI. I k-slow was driven by depolarization in the second part of the ISI and showed a progressive increase. It should be noted that I Ca-HVA activated during the spike upstroke. During pacemaking, intracellular Ca 2þ returned to basal levels within an ISI, but during high-frequency bursts Ca 2þ could accumulate and thereby intensify and protract I K-AHP activation. Spike-frequency adaptation during discharge elicited by current injection was caused by I K-AHP and I K-slow . Rebound excitation was determined by I h deactivation and intensified by I Ca-LVA . Certain currents accurately regulated the responses to current injection (Figure 7). During depolarization, I K-A delayed the first spike and I Na-r enhanced the initial doublet. As well as during depolarization, I K-A contributed to regulate the delay of rebound excitation (data not shown). The WCR from a Golgi cell shows sagging inward rectification in response to hyperpolarizing current injection and, at the end of the hyperpolarization, rebound excitation with an early and a protracted phase of intensified firing. Simulations of this specific experiment (parameters were adjusted to mimic the experimental trace in detail: leakage reversal potential À64.5 mV, I K-slow 90%, I h-HCN1 60 %, I h-HCN2 12.5 % of their value in the canonical model) show that the model could faithfully reproduce sagging inward rectification and rebound excitation (gray trace). The step-current protocol is shown at the bottom. (B) The plots report the time to first spike and the first ISI during rebound excitation as a function of sag amplitude in 10 Golgi cells (black circles, the larger one corresponding to the trace in A) stimulated using the protocol shown in A. Note that the canonical model (large gray dots) did not significantly differ from the data (p > 0.13, t-test). Even the model adjusted to match the specific trace shown in A (large gray circle, see details reported in A) falls within the experimental data scatter.

Robustness and the canonical neuron
The simulations we have presented till now lead to two conclusions. (i) Although the currents introduced into the model (on the basis of experimental observations) controlled characteristic aspects of electroresponsiveness (see also Solinas et al., 2007), they could have multiple and redundant effects. (ii) The canonical Golgi cell model is one that succeeds in approximating the behavior of the experimental population of Golgi cells. The robustness of this model should be tested by checking its ability to reproduce the variability of behaviors measured in different cells when varying model parameters. We therefore measured robustness of the model by assessing its ability to maintain typical Golgi cell characteristics when key parameters were varied. This aspect is shown in Figure 8, in which certain current densities or the time constant of [Ca 2þ ] i are modified in turn. The analysis in Figure 8 confirms the relative importance of different currents in generating Golgi cell behavior reported in Figure 5. The model was quite sensitive to changes in I Na-p conductance, followed by changes in I K-slow . Calcium dynamics, determined by the conductance of I Ca-HVA and the decay constant b Ca were also relatively important, while the model was quite insensitive to changes in I h conductance. An example was presented in Figure 3, in which some parameters were modified to adapt the canonical model to mimic an individual Golgi cell response in detail. Notice that a more extended robustness assessment (Achard and De Schutter, 2006;Druckmann et al., 2007) would imply also testing the effect of conjoint variations of free parameters. This is a complex task that may be addressed in future research.
Coupling of subthreshold membrane potential oscillations with spikes Given the good response of the model over a wide range of experimental behaviors, it could be used to predict the mechanism of pacemaking. In several neurons (e.g., in amygdaloid neurons (Papae and Driesang, 1998), in enthorinal cortex neurons (Alonso and Lliná s, 1989;Dickson et al., 2000a), or in inferior olivary neurons (Lliná s and Yarom, 1986)) pacemaking is clearly related to subthreshold membrane potential oscillations. However, these could not be isolated in WCRs by gradually hyperpolarizing Golgi cells with steady currents. This situation clearly resembles that observed in subthalamic neurons, in which pacemaking was extremely regular and subthreshold oscillations could not be decoupled from spikes (Bevan and Wilson, 1999). In the model, after switching-off spike-related mechanisms (i.e., I Na-t , I K-V , and I K-C ), subthreshold oscillations (30 mV peak-to-peak centered around threshold) emerged at nearly the same frequency (5 vs. 4.46 Hz) of spontaneous firing ( Figure 9A). Several putative mechanisms might be responsible for these subthreshold oscillations in Golgi cells, exploiting the specific properties of the two depolarizing currents, I Na-p and I h , and of the two Simulations show that the model could faithfully reproduce this behavior (gray traces). Both the real cell and the model were held at À70 mV by injecting a hyperpolarizing current to prevent spontaneous firing (the step-current protocol is shown at the bottom). (B) In the plots, the response of the model (solid gray line) to current injection is compared to the data collected in experimental recordings (dots; n ¼ 10 Golgi cells). Simulations (gray lines) show that the model could appropriately fit the experimental measurements of : -first spike delay versus injected current (left: p > 0.1, t-test) -instantaneous and steady-state frequency (the inverse of the first and last ISI, experimental data correspond to black dots and diamonds, modeling data correspond to solid and dashed gray lines, respectively) versus injected current (middle plot: p > 0.1, t-test) -the adaptation factor (last ISI/first ISI) versus initial firing rate (p > 0.9, t-test).
slow repolarizing currents, I K-slow and I K-AHP . The role of each of these currents was considered in turn.
Blocking either I Na-p or I h modified the subthreshold oscillations although in quite a different way. By blocking I h , the oscillation slowed down (from 4.46 to 3.3 Hz) and the original frequency could be restored by injecting a tonic depolarizing current ( Figure 9A). Conversely, following 50% I Na-p reduction, the oscillation disappeared. In this case, a large tonic current injection compensated the partial block only to a very small degree ( Figure 9B). Concerning the repolarizing feedback, blocking either I K-slow or I K-AHP increased oscillation amplitude but with opposite effects on frequency (I K-slow block reducing it to 2.9 Hz and I K-AHP block increasing it to 4.5 Hz), while blocking both I K-slow and I K-AHP led to a cessation of oscillation ( Figure 9B).
The results reported in Figure 9A, B suggest that I Na-p and I h played distinct roles in Golgi cell subthreshold oscillations. The current plots in Figure 9C and the phase diagram in Figure 8D show that I Na-p underwent pronounced modulation during the cycle, while I h did not. Thus, while the brisk I Na-p activation was able to drive the cycle, the slow activation of I h provided a nearly tonic depolarization setting the cell into the appropriate membrane potential range. Concerning feedback repolarization, both I K-slow and I K-AHP , which also showed pronounced dynamics during the cycle ( Figure 9C, D), contributed to the feedback repolarization. The differential action of the two K þ currents on cycle shape and frequency reflected their different voltage dependence and kinetics. It should be noted that the subthreshold oscillation reached potentials normally visited by spikes (about À30 mV, cf. Figure 2A), so that I K-AHP could be activated causing a rapid and intense repolarization.
In accordance with the pharmacological tests reported in Figure 5 and the mechanisms reported in Figure 6, simulations shown in Figure 9 suggested that pacemaking was driven by subthreshold oscillations tightly coupled to spikes.

DISCUSSION
In this paper, we show that the theta-frequency pacemaking of cerebellar Golgi cells observed in vitro  can be explained by a computational model based on a mechanistic reconstruction of the underlying ionic processes. The model predicted that pacemaking depended on the actions of four ionic currents, I h , I Na-p , I K-AHP , and I K-slow : I h brought membrane potential into the pacemaker region where the I Na-p /I K-AHP /I K-slow interaction generated pacemaking. Moreover, following hyperpolarization, I h and I Ca-LVA caused rebound excitation. In response to a depolarization, I K-A helped setting the delay for spike activation and I Na-r (Afshari et al., 2004;Magistretti et al., 2006) enhanced the initial doublet, whose separation from the rest of the discharge was sharpened by adaptation caused by I K-AHP and I K-slow . Calcium-dependent regulation of ionic currents played also a critical role. By being coupled to I K-C , I Ca-HVA enhanced the fast phase of spike AHP Figure 5. Pharmacology of pacemaking. The model predicted the effect of ionic channel blockage  on autorhythmic firing. (A) In these simulations, the control trace shows autorhythmic firing (same as in Figure 1A, firing rate 5 Hz, CV ISI ¼ 0.0027), a 50% block of I h (mimicking the effect of ZD7288) reduces the pacemaker frequency by 27%, and the partial block of I Na-p prevents spontaneous firing. (B) In order to simulate retigabine effects on I K-slow reported by Forti et al., (2006), I K-slow was modified by shifting its activation by À6.5 mV (from Tatulian et al., 2001). The model was then stimulated with 1-second current steps (100 pA) from -70 mV. The traces show that, following I K-slow changes, the model fires less action potentials than in control. The model response to current injection is shown for control (full squares) and I K-slow -shift conditions (empty squares). (C) The effect of membrane noise on firing was simulated by injecting Gaussian noise into the canonical model. In control conditions, the mean firing rate was unaffected, but CV ISI was increased from 0.002 to 0.16. When I K-AHP was blocked, the mean firing rate did not significantly change (from 5 to 5.7 Hz, p > 0.05, t-test). Conversely, CV ISI changed from 0.16 to 0.36 (cf. Forti et al., 2006). The traces show that pacemaker irregularity is markedly increased by I K-AHP blockage (the dashed lines indicate the À70 mV level).
thereby resetting the spiking mechanism and sustaining high-frequency discharge. By being coupled to I K-AHP , I Ca-HVA enhanced the slow phase of spike AHP and spike-frequency adaptation during repetitive discharge.
The numerous currents revealed experimentally allowed the model to simultaneously control multiple aspects of electroresponsiveness, ranging from pacemaking to various response patterns elicited by depolarization and hyperpolarization. Each current could be attributed a characteristic role, as further explained in the following paper (Solinas et al., 2007), despite the apparent redundancy of some mechanisms (e.g., of I K-AHP and I K-slow , which are both slow-repolarizing currents, or I Na-p and I h , which could both sustain pacemaking). Moreover, the model allowed to define a typical, or ''canonical,'' Golgi cell electroresponsive pattern covering an extended experimental set of recordings . Apparently, histochemical differences (Geurts et al., 2001;Simat et al., 2007) did not remarkably impact on basic electroresponsiveness. Although the present model could be improved by more detailed analysis of channel gating and calcium dynamics in Golgi cells, it is a clear advancement over previously published models (Maex and De Schutter, 1998). Pacemaking in the model was sustained by subthreshold oscillations caused by I Na-p coupling with slow repolarizing currents, resembling mechanisms observed in some basal ganglia neurons (see Surmeier et al., 2005 for review). The I Ca-HVA -I K-AHP system influenced the oscillations but, unlike in a model of dopaminergic neurons of the substantia nigra , it was not essential for pacemaking. Although a careful experimental characterization of Na þ and Ca 2þ currents in Golgi  Effect of I Na-r on the first two spikes elicited by a depolarizing current step. The model was hyperpolarized to À70 mV by tonic current injection and a step current injection of 600 pA elicited an initial firing rate of 200 Hz. I Na-r block (middle trace) increased the first ISI (2 ms increase) reducing the initial firing rate to 142 Hz. I Na-r did not affect the 1st spike latency but its activation after the 1st spike could accelerate membrane depolarization and generation of the 2nd spike. (B) The model was hyperpolarized to À70 mV by tonic current injection and then injected with a step current of 100 pA. The 1st spike latency was 18 ms and was reduced by I KA blockage to 12 ms. The lower traces show the stimulus protocols. cells remains to be done, we point out that pharmacological block of SK channels generating I K-AHP did not arrest pacemaking in vitro . I h in our model proved critical to keep the oscillatory mechanisms in the appropriate membrane potential range for activation. This does not exclude that specific regulation of I h reversal potential might enable I h to take more actively part to oscillatory dynamics, as proposed for subthreshold oscillations in entorhinal neurons (Dickson et al., 2000a) and for pacemaking in cholinergic interneurons (Bennett et al., 2000). Finally, I h and I Ca-LVA did not give raise to any low-threshold pacemaking McCormick and Huguenard, 1992;McCormick and Pape, 1990) in Golgi cells but rather played a substantial role in rebound excitation, as also proposed for subthalamic neurons (Chan et al., 2004;Pennartz et al., 1997).
It may appear surprising that Golgi cell pacemaking was not associated with evident subthreshold oscillations, as it has been reported in, for example, amygdaloid neurons (Pape and Driesang, 1998), enthorinal cortex neurons (Alonso and Lliná s, 1989;Dickson et al., 2000a) or inferior olivary neurons (Lliná s and Yarom, 1986). The model predicted that subthreshold oscillations in Golgi cells are tightly coupled to spikes, resembling the case of subthalamic neurons (Bevan and Wilson, 1999). Since pacemaking in the model maintained a similar frequency with or without spikes, the oscillation controlled the pace. Coupling of pacemaker oscillations to spikes was reinforced and regularized by I K-AHP . This effect became clear after introducing noise in the model to mimic in vitro recordings, which are affected by stochastic channel gating. In the noisy regime, pacemaking remained rather regular until I K-AHP was switched off, thereby imitating the effect of I K-AHP blockage by apamine , for a similar effect see also subthalamic nuclear neurons: Bevan and Wilson, 1999;Hallworth et al., 2003). The regularizing action of I K-AHP occurs therefore through a repeated phase-reset of the pacemaker cycle, a concept fully developed in the continuation of this work (Solinas et al., 2007).
The intrinsic excitable properties of Golgi cells could have relevant implications for cerebellar signal processing. The basal discharge may become irregular in vivo due to synaptic noise suggesting that, similar to what was observed in a Purkinje cell model (De Schutter and Bower, 1994;Jaeger et al., 1997, see also Shin et al., 2007), in vivo Golgi cell spiking is still sustained by intrinsic pacemaking (see Solinas et al., 2007). The Golgi cell responses to protracted depolarizing pulses in vitro is characterized (after the initial burst) by a sustained discharge with a frequency depending on the intensity of the input. The modulation of response frequency would allow the Golgi cell to follow the temporal evolution of sensory-motor patterns (Miles et al., 1980). This confirms that the shortlasting Golgi cells spiking responses to long duration stimuli observed in vivo (Tahon et al., 2005) are due to depression of afferent input, as opposed to intrinsic mechanisms. Moreover, the presence of mechanisms enhancing the initial spike burst and regulating its initiation may be relevant for exerting precisely timed inhibitory effects on granule cells. An extended investigation of the effect of the ionic channel complement is presented in the companion paper (Solinas et al., 2007) showing that Figure 8. Robustness of the model. Tests of robustness of the model against changes in individual parameters. Polar plot: The maximum conductance of I K-slow , I Na-p , I KA , I K-AHP , I Ca-HVA , and I h , as well as b Ca , were modified in turn. The changes covered (in 30 steps) the range from total block (center of the polar plot) to threefold increase (outer solid circle) in the canonical model value (thick gray line). Model simulations were obtained using the protocols described in Figures 2-4. The thin black lines delimiting the light gray area show the minimum and maximum parameter changes (percentage of canonical model values), which maintained the model within the boundaries of Golgi cell characteristic response. Lower panels: Results were considered to conform to Golgi cell behavior on the base of five criteria. (i) The pacemaker frequency had to fall in the range 0.5-9 Hz. (ii) The first spike latency triggered by the rebound excitation had to fall in the range 30-100 ms. (iii-v) The response to current injection had to show first spike latency, initial firing rate, and spike frequency adaptation within the boundaries reported in the panels. Thin black lines are boundaries derived from the envelope of data points in Figure 4B and the thick gray line shows the canonical model behavior.
I K-AHP and I K-slow appear primarily involved in causing phase-reset and resonance, respectively. The model has therefore the ability to predict Golgi cell responses in various operating condition and appears as a useful tool for large-scale simulations of the cerebellar network. Figure 9. Subthreshold oscillations and spike coupling. In order to investigate the mechanisms of subthreshold oscillations, simulations were performed after blocking the currents generating the spike (I Na-t , I K-V , and I K-C ). (A) The simulated traces show that preventing action potentials lead to the emergence of subthreshold pacemaker oscillations slightly slower than in the spiking regime (4.46 vs. 5 Hz). Upon 100% I h blockage, the oscillation frequency is reduced (3.3 Hz) but the effect is reversed with the injection of a tonic depolarizing current (dashed trace). (B) The simulated traces (all voltage-dependent currents blocked except for I Na-p , I K-AHP , I Ca-HVA , and I K-slow ) show that subthreshold pacemaker oscillations are interrupted by blocking I Na-p and that only marginal recovery is obtained by injecting a tonic depolarizing current (dashed trace). Either blocking I K-AHP or I K-slow deforms the oscillation trajectory, while their simultaneous block causes the oscillation to cease. (C) The voltage and current traces are shown as a function of time during a full oscillation cycle (all voltage-dependent currents blocked except for I Na-p , I K-AHP , I Ca-HVA , and I K-slow ). (D) Schematic showing the cascade of processes that generate pacemaking in the Golgi model.