Several stimulation protocols were designed to reproduce the basic patterns used experimentally.

The model response was recorded either in voltage clamp or in current-clamp at the soma, and the voltage-dependent ionic conductances were modified in some cases to test their impact on BPR generation. The membrane voltage, the calcium concentration, the ionic currents from each model section were recorded and saved in a nested MATLAB file. The model response properties were analyzed using MATLAB routines (MathWorks, Natick, MA, USA) and custom made Python scripts, using data recorded during each simulations. Each simulation data file contained information about the membrane voltage, ionic channels and synaptic currents recorded in specific location throughout the models.

Biophysical analysis of the model was carried out according to standard theory (Jack et al., 1975). The model electrotonic properties were analyzed using NEURON functions yielding the input impedance, Z_in, for each model section and the corresponding signal attenuation between e.g., a dendritic sections and the soma, A = V_dend/V_soma. This allowed to determine the electrotonic distance L = ln(A). This definition coincided with the more classical one (L = anatomical distance/length constant) defined when a neuron can be reduced through the 3/2 power branching rule to an equivalent cylinder. The NEURON definition proves particularly useful since PCs do not follow the 3/2 power branching rule (Hines and Carnevale, 2001; Hines et al., 2007, 2009).

Since the PC was continuously pacing and input bursts were applied randomly with respect to ongoing spike discharge, repeated PC responses to the inputs differed one from each other. This behavior was represented using raster plots and peri-stimulus time histograms (PSTH), that were constructed using MATLAB routines and were subsequently analyzed to extract the average properties of PC responses.

A first question was whether the PC model was able to generate BPRs following activation of excitatory and inhibitory synapses (Cao et al., 2012; Herzfeld et al., 2015) and whether BPRs were different between Z+ and Z− PCs. The PC responses to input bursts generated by Granule cells (Rancz et al., 2007; van Beugen et al., 2013; Powell et al., 2015; Wilms and Häusser, 2015; Delvendahl and Hallermann, 2016) were tested in Z+ and Z− PC models in three stereotyped functional cases corresponding to the definitions of type I, type II and type III responses reported in vivo (De Zeeuw et al., 2011). Type I responses were driven by pure Granule cell excitatory inputs (100 synapses), type II responses were driven by both granule cell excitatory inputs (100 synapses) and SCs inhibitory inputs (25 synapses), type III responses were driven by SC inhibitory inputs only (25 synapses; Figure 2). In all cases, molecular layer interneuron inhibition was delayed by 4 ms to account for delays accumulated along the synaptic chain (Eccles et al., 1967; Ramakrishnan et al., 2016). In this set of simulations, a random distribution of active synapses was used and maintained the same with all the different stimulation patterns (Figure 2A).

The BPRs mechanism was analyzed by tracking transmembrane currents and intracellular axial currents in different cellular compartments using the Z+ PC type I model, in which the intrinsic ability to generate BPRs was evident (Figure 3A). In the PC model, APs arise from the AIS (Masoli et al., 2015) and, not surprisingly, the pause is characterized by a protracted decrease of the AIS transmembrane current (I_AIS; Jack et al., 1975). In fact, this makes the pause appearing as an interruption of pacemaking. Pacemaking in PCs is sustained by the persistent Na currents generated mostly in AIS and by Ca currents generated mostly in the dendrites and transmitted to the AIS through the internal resistance (Llinás and Sugimori, 1980a,b), providing two candidate mechanisms for pause generation.

Na currents in the AIS may be inactivated during the burst and then take time to recover, reducing the depolarizing drive. However, there was no remarkable inactivation of Na currents during the bursts that might prevent reactivation of the pacemaker. Therefore, this first mechanism could be excluded.

Ca channel activation and the consequent Ca entry into the dendrites may activate KCa currents overtaking inward currents and generating a protracted repolarizing drive. Actually, the current transmitted from the dendrite to soma (I_dend-soma) and from soma to AIS (I_soma-AIS) was reduced during the pause compared to pacemaking regime. Therefore, reduced depolarizing current transmission from the somato-dendritic compartment to AIS appears as the main responsible of the pause.

The ionic nature of the mechanism was confirmed by specific manipulations of the Ca/KCa ionic mechanisms. When either dendritic KCa conductances or Ca conductances were set to zero, the BPR was altered with a marginal reduction of burst duration (from 0% to −16.7%) but a much more dramatic reduction of the pause (from −20% to −80%) (Figure 3B). It should be noted that all Ca channels proved critical for BPR, especially Cav2.1 and Cav3.1 (Cav3.2 switch-off blocked pacemaking; Masoli et al., 2015) as well as all KCa channels, including KCa1.1, KCa2.2, KCa3.1.

These observations show that making BPRs is an intrinsic property of the PC model, which depends on the generation of large KCa currents in the dendrites. In Z− PCs, the mechanism were the same as in Z+ PCs except that TRP channels injected a constant inward current through the dendrites counterbalancing KCa and making the pause more difficult to elicit (not shown).

A central issue in cerebellar physiology is how PCs transform signals coming from granule cells into specific outputs. The spike patterns emitted by granule cells in response to punctuate stimulation consist of short bursts composed of spikes with variable number and ISI (Chadderton et al., 2004; Rancz et al., 2007; D’Angelo and De Zeeuw, 2009). We have therefore used the PC model to simulate the impact on BPR of stereotyped patterns composed by a fixed number of pulses at different frequencies. Concerning burst length, the PC model showed an almost linear input/output relationship at frequencies ranging from 100 Hz to 500 Hz (Figure 4A, left). The pause also showed a slight dependency on the length of the input burst and this eventually caused the burst/pause ratio to reliably report the duration of the input burst (Figure 4A, right). There are two additional noteworthy properties. First, the linearity of burst and burst/pause coding was maintained almost independently from the input frequency. Secondly, the Ranvier nodes in the PC axon filtered the highest frequencies (Masoli et al., 2015) so that the burst was more reliably transmitted at the lower input frequencies (Figure 4A, bottom).

Another modality of granule cell response, largely investigated in the vestibulo-cerebellum, is to generate spike trains in response to prolonged mossy fiber inputs (Arenz et al., 2008). When such stimuli were used, the PC model followed the input frequency (Figure 4B, left), although the frequency dependence was rather weak (as also noted for the BPR in Figure 4A). However, there was a steep relationship between the output frequency and the number of active synapses (Figure 4B, right). Therefore, the PC model was more efficient in detecting the intensity of the granule cell input rather than the frequency at individual synapses.

The recoding of input into output spike patterns in PCs is thought to depend on how long-term synaptic plasticity modifies PC synaptic responsiveness. While classical pf-PC LTD has postsynaptic expression and simply causes a scale-down of postsynaptic currents, a presynaptic change in release probability (P) would modify neurotransmission dynamics (Tsodyks and Markram, 1997; Nieus et al., 2006). These can shift from the typical short-term facilitation at low P to short-term depression at high P. These changes can therefore differentially affect the PC BPR.

At different P values (range 0.1 and 0.9), the overall response pattern in PC output bursts did not change remarkably. In all cases, the output spikes followed the input spikes quite closely during the burst, then one or two extra spikes were generated and the pause occurred. However, at a closer inspection of raster plots, the precision of action potential emission changed. By increasing P (from P = 0.13 to P = 0.91) the PC model responses were characterized by a greater precision, with decreased variability in burst spike pattern and pause length (see also Figure 2C). Increased precision derived from the decreased paired pulse ratio (A2/A1, where A1 and A2 are the amplitudes of the first and second response in a pair). In this way, the first EPSP in the train became more precisely aligned with the input burst leading to a repeatable spike generation independent from the relative phase of background firing activity (see also Figure 8 below). Therefore, the highest precision was obtained when the pf-PC synapse expressed presynaptic pf-PC LTP. This suggests that release probability controls fine regulation of spike timing on millisecond scale.

In the case of a postsynaptic change in synaptic conductance G (−30% to 30% with respect to control), the only remarkable effect was a reduction in burst spikes precision when G was reduced, while no appreciable changes were observe in pause precision. The spike precision decrease at low G was likely related to a greater influence of previous spikes on burst initiation. Therefore, in the present conditions, precision could be more effectively tuned by P than G.

The PC shows complex branching (Nedelescu and Abdelhack, 2013) and active electroresponsiveness in the dendrites, so that the specific location of afferent synapses may influence the response pattern. The model was exploited to redirect to specific dendritic sectors (numbered I–IV) the same stimuli that were used before for synapses randomly distributed over the whole dendritic tree (see Figures 1–6).

The PC response was different depending on the stimulated sector (Figure 6A). When excitatory pf and inhibitory stimulations were delivered to dendritic sector I or IV, the PC BPRs were similar to those obtained using a random distribution of synapses. However, when excitatory stimulation was delivered to sector II or III, the pause was much longer than usual (over 250 ms in the Z+ model). When aa was substituted to pf stimulation, the responses were similar except for sector I, which showed a longer pause with aa then pf stimulation (the reason of this will be explained below). The difference between sectors was poorly sensitive to the location of inhibition, and BPRs did not change remarkably when inhibition was moved to sectors different from the one that was excited by pf synapses (Figure 6B). Similar response properties were observed in the Z− model, although pauses were shorter (data not shown).

In order to understand the mechanism differentiating responses among dendritic sectors, we compared sector I to sector II, which showed a remarkably different pause lengths. The impact of dendritic structure was considered first (Figure 7A). The local currents generated by synaptic activation (that were identical in the two sectors) caused a stronger depolarization in sector II than sector I. This reflected the different Z_in, that was about twice as large for sector II than sector I. Accordingly, this determined different voltage attenuation profiles, so that sector I depolarization started from a higher level and then decayed over a longer distance while approaching the soma. It should be noted that, in these simulations, we used single-synapse stimuli causing small depolarizations from a hyperpolarized membrane potential, so that voltage-dependent currents were not remarkably activated.

Then we considered the activation of voltage-dependent currents in the dendrites (Anwar et al., 2014) by delivering the appropriate multisynaptic stimulation pattern to the PC model in pacemaking regime. The higher synaptic depolarization in sector II than in sector I resulted in a stronger voltage-dependent activation of Ca channels in the former than in the latter (Figure 7B). In sector II, LVA currents amplified the EPSPs leading membrane potential to raise enough to activate also the HVA currents, thus causing a regenerative calcium spike. Consequently, a large raise in intracellular calcium activated the KCa system. This effect eventually influenced the spike-generating mechanisms in the AIS and regulated the pause, that was longer in sector II than sector I.

The non-linear nature of voltage-dependent mechanisms of BPR and their topographical nature bear about two main consequences.

A last question is how responses generated in different dendritic sectors can integrate to generate BPRs.

This investigation was made possible by the use of an advanced model of the PC (Masoli et al., 2015), which was developed from earlier ones (De Schutter and Bower, 1994a,b; Rapp et al., 1994; Roth and Häusser, 2001; Steuber et al., 2007; Masoli et al., 2015) by including novel electrophysiological features and was validated against a large set of recent experimental data. The model was implemented with ionotropic synapses coming from granule cells (both aa and parallel fibers) and SCs. Additional synaptic mechanisms that may affect the response to specific input patterns (e.g., see Barbour et al., 1994; Takahashi et al., 1995; Tabata et al., 2005; Blot and Barbour, 2014) remain to be assessed.

The pf synaptic transmission was calibrated to match responses evoked by random input patterns (Dittman et al., 2000), ensuring that pf synapses could precisely reproduce the temporal dynamics of short-term synaptic plasticity. The aa were made identical to pf synapses but were located on distal rather than proximal dendrites. A potentiation of the aa (Sims and Hartell, 2005, 2006; Walter et al., 2009) was required to counterbalance the longer electronic distance (L = 1.4 vs. L = 0.4) from soma (Roth and Häusser, 2001), thus ensuring the reported functional equivalence of aa and pf synapses (Walter et al., 2009).

It should be noted that in no case the model generated bistable switching between up-down states (Loewenstein et al., 2005), which may require additional mechanisms or specific modulation of receptors and ionic channels.

Model simulations brought about a number of testable predictions.

These observations suggest that BPRs represents a flexible coding strategy accounting for the timing, duration and intensity of input granule cell spike patterns and that can be modulated by the spatiotemporal distribution of the input, by its intensity, by synaptic plasticity and by inhibitory interneuron activity.

Model simulations also helped hypothesizing how PC dendrites process incoming synaptic inputs through BPRs. Dendritic Ca channels, primed by pf or aa EPSPs, would generate a depolarizing current flowing through the dendrites and activate a spike burst in the AIS. The corresponding raise in [Ca]_i would activate KCa channels causing a protracted hyperpolarization, interrupting pacemaking and generating the pause. The whole process turned out to be very sensitive to the location of synaptic inputs on the dendrites. Actually, EPSPs amplitude was larger at longer electrotonic distances and this enhanced the threshold crossing for generating local Ca spikes, that could then spread around invading neighboring dendritic regions. Thus, the combined effect of dendritic structure and voltage-dependent ionic channels generated non-linear interactions sculpting the spatio-temporal profile of local [Ca]_i and BPRs. It is tempting to speculate that this would eventually extend the plastic and computational properties of PCs beyond the linear perceptron hypothesized previously (Brunel et al., 2004).

A peculiarity of the PC dendrites is the enrichment in Ca channels (Llinás and Sugimori, 1980a,b), which include both HVA (Cav2.1) and multiple LVA (Cav3.1, Cav3.2, Cav3.3) subtypes (Ly et al., 2016). While Cav2.1 (the classical PC P-type channel) was known to regulate KCa1.1 (BK) channel activation and therefore the fast action potential AHP, the role of LVA channels remained uncertain. The present simulations show that LVA Ca channels (including at least cav3.1 and Cav3.3) are critical for BPR, since they are activated during the burst and then play a local but key role in activating the KCa2.2 (SK2) and KCa3.1 (SK4) channels, thereby regulating the pause. The relevant role of LVA calcium channels in EPSP generation by input bursts has recently been reported experimentally (Ly et al., 2016).

It should be noted that, when simulating intrinsic electroresponsiveness with somatic recordings, LVA Ca channels and KCa channels proved either noncritical or subcritical (Masoli et al., 2015). Therefore, in light of the present simulations, the tuning procedures of maximum ionic conductances in neuronal models should be targeted toward features addressing synaptic response pattern (Marasco et al., 2013) rather than just responses generated by somatic current injection (Druckmann et al., 2007, 2008; Masoli et al., 2015, 2017). Since some of these channels are located in the spines (specifically Cav3.1), a future assessment of the impact of spines in the PC model is warranted.