Firing rate-dependent phase responses of Purkinje cells support transient oscillations

Both spike rate and timing can transmit information in the brain. Phase response curves (PRCs) quantify how a neuron transforms input to output by spike timing. PRCs exhibit strong firing-rate adaptation, but its mechanism and relevance for network output are poorly understood. Using our Purkinje cell (PC) model, we demonstrate that the rate adaptation is caused by rate-dependent subthreshold membrane potentials efficiently regulating the activation of Na+ channels. Then, we use a realistic PC network model to examine how rate-dependent responses synchronize spikes in the scenario of reciprocal inhibition-caused high-frequency oscillations. The changes in PRC cause oscillations and spike correlations only at high firing rates. The causal role of the PRC is confirmed using a simpler coupled oscillator network model. This mechanism enables transient oscillations between fast-spiking neurons that thereby form PC assemblies. Our work demonstrates that rate adaptation of PRCs can spatio-temporally organize the PC input to cerebellar nuclei.


Introduction 21
The propensity of neurons to fire synchronously depends on the interaction between cellular 22 and network properties (Ermentrout et al., 2001). and shows significant adaptation with spike rate. It was theoretically predicted that PRCs 29 decrease at high firing rates in pyramidal neurons (Gutkin et al., 2005), which unfortunately 30 did not match later experimental observations showing an increased PRC peak at higher rates 31 (Tsubo et al., 2007). Similarly, for PRCs in Purkinje Cells (PCs), the responses to weak stimuli 32 at low spiking rates are small and surprisingly flat. With increased rates, responses in later 33 phases become phase-dependent, with onset-phases left-shifted and gradually increasing peak 34 amplitudes, which has never been theoretically replicated or explained (Couto et al., 2015;35 Phoka et al., 2010), nor has its effect in synchronizing spike outputs been explored. 36 On the circuit level, reciprocal inhibition causing high frequency oscillations has been 37 observed in many regions of the brain such as cerebellum and hippocampus (Bartos et  regulate the oscillations is still poorly understood. Furthermore, the functional importance of 40 oscillations in information transmission will be largely determined by their spatio-temporal 41 scale, which is difficult to predict given the hard-wired inhibitory connections. Since the PRC 42 is spike rate-dependent, it is interesting to explore whether this cellular property can 43 dynamically regulate the spatial range of oscillations based on spike rate changes. 44 To examine the mechanism of rate-dependent PRCs, we use our physiologically detailed PC 45 model (Zang et al., 2018) and a simple pyramidal neuron model to explore the rate adaptation 46 of PRCs. By analyzing simulation data and in vitro experimental data (Rancz and Hausser,47 2010), we identify that rate-dependent subthreshold membrane potentials can modulate the 48 activation of Na + channels to shape neuronal PRC profiles. We also build a PC network model 49 connected by inhibitory axon collaterals to simulate high-frequency oscillations (de Solages et 50 al., 2008; Witter et al., 2016). Rate adaptation of PRCs can increase the power of oscillations 51 to link the rate with spike timing. Moreover, firing irregularity and network connectivity also 52 regulate the oscillation level. The combination of these factors enables PC spikes uncorrelated 53 at low basal rates to become transiently correlated in a subgroup of cells at high rates. 54   To answer the aforementioned questions, we examine how spike properties vary with rate 90 and find that the facilitation of Na + currents relative to K + currents, due to elevated subthreshold 91 membrane potentials at high rates, underlies the rate adaptation of PRCs. After each spike, 92

PRC Exhibits Rate Adaptation in PCs
there is a pronounced after-hyperpolarization (AHP) caused by the large conductance Ca 2+ -93 activated K + current, and subsequently the membrane potential gradually depolarizes due to 94 intrinsic Na + currents and dendritic axial current (Zang et al., 2018). As confirmed by re-95 analyzing in vitro somatic membrane potential recordings (shared by Ede Rancz and Michael 96 Häusser (Rancz and Hausser, 2010)), subthreshold membrane potential levels are significantly 97 elevated at high rates, but spike thresholds rise only slightly with rate ( Fig. 2A). This means 98 the ISI phase where Na + activation threshold (~ -55 mV for 0.5% activation in PCs (Khaliq et 99 al., 2003;Zang et al., 2018)) is crossed shifts left and larger phase ranges of membrane 100 potentials are above the threshold at high rates (Fig. 2B). 101 During early phases of all rates, membrane potentials are distant from the Na + activation 102 threshold ( Fig. 2A,B). The depolarizations to weak stimuli fail to activate sufficient Na + 103 channels to speed up voltage trajectories, and phase advances are caused by just the passive 104 depolarizations (Fig. 2C). Consequently, phase advances in early phases are small and flat (or 105 phase independent). At later phases, membrane potentials gradually approach and surpass the 106 Na + activation threshold. Stimulus-evoked depolarizations activate more Na + channels to speed 107 up trajectories in return. Therefore, phase advances become large and phase-(actually voltage-) 108 dependent. At low rates, membrane potentials are below the Na + activation threshold during 109 nearly the entire interspike period (Fig. 2B). Responses are thus generally phase-independent. 110 At high rates, onset-phases of phase-dependent responses parallel the left shifts of Na + 111 activation threshold-corresponding phases, due to elevated subthreshold membrane potentials. 112 Because high rate-corresponding elevated membrane potentials have larger slopes at the foot 113 of the Na + activation curve, a same ∆V (passive depolarization) activates more Na + channels 114 and consequently causes larger PRC peaks at high rates (Fig. 2C,D). Under all conditions 115 (except phase = 0.2, 162 Hz), stimulus-evoked depolarizations also increase outward currents, 116 but this increase is smaller than that of inward currents (mainly Na + ) due to the high activation 117 threshold of K + currents (mainly Kv3) in PCs (Martina et al., 2003;Zang et al., 2018). As the 118 stimulus becomes stronger, its triggered passive depolarization increases and the required pre-119 stimulus membrane potential (phase) to reach Na + activation threshold is lowered (left shifted). 120 Thus, increasing the stimulus amplitude not only increases the PRC peaks, but also shifts the 121 onset-phases of phase-dependent responses to the left (Fig. S1C PRCs compared to PCs (Fig. 3A). In the Traub model, responses become smaller and relatively 140 phase-independent at high rates. In contrast to PCs, subthreshold membrane potentials are 141 significantly lower at high rates due to the accumulation of delayed rectifier K + current (kdr, 142 Fig. 3B,C), which has a low activation threshold and large conductance. The lower 143 subthreshold membrane potentials are far below the Na + -activation threshold, making 144 responses to weak stimuli passive at high rates. Accordingly, PRCs in the model become 145 smaller and relatively phase-independent at high rates, this was not confirmed in more recent 146 experimental recordings (Tsubo et al., 2007). We minimally modified the Traub model by 147 reducing the conductance of the kdr current, raising its activation threshold and increasing the 148 AHP current (details in Methods) ( high-frequency oscillations and each PC fires independently, as evidenced by the flat power 196 spectrum (Fig. 4D). High-frequency oscillations and their rate-dependent changes are also 197 reflected in the average normalized cross-correlograms (CCGs) between PC pairs (Fig. 4E). 198 When PCs fire at 70 Hz and 116 Hz, in addition to positive central peaks, two significant side 199 peaks can be observed in the CCGs, suggesting correlated spikes with 0 ms time lag and ~ 6 200 ms time lag. Amplitudes of the peaks increase with cellular rates and disappear when the rate 201 is low (10 Hz). Notice that the computed CCG shows 'excess' correlation, which is computed 202 by the raw correlation minus the shift predictor (Heck et al., 2007;Smith and Kohn, 2008). 203

204
In the previous section, the variation of cellular rates was driven by synaptic input to 205 demonstrate the rate-adaptation of high-frequency oscillations. However, it is difficult to 206 differentiate the relative contribution of PRC shapes and firing irregularity (measured by the 207 CV of ISIs) since they covary with rate (Fig. 4D). Therefore, cellular rates were systematically 208 varied by dynamic current injections, which were approximated by the Ornstein-Uhlenbeck 209 (OU) process (Destexhe et al., 2001). This simulation protocol also causes the formation of 210 high-frequency oscillations (Fig. S2). When PCs fire with low to moderate CV of ISIs, they 211 show loose spike-to-spike synchrony at high rates, but not at low rates. With high CV of ISIs, 212 spikes are jittered and spike-to-spike synchrony is disrupted. High-frequency oscillations were 213 never observed under the condition of low cellular rates. Both irregular spiking (high CV of ISIs) in A or low cellular rates in B decorrelate network 217 output in the forms of reduced peaks of power spectrums (left) and CCGs (right). In A, the 218 cellular rate is ~ 141 Hz. In B, the CV ISI is ~ 0.45. Both small conductance (cond) of inhibitory 219 synapses in C and a short connection radius (D) decorrelate network output in the forms of 220 reduced peaks of power spectrums (left) and CCGs (right). In A&B, the cond is 1 nS and radius 221 is 5. In C&D, the cellular rate is ~ 151 Hz and the CV ISI is ~ 0.45. 222 Both spiking irregularity and rates of PCs covary with cerebellum-associated behaviors 223 (Chen et al., 2016). Our results show that small spiking irregularity supports high-frequency 224 oscillations. However, further increasing spiking irregularity reduces the power of high-225 frequency oscillations and makes the power spectrum flatter when rates are the same (141 Hz) 226 (Fig. 5A). In average normalized CCGs, both central and side peaks decrease with increased 227 firing irregularity. Both results suggest reduced correlation when PCs fire very irregularly. 228 Next, we explore how rate-dependent PRCs regulate network output. The power of high-229 frequency oscillations increases and the power peak becomes sharper with large, broad PRCs 230 at high rates (Fig. 5B). Peaks in the average CCGs also increase, suggesting more correlated 231 spike outputs at high rates. In Fig. S3, we demonstrate that the PRC peak amplitude at high 232 cellular rates controls the oscillation power. 2016), the network generates robust high frequency oscillations (Fig. 5C, D). In addition, we 238 find that increasing the conductance of inhibitory synapses or their connection radius increases 239 the power of high-frequency oscillations and make the power spectrum sharper. The increased 240 oscillation power due to connectivity properties is also captured by the larger peaks in CCGs.

241
Both effects can be explained by larger phase responses due to larger inputs (synaptic 242 connections, Fig. S1C). Together, our simulation data suggest that the correlation between PC 243 spikes is strong under conditions of low to moderate spiking irregularity, high cellular rate, 244 high synaptic conductance, and large connection radius. 245 Transient Correlations Are Robust to Heterogeneous Spike rates 246 Although inhibitory connections loosely synchronize spike output and cause oscillations, 247 their functional role will depend on the following answers: When do they occur? How fast can 248 they be achieved? Are they robust to heterogeneous cellular spike rates? We have previously 249 simulated networks with a range of homogeneous stable cellular rates. Here, we first test 250 whether rate-dependent synchrony still holds when population rates change dynamically.

251
Population rates of the network approximate the half-positive cycle of a 1 Hz sine wave (peak 252 ~ 140 Hz) with the duration of each trial being 1 sec (Fig. 6A). We computed shuffle-corrected, 253 normalized joint peristimulus time histograms (JPSTHs) to reflect the dynamic synchrony 254 (Aertsen et al., 1989) (Fig. S4A). The main and the third diagonal of the JPSTH matrix, 255 corresponding to 0 ms time lag correlation and 6 ms time lag correlation respectively, are 256 plotted to show the dynamic synchrony at transiently increased rates (bin size is 2 ms, Fig. 6B rates of increased-rate cells and decreased-rate cells are shown in Fig. S4B). They were changes. From C to E, the number of decreased rate cells increases from 10 to 30. 281 larger PRCs peaks at high rates (Fig. 3) as in experiments (Tsubo et al., 2007). To the best of 320 our knowledge, decreased neuronal PRCs at high firing rates haven't been experimentally 321 observed yet. 322

323
Whether spike timing carries critical information for cerebellar function is still controversial.

324
For smooth pursuit eye movement in monkeys, the movement was reported to be coded just by increase the oscillation level at high rates, with no need to increase afferent input correlation. 339 Note that rate-related synchrony can also be achieved via common synaptic inputs (Heck et al., conditions, when the PRC is constantly 0 (equivalent to disconnection), no correlations can be 356 achieved (Figs. 4-6). Additionally, the oscillation power increased by the application of WIN 357 55,212-2, which was intended to suppress background excitatory and inhibitory synapses (de 358 Solages et al., 2008). The increased power could be due to more regular spiking after inhibiting 359 the activity of background synapses (Fig. 5A). However, it could also be caused by increased 360 spike rates (Fig. 5B), because this agent also blocks P/Q type Ca 2+ channels and consequently 361 P/Q type Ca 2+ -activated K + currents, to increase spike rates (Fisyunov et al., 2006). Similarly, 362 enhanced oscillations have also been observed in calcium-binding protein gene KO mice, 363 which accompany significantly higher simple spike rates (Cheron et al., 2004). 364

365
We built a link between oscillation level and firing irregularity (Fig. 5A) Sarnaik and Raman, 2018), optimal cerebellar function seems to occur between excessive 371 asynchrony and synchrony (Shakkottai, 2014). It has been shown that synchronized simple 372 spikes are time-locked to reaching-grasping movements in rats (Heck et al., 2007). Before and 373 after such movements, synchronized simple spikes were not observed. High-frequency 374 oscillations time-locked to lever-pressing in rats have been reported (Groth and Sahin, 2015). 375 Additionally, in humans high-gamma oscillations in the cerebellar cortex significantly increase 376 during a reaching task (Carver et al., 2019). Nonetheless, global and rhythmic increased 377 synchrony (Cheron et al., 2004) may abrogate normal separation of timing signals to different 378 muscle groups (for example agonist and antagonist muscles), causing impaired motor 379 coordination, such as dystonia (Shakkottai, 2014). In our model, PRCs are quantitatively close 380 to experimental data (Fig. 1). When cell and network parameters fall within physiological 381 ranges (Figs. 5,6), the network shows very weak oscillations at low basal cellular rates, but the 382 PC ensemble can dynamically increase the correlation level within a subgroup of PCs with 383 increased rates (Person and Raman, 2011) (Fig. 6). Local gamma oscillations have been shown 384 to selectively route input information in a cortical circuit model (Palmigiano et al., 2017 Exp2Syn. G = weight * (exp(-t/τ2) -exp(-t/τ1)), with τ1 = 0.5 ms (rise time) and τ2 = 3 ms (decay 447 time). The reversal potential of the IPSC was set at -85 mV (Watt et al., 2009 frequency oscillations, we varied the cellular rates in two paradigms. In the first (Fig. 4)  existing data support decreased PRC at high cellular rates. In our model, in a physiological 467 range, it is not available either. We therefore reduce the phase response by halving the input 468 current (synaptic conductance) to achieve a smaller response at high firing rates (Fig. S3). We 469 also explored the effect of connection radius with the values of 3, 5, 7 and10 in Fig. 5D. The 470 conductance of inhibitory synapses was tested with the values of 0.75, 1.0, 1.25 and 1.5 ns in 471 Fig. 5C. To test the spatio-temporally increased correlation, we randomly distributed extra 10 472 -30 PCs with decreased cellular rates into the original network (Fig. 6), including 100 473 increased-rate cells. Their mean population firing rates are shown in Fig. S4B. 474

Data analysis 475
The power spectrum of the spike trains of the network was estimated by Welch's method, 476 which calculates the average of the spectra of windowed segments (window size 128 points). 477 In each trial under each specific stimulus condition, the length of the signal was 2 sec, with a 478 time resolution of 1 msec. The final result was the average of 14 trials. 479 To compute the cross-correlogram (CCG) under each specific stimulus condition, we first 480 computed pairwise correlations between the spike trains of two neurons and then corrected 481 them by shift predictors, which removed the 'chance correlations' due to rate changes. Then 482 correlations were divided by the triangular function Θ( ) = − | | and D < F . T was the 483 duration of each trial and was the time lag. Θ( ) corrects for the degree of overlap between 484 two spike trains for each time lag . < was the mean firing rate of neuron i (Kohn and Smith, 485 2005). Finally, the CCGs between all pairs in the network were averaged to reflect the 486 population level spike correlations. Thus, similar with previous work (Heck et al., 2007), the 487 computed CCGs reflect the 'excess' correlation caused by axon collaterals in our work. To 488 measure the dynamic correlation over the time course of the stimulus, we computed JPSTHs. 489 Here we used a 2-ms time bins due to the narrow central peak and side peaks. Larger time bins 490 annihilated the positive peaks due to the significant negative correlations in paired spikes. 491 Therefore, we simulated 1992 trials to compute the raw JPSTH between PC pairs. Similar with 492 CCGs, the raw JPSTH was corrected by subtracting the shift predictor (cross-product matrix 493 of individual peri-event time histograms) to remove the coincident spikes due to rate changes 494 and co-stimulus. The corrected JPSTH was then normalized by the squared root of product of 495 each neuron's PSTH standard deviations (Aertsen et al., 1989). The corrected matrix values 496 become correlation coefficients, with values between -1 and +1. Finally, all pair-wise JPSTHs 497 were averaged to reflect the population level dynamic correlations. 498