A multi-compartmental model representative of the main morphological components of a UBC was constructed using the NEURON simulator (NEURON version 7.3; Hines and Carnevale, 2001). The model consisted of soma, dendrite, initial segment and axon compartments generating a morpho-electrical equivalent of the UBC (Figure 1A). The voltage- and Ca²⁺-dependent mechanisms (see below) were distributed among the compartments. With this approach, the model reproduced satisfactorily basic aspects of UBC electroresponsiveness elicited by somatic current injection. It should also be noted that, although the existence of UBC subtypes has been suggested based on histochemical analysis, the basic electroresponsive properties were homogeneous in a large majority of UBCs (Locatelli et al., 2013). Thus, we have reconstructed a canonical UBC model, which simulates the typical electrophysiological behavior of UBCs and lies within the scatter of physiological parameters values. The model was coupled to a recording electrode to reproduce realistic current-clamp conditions and was calibrated on the cellular response corrected for −10 mV liquid-junction potential.

The model ionic conductances were non-homogenously distributed over different compartments (Figure 1B; Table 1). Gating kinetics were corrected using a Q₁₀ = 5 for I_CaLVA activation (Destexhe et al., 1994) and a Q₁₀ = 3 for all the other currents according to the relation Q^{(Tsim−Texp)/10}₁₀ (Gutfreund et al., 1995) to account for temperature differences between experimental recordings and the model (30°C). Nernst equilibrium potentials were pre-calculated from ionic concentrations used in current-clamp recordings and maintained fixed, except for the Ca²⁺ equilibrium potential, which was updated during simulation according to the Goldman-Hodgkin-Katz (GHK) equation. The maximum ionic conductances were regulated to fit the UBC responses to various stimulations (Traub and Llinás, 1979; Traub et al., 1991; Vanier and Bower, 1999; Achard and De Schutter, 2006; Druckmann et al., 2007, 2008, 2011, 2013; Solinas et al., 2007a,b). Ionic channel gating was modeled following the Hodgkin and Huxley formulation or Markov-chains for multi-state transitions using mathematical methods reported previously (D'Angelo et al., 2001; Nieus et al., 2006; Solinas et al., 2007a,b) and are reported in Table 1. Membrane voltage was obtained as the time integral of the equation (Yamada et al., 1989).

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. The intracellular Ca²⁺ concentration, [Ca²⁺], was calculated through the equation, (2)d[Ca2+]dt=−ICa(2F ∗ Ad)−βCa([Ca2+]−[Ca2+]0)

Where d is the depth of a shell adjacent to the cell surface of area α, β_Ca determines the loss of Ca²⁺ ions from the shell approximating the effect of fluxes, ionic pumps, diffusion, and buffers (Traub and Llinás, 1979; McCormick and Huguenard, 1992; De Schutter and Smolen, 1998), and [Ca²⁺]₀ is resting [Ca²⁺].

UBC recordings were carried out as reported previously (Locatelli et al., 2013). Whole-cell patch-clamp recordings were performed from 30 UBCs in the internal granular layer of lobule X in acute rat cerebellar slices of P18-P25 wistar rats. Out of these, 12 UBCs were used to study intrinsic electroresponsiveness and 18 to study the LOR. In addition, to estimate the passive properties and the resting membrane potential, data were also re-analyzed from a population of 51 UBCs reported by Locatelli et al. (2013) for a total of 81 measurements.

In brief, the patch pipettes were pulled with a horizontal puller (Sutter Instruments, Novato, CA, USA) from thick-walled borosilicate glass capillaries (Hingelberg, Germany) and had a resistance of 7–10 MΩ when filled with the intracellular solution (in mM): 126 potassium gluconate, 4 NaCl, 15 glucose, 5 Hepes, 1 MgSO4.7H2O, 0.1 BAPTA, 3 ATP, 100 μm GTP; pH adjusted to 7.2 with KOH (in solution at 0.1 mM). The BAPTA-Ca²⁺ buffer was prepared as explained previously (D'Angelo et al., 1995; Gall et al., 2003). To minimize pipette tip capacitance, pipettes were coated with Sylgard (Dow Chemical, Midland, MI, USA) and the bath level was kept as low as possible. After careful pipette capacitance cancelation (D'Angelo et al., 1995), −10 mV 250 ms voltage step were applied to the cell from the holding potential of −70 mV and the corresponding current transients were recorded (sampling frequency at 20–40 kHz, low-pass filtering 5 kHz). The current transient decayed multi-exponentially, probably reflecting charging of compartments corresponding to the UBC soma, dendritic shaft, brush, and axon. Tri-exponential fittings allowed us to estimate the electrode series resistance (R_s) the cell input capacitance (C_in) and the input resistance (R_in) using the equation: (9)I(t)=A1Xexp(k−t/τ1)+A2Xexp(k−t/τ2)             + A3Xexp(k−t/τ3)+C

Series resistance (R_s) was calculated as R_s = τ_VC/C_in, where τ_VC = τ₁ is the decay time constant of the current transient related to somatic charging. In our extended UBC sample, we measured C_in = 16.9 ± 0.67 pF, R_in = 0.87 ± 0.07 GΩ, R_s = 22 ± 0.96 MΩ (mean ± MSE, n = 81 for all parameters).

All current-clamp recordings were performed using the fast current-clamp mode of the amplifier to accelerate membrane charging (D'Angelo et al., 1995, 1998; Prestori et al., 2008). The resting membrane potential was measured over 50 ms in 12 consecutive traces. Intrinsic excitability was investigated by setting resting membrane potential at −80 mV and injecting 800 ms current steps (from −16 to 48 pA in 4–8 pA increments). Current steps were applied from the holding potential of −80 mV and the corresponding voltage deflections were recorded (sampling frequency at 20–40 kHz, low-pass filtering 5 kHz).

Electrical stimulation was performed by placing a bipolar tungsten electrode over the mossy fiber bundle in the granular layer (D'Angelo et al., 1993, 1995). The stimuli consisted of voltage pulses (0.25 ms, 5–15 V) organized in 100 Hz trains of varying intensity and length.

UBCs (Figure 1A) were recorded from the granular layer of lamina X of the rat vestibular cerebellum (Mugnaini and Floris, 1994) using patch-clamp techniques in acute slices and were modeled using the NEURON simulator. In order to investigate whether the proposed LOR mechanisms were consistent with current knowledge on UBC intrinsic electroresponsiveness and whether the prediction of LOR generating mechanisms was correct (Locatelli et al., 2013), we have developed a detailed UBC computational model, in which synaptic transmission was coupled to H- and TRP-channels through a cytoplasmic modulatory mechanism.

The model was constructed respecting the morphological constraints and passive properties (Figure 1) and using the available electrophysiological data on ionic channels of UBCs (Table 1). Then the model was calibrated through a multi-parametric comparison with current-clamp response to current injection (Figures 2, 3). Finally, the model was extended through a mechanism coupling synaptic activity to ionic channel gating and became able to generate the LOR (Figure 4). The robustness of this procedure was evaluated by performing various tests allowing model predictions to be compared to native UBC responses (Figures 5–8).

The UBC model (Figure 1B) was composed of 4 main compartments (brush, dendritic shaft, soma, axonal initial segment, axon). The model was endowed with 8 voltage-gated or Ca²⁺-gated ionic currents (I_Na, I_KV, I_KA, I_Kslow, I_KCa, I_CaLVA, I_CaHVA, I_H) and with a voltage-independent leakage current (I_TRP). The kinetics of individual ionic currents were fitted to UBC data (when directly available) or adapted from previous cerebellar neuron models (D'Angelo et al., 2001; Anwar et al., 2012). The density of ionic channels in different model compartments was pre-set according to the known channel localization (e.g., Na channels were more expressed in the AIS, CaLVA channels were more abundant in the brush than in the soma while CaHVA channels were equally distributed in soma and brush; (Khaliq and Raman, 2006; Diana et al., 2007) and subsequent adjustments to the current density were made to properly fit the UBC response to current injection (cf. Solinas et al., 2007a,b; see below).

The calibration of UBC passive properties was done with respect to electrophysiological measurements (Figure 1C). After having reconstructed the UBC compartmental morphology, the injection of voltage step into the model yielded input capacitance and resistance values falling in the middle of the experimental data distribution (the average experimental values were 16.3 pF and 0.97 GΩ, the average model values were 16.7 pF and 1 GΩ).

In current-clamp whole-cell recordings, UBCs showed a typical electroresponsive pattern comprising the following elements (Figure 1A):

UBCs were silent at rest, with a resting potential of −66.7 ± 1.8 mV (n = 81) (cf. −67.6 ± 0.9 mV in Diana et al., 2007; −67.5 ± 1.4 mV in Locatelli et al., 2013).

These electroresponsive properties were common to those reported in the previous experimental investigations (Diana et al., 2007; Birnstiel et al., 2009; Locatelli et al., 2013), in which the implication of specific ionic channels was examined (see below).

Following mossy fiber bundle stimulation in the presence of GABA and glutamate receptor antagonists, UBCs showed a LOR composed of a depolarizing ramp followed by a LTS surmounted by a high-frequency spike burst followed by a prolonged tonic discharge (Figure 2B1). The LOR corresponded to the same phenomenon reported by Locatelli et al. (2013), who suggested that the LOR depended on a metabotropic mechanism involving cAMP increase in the cytoplasm, leading to the H- and TRP-channel modulation. Notably, the delay, duration and frequency of the LOR were apparently related to the intensity and number of pulses used for mossy fiber bundle stimulation (Figure 2B1). Higher intensity or number of pulses caused shorter delay, higher frequency, longer duration of the LOR (quantitative data are given below in Figures 6, 7).

The UBC model was parameterized against responses to step-current injection (Figure 2B1). The parameters related to spikes elicited by depolarizing current steps (e.g., using 16 pA from a holding potentials of −80 mV) included spike delay (30.3 ms), spike overshoot (36.7 mV), instantaneous frequency (176 Hz), steady state frequency (126 Hz), and frequency ratio (0.715). The parameters related to hyperpolarizing current steps included sag amplitude (e.g., −17.3 mV with a −20 pA step from −80 mV) and current threshold for rebound bursting (e.g., −20 pA from −80 mV). The match of model to the experimental data was statistically significant (Student's t-test, p < 0.01) for all the parameters. This model, once endowed with coupling mechanisms regulating TRP and H channels, predicted the LOR following synaptic stimulation (Figure 2B2).

It should be noted that what has been developed in this paper is a canonical model matching the properties of an ideal UBC. Variants in terms of burst length, frequency ratio, or rebound depolarization observed in some UBCs could all be obtained through slight changes of maximum current densities for specific ionic channels (data not shown).

In the model, as in real UBCs, a small depolarizing current (8–16 pA) from a negative potential (−80 mV) was sufficient to trigger the LTS and a short spike burst, while a larger depolarizing current (>16 pA) was required to drive the model into tonic discharge (Figures 2B1, 3A; cf. Locatelli et al., 2013). However, when depolarizing step-current injection was applied from resting potential (−67.5 mV; Diana et al., 2007; Locatelli et al., 2013; Figure 2B1, Left), the response was characterized by a tonic discharge almost lacking the initial spike bust (Figure 2B1, Right). Even by using a net current injection stronger than that used from −80 mV, the instantaneous spike frequency was apparently lower. While a classical set of currents—including I_Na, I_KV, I_KA, I_CaHVA, and I_KCa—was sufficient to determine the spike shape as well as spike frequency and its variations (D'Angelo et al., 2001), specific ionic channels were required to determine these characteristic properties of UBC electroresponsiveness.

The high responsiveness of the UBC model to small current injection from negative holding potentials was due not only to the relatively high input resistance but also to activation of the CaLVA current. CaLVA was the most prominent current in this functional regimen and regeneratively amplified the effect of current injection causing the LTS. CaLVA currents have a voltage-dependent inactivation, reflected into the marked weakening of bursting elicited from more depolarizes membrane potentials (Diana et al., 2007; Birnstiel et al., 2009; Locatelli et al., 2013). The density of CaLVA currents was derived from voltage-clamp recordings (Diana et al., 2007) and, together with a precise setting of input resistance, provided a major constraint for calibrating the whole electroresponsive mechanism. Once the LTS was generated, the depolarization gated the other ionic channels so that various currents contributed to generate the spike bursts and control its evolution. As in real UBCs, burst spike frequency could raise up to 200 Hz along with spike amplitude adaptation. The burst was probably enhanced by the resurgent Na current (Raman and Bean, 1997; Afshari et al., 2004; Dover et al., 2010) and the first spike delay was regulated by I_KA.

With depolarizing current injections <16 pA, the burst was driven by the LTS and terminated after a time compatible with CaLVA channel inactivation. With depolarizing current injections >16 pA, spike discharge protracted beyond the burst through the intervention of another set of ionic channels: the CaHVA current and the persistent Na current sustained a depolarizing plateau maintaining spike discharge. The model predicted that an appropriate balance between the CaHVA and a slow K current, Kslow, was required to regulate the intensity and duration of this protracted spike discharge. Finally, during depolarization, the H current was deactivated and therefore reduced the total input conductance favoring depolarization, but its contribution to the total UBC current was limited by the decreased driving force.

As in real UBCs, when injected with hyperpolarizing current steps, the model generated sagging inward rectification followed by rebound excitation after the end of the pulse. The rebound could be consisted of a depolarizing ramp that could prime a LTS and a spike burst when the hyperpolarizing was more intense and protracted (Figure 3C). The rebound was driven by H current deactivation and amplified by CaLVA current activation (see Figure 6; Russo et al., 2007; Locatelli et al., 2013). The H current was derived from voltage-clamp recordings (Locatelli et al., 2013) and provided a critical constraint for model calibration.

By using the model, we tested the hypothesis that synaptic modulation of H- and TRP-currents through an intracellular cascade was necessary and sufficient to generate the LOR (Locatelli et al., 2013). Therefore, if the UBC model was correctly constructed, it should be able to generate the LOR by coupling synaptic transmission to H- and TRP-current gating through an appropriate second messenger pathway. The second messenger pathway was implemented as a generic reaction scheme, in which (1) mossy fiber bundle stimulation activated the extra-synaptic receptor X to X^*, (2) X^* increase above 60% of its maximum value initiated conversion of a second messenger Y to Y^*, and (3) Y^* shifted the voltage-dependence of I_H and the maximum conductance of I_TRP (see Methods; DiFrancesco and Mangoni, 1994; Locatelli et al., 2013; Figure 4). The increased Y^* concentration eventually caused a proportionate increase in TRP maximum conductance and a proportionate shift of I_H activation toward positive membrane potentials. TRP depolarized the membrane sufficiently to generate the LOR, while I_H transiently unbalanced the sub-threshold current generating a slow depolarizing ramp tuning the delay of LOR initiation (Figure 4). It is worth noting that, after setting model intrinsic electroresponsiveness, no further changes was required to obtain the LOR, which is therefore a true model prediction.

The model could reproduce the fundamental properties of the LOR, which are observed in WCR recordings (this work and Locatelli et al., 2013). (i) The LOR was dependent on the intensity and number of synaptic pulses (cf. Figure 2). (ii) The TRP current generated the major depolarizing drive for LOR generation, while the H-current regulated LOR delay by controlling the slope of the depolarizing ramp. Blocking H- and TRP-currents together abolished the LOR (Figures 4, 5). (iii) The Ca²⁺ currents enhanced LOR spike frequency, with a more evident contribution given by LVA than HVA currents. However, Ca²⁺ currents did not regulate the depolarizing ramp nor they were required to generate the LOR itself (Figure 5). (iv) The LOR was poorly affected by holding potential and could be generated either from −80 or −67.5 mV (data not shown). The reliability of model predictions about UBC electroresponsiveness and LOR generation were assessed through a quantitative comparison with experimental data.

I_H provided a critical constraint for the model and the effect of pharmacological block by ZD7288 and Cs⁺ has been reported (Diana et al., 2007; Russo et al., 2007; Locatelli et al., 2013). I_H switch-off in the model caused hyperpolarization of resting membrane potential (Figure 5) and eliminated the sag during hyperpolarizing current steps (Figure 5, inset). The resting membrane potential shift (−4 mV) and differences between the maximum voltage deflection and the voltage attained at the end of the current step (from 3.96 to 0.13 mV) and the subsequent rebound depolarization (from 2.54 to 0.15 mV) were consistent with previous experimental data (Locatelli et al., 2013). Finally, the LOR delay increased by 12.8% and the slope of the depolarizing ramp and its duration decreased by 12.2 and 24.9%, respectively, again reflecting experimental data (Locatelli et al., 2013). Therefore, in the model the I_H regulated the delay of the late-onset burst by controlling the depolarizing ramp similar as in whole-cell recordings.

The TRP channels provided a background leakage current, which increased during the LOR. I_TRP switch-off in the model prevented full LOR generation leaving only a subthreshold depolarization driven by the H current (Figure 5). The slope of the depolarizing ramp decreased by 69%, consistent with 49% decrease measured experimentally following channel blockage with SKF96365 (Locatelli et al., 2013).

Pharmacological blockage of CaLVA channels (with mibefradil) and CaHVA channels (with nimodipine) was reported to modify spike firing in the LOR but not to prevent LOR generation (Locatelli et al., 2013). In the model, CaLVA switch-off did not prevent the LOR but reduced the number of spikes by 20% and the instantaneous frequency by 38%. Combined CaHVA and CaLVA currents switch-off did not prevent the LOR but reduced the number of spikes by 80% and the firing frequency by 88% (Figure 5). These effects in the model are comparable to those of pharmacological blockage reported experimentally.

As explained above (see also Methods), reference values were obtained from UBC voltage-clamp data on the H- LVA- HVA- and TRP-currents and were then used to pre-set the corresponding maximum conductance values in the model. The other maximum conductance values were set respecting proportions evaluated in previous models [for example the Gmax (I_Na)/G_max (I_KV) ratio, the G_max (I_KCa)/G_max(I_CaHVA] ratio, (Gabbiani et al., 1994; D'Angelo et al., 2001; Solinas et al., 2007b). Nonetheless, the maximum ionic conductances in the model were actually free parameters, as they were tuned to match UBC electroresponsiveness. A robustness analysis was then performed in order to estimate the confidence intervals of the maximum ionic conductances still allowing to generate existing UBCs. The range of variations allowed for each maximum conductance were used to maintain the model within the experimental values of first spike delay, instantaneous frequency, steady-state frequency, frequency ratio and spike overshoot (Figure 6). Finally the model predicted, through the cytoplasmic coupling mechanism, the experimental values measured in the LOR for the parameters indicated above as well as for LOR duration.

The experimental parameters measured in the UBCs covered in this study were used to set the range limits for parameter variation in the model. Intrinsic excitability was characterized by spike delay (31.1–93.4 ms, n = 12), spike overshoot (15.0–45.1 mV, n = 12), instantaneous frequency (95.9–202.7 Hz, n = 12), steady-state frequency (41.2–141.8 Hz, n = 12), and frequency ratio (0.61–1.3, n = 12). The LOR was characterized by spike delay (98–1000 ms, n = 18), spike overshoot (10.2–44.5 mV, n = 18), instantaneous frequency (92–209 Hz, n = 18), steady-state frequency (19.5–60 Hz, n = 18), frequency ratio (0.1–0.6, n = 18), and burst duration (100–2500 ms, n = 18). Simulation results were accepted as possible UBC responses when falling within these limits.

Both concerning responses to current injection (Figure 6A) and the LOR (Figure 6B), the modulation of maximum ionic conductances led within boundaries spanning ≤ ± 20%. By combining the two sets of boundaries obtained for current injection and LOR gave the compound boundaries (Figure 6C). This robustness analysis showed that the canonical model was well balanced and could use the same ionic channel settings to faithfully explain both electroresponsiveness and LOR.

The model allowed to predict burst delay, duration and frequency of LOR in relation to different combinations of stimulus duration and intensity (Figure 7A). The response showed opposite changes for delay (which decreased with stimulus duration and intensity) compared to duration and frequency (which increased with stimulus duration and intensity). If the simulated response space was correct, the experimental data points obtained changing stimulus duration should actually correspond to the same stimulus intensity. Actually, for all parameters, the experimental data points fell on nearly horizontal lines (Figure 7A). A quantitative assessment of the matching between experiments and model was obtained by calculating the slope of the duration/intensity relationships (Figure 7B) and by assessing the variability of the mean predicted stimulus intensity (Figure 7C). The slope of the duration/intensity relationships was 0.0005 ± 0.002 and was statistically indistinguishable from 0 (p < 0.001, n = 7, t-test). The variability of the predicted stimulus intensity was <10%. Therefore, the model could effectively predict LOR properties while varying the stimulus pattern.

The explanation of why, in the LOR stimulus/response space, similar parameter values were determined by varying either stimulus duration or intensity (cf. Figure 8A) had to be searched in their control over the common factor, Y^*. In turn, Y^* controlled membrane channel gating and therefore the LOR (Figure 8B).

The effect of Y^* on LOR parameters (burst delay, duration and frequency) for all possible combinations of stimulus duration and intensity used in Figure 7 are shown in Figure 8A. In particular, while increasing Y^*, burst delay decreased almost exponentially, burst frequency increased almost exponentially, and duration increased almost linearly (Figure 8A; the presence of multiple points at each Y^* value reflected the different Y^* kinetics for different combinations of stimulus intensity and duration).

A mechanistic explanation of the LOR properties was obtained by correlating Y^* with the underlying ionic currents (Figure 8B). It turned out that the Y^* dependence of I_Na, I_KV, I_Kslow, I_CaHVA, I_KCa, and I_TRP was similar to that of LOR duration. The Y^* dependence of I_H and I_KA was similar to that of LOR frequency and opposite to that of LOR delay. The Y^* dependence of I_CaLVA was inverted compared to any others. This result implies a complex interaction of all currents in regulating LOR properties and a primary role of I_H in controlling LOR delay.

The present model provides a biologically realistic reconstruction of UBC electroresponsiveness (Diana et al., 2007; Russo et al., 2007; Birnstiel et al., 2009; Locatelli et al., 2013) based on the representation of neuronal geometry, passive properties and on the subsequent insertion of appropriate ionic and intracellular mechanisms. The H, TRP, and Ca²⁺ conductances provided the critical constraints for UBC modeling and had density, gating and kinetics derived from UBC recordings (Diana et al., 2007; Russo et al., 2007; Birnstiel et al., 2009; Locatelli et al., 2013). These core mechanisms were associated with a spike generation mechanism sustained by Na⁺ and K⁺ conductances. The CaLVA and CaHVA channels were distributed across soma, dendrite and brush in order to generate calcium currents proportional to those determined experimentally (Diana et al., 2007; Birnstiel et al., 2009), while the Na channels were localized in soma, initial segment and axon. The localization of ionic channels had an impact on model electrogenesis. On one side, the CaLVA current was larger in the brush than soma, while the CaHVA current was equally distributed, as indicated by calcium current measurements (Diana et al., 2007; cf. Figure 3C). The high density of CaLVA channel in the brush, coupled with limited current flow to the soma through the dendrite, enhanced the local auto-regenerative LTS depolarization. This would become particularly effective when synaptic currents are injected into the brush. On the other side, the Na⁺ current was concentrated in the soma, initial segment and axon, conferring higher excitability and normalizing spike and burst discharge parameters (not shown). The precise localization of Na⁺ channels in UBCs remains to be determined by immunohistochemistry.

The resulting model was a canonical (or prototypical) representation of the UBC, such that ±20% variations of ionic channel maximum conductances still allowed model responses (first spike delay, initial firing rate, steady-stated firing rate, frequency ratio, spike overshoot) to remain within the range of parameters recorded experimentally. This variation range is normally observed while measuring ionic currents in neurons (Bhalla and Bower, 1993) and has also been reported in UBCs for H-, CaLVA, and CaHVA currents (Diana et al., 2007; Birnstiel et al., 2009; Locatelli et al., 2013). Most importantly, the model allowed to predict the LOR without any further modifications, provided that the cytoplasmic coupling mechanism was properly implemented.

Beside its effectiveness, some aspects of the UBC model were admittedly simplified. (i) The H-current kinetics observed in UBC recordings (Locatelli et al., 2013) suggest that H-current gating could reflect differential expression and assembly of subunits. This could contribute to tune UBCs toward specific LOR kinetics. (ii) The Na current was implemented using a 13-state stochastic model for Nav1.6 channels of granule cells generating transient, persistent and resurgent current components (Raman and Bean, 2001; Magistretti et al., 2006; Dover et al., 2010). While these three components are present in all cerebellar neurons (Afshari et al., 2004; Khaliq and Raman, 2006), the resurgent current may be proportionally higher in UBCs than granule cells (Afshari et al., 2004). This could further promote the initial spike burst primed by the LTS. (iii) The regulation of duration and frequency of tonic discharge following the LTS required that a slow K current, Kslow(D'Angelo et al., 2001), was added to the model. Moreover, the regulation of first-spike delay required that a fast-activating K current, KA(D'Angelo et al., 2001), was added to the model. These current have not been reported in UBCs yet and are therefore model predictions, which remain to be proved experimentally. (iv) Intracellular calcium regulation was modeled using a simplified scheme useful to couple calcium entering through CaHVA channels to BK-type channel gating and the effect of more complex schemes remains to be explored.

In a further evolution, the model may benefit of a detailed morphological reconstruction of UBC brush and axon along with proteomic and channelomic analysis directly defining the molecular properties of ionic channel subtypes and of the cytoplasmic cascade coupling membrane receptors to ionic channels. It will also be of interest to determine under which circumstances routines based on genetic algorithms will be able to solve the complex optimization problem of the UBC model (Druckmann et al., 2007, 2008, 2011, 2013).

A quite relevant property of the UBC model, as well as of real UBCs, was that of generating a variety of burst responses, which were controlled by the intervention of H-, TRP-, and Ca²⁺- currents. The model showed that the interaction of H- TRP-, and Ca²⁺- currents emerged in three conditions.

These mechanisms, especially the first two, proved markedly voltage-dependent due to the specific voltage-dependent inactivation of CaLVA channels. It is therefore predicted that alternating cycles of depolarization and hyperpolarization involving the granular layer inhibitory circuit (Vos et al., 1999) could modify the ability of the UBCs to generate spike bursts substantially.

Rather than reproducing the intracellular biochemical pathway (which remains largely undetermined) we have developed a generic coupling mechanism including receptors activated by synaptic activity and leading to the production of a second messenger (Y^*), possibly corresponding to ATP/cAMP conversion. The level of this second messenger caused the current influx (~-16 pA) generating the depolarizing ramp and the LOR, with a delay depending on the second messenger concentration and ranging up to several hundreds of milliseconds.

Thus, the model provided substantial support to the existence of an intracellular mechanism coupling membrane receptor activation to H- and TRP-channels and causing the LOR in UBCs. Modeling suggests that I_H was responsible for controlling LOR delay and frequency, that I_CaLVA and I_KA were driven into the control of LOR delay and frequency through their time- and voltage-dependent inactivation, that I_Na, I_CaHVA and I_TRP currents drove the depolarization during the burst and determined burst duration and activation of I_KV, I_Kslow, and I_KCa, which actually regulated the frequency and terminated the burst.

Through the LOR, UBCs can translate the intensity of the input (coded as the number of active fibers and duration of their discharge) into output bursts with different delay, duration and frequency. This implements a time-code that could reverberate through the network along the chains made by UBCs with other UBCs and granule cells (Nunzi and Mugnaini, 2000; Nunzi et al., 2001). The LOR could be modulated by local network activity causing various patterns of UBC excitation and inhibition (Ruigrok et al., 2011; Rousseau et al., 2012). Recently, a frequency-dependent regulation of AMPA receptor-mediated after responses has also been reported (van Dorp and De Zeeuw, 2014) and has been suggested to provide a system to generate protracted responses in UBCs. The implementation of these AMPA receptor-dependent mechanisms into the UBC model would allow to test its actual impact on UBC responses and its relationship with the LOR.

It has been proposed that UBCs are essential for shifting and converting the phase of mossy fiber activity that relays information from the vestibular apparatus, eyes, or neck, and that their characteristic cellular properties are particularly relevant for controlling and consolidating motor learning in paradigms such as VOR phase reversal (Gao et al., 2012). Actually, the LOR could help implementing the “velocity storage” system (Dai et al., 2007), which allows the vestibulo-cerebellum to transform the head velocity signal into commands for ocular motor neurons controlling the slow phase of nystagmus (Goldberg et al., 1984; Okada et al., 1999). Patterns of activity compatible with the intervention of the LOR in UBCs have been reported in the electric organ (a cerebellar-like structure) of mormyrid fishes during the prediction of the sensory consequences of motor acts (Kennedy et al., 2014). By involving intracellular cascades, the LOR could be modified by receptor systems—involving e.g., acetylcholine, serotonin, and noradrenaline—known to modulate H and TRP channels, thus correlating cerebellar computation with the brain functional state.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.