Using subthreshold events to characterize the functional architecture of electrically coupled networks

The electrical connectivity in the inferior olive (IO) nucleus plays an important role in generating well-timed spiking activity. Here we combined electrophysiological and computational approaches to assess the functional organization of mice IO nucleus. Spontaneous fast and slow subthreshold events were commonly encountered during in vitro recordings. We show that the fast events represent a regenerative response in unique excitable spine-like structures in the axon hillock, whereas the slow events reflect the electrical connectivity between neurons (‘spikelets’). Recordings from cell pairs revealed the synchronized occurrence of distinct groups of spikelets; their rate and distribution enabled an accurate estimation of the number of connected cells and is suggestive of a clustered organization. This study thus provides a new perspective on the functional and structural organization of the olivary nucleus, insights into two different subthreshold non-synaptic events, and a novel experimental and theoretical approach to the study of electrically-coupled networks.


Introduction 41
In recent years research has confirmed that electrically coupled neural networks are found in 42 every major region of the central nervous system (Condorelli et  Marder, Gutierrez and Nusbaum, 2016; Coulon and Landisman, 2017). It thus stands to reason 51 that the functional architecture of these networks must undergo continuous modification to 52 meet system demands. This underscores the urgent need to determine the functional state of a 53 network and associate it with the corresponding brain states. Since anatomical information is 54 insufficient, this can only be done using physiological parameters that capture the functional 55 architecture of a network at any given time. 56 The inferior olive electrically coupled network, which was among the first networks to be studied 57 in the mammalian brain, provides primary excitatory input to the cerebellar cortex (Eccles,Llinás 58 and Sasaki, 1966). There is a general consensus that the function of this network is to generate 59 synchronous activity of the olivary neurons, which provide temporal information for either 60 learning processes, motor execution, sensory predictions or expectations (Llinas and Sasaki,61 1989; Lou and Bloedel, 1992  The relatively broad range of spikelet parameters (Figure 1) can be attributed to a wide range of 164 coupling strengths, different locations of the gap-junctions along the dendritic structure or 165 different durations and shapes of the pre-junctional action potential which is a well-known 166 feature of olivary action potentials (Llinás andYarom, 1981a, 1981b). We first examined the 167 effect of coupling strength by calculating the ratio of the amplitudes of the pre-junctional action 168 potential to the post-junctional spikelet, and compared it to the coupling coefficient measured 169 by direct current injection (see Methods). As shown in Figure 2C, there was a significant positive 170 correlation (with a slope of 0.134; R 2 =0.614, p<0.0001; Pearson correlation). Next, we examined 171 the effect of the shape of the pre-junctional action potential on the spikelet parameters. To that 172 end, we partially blocked the voltage dependent potassium current by adding TEA (10mM) to the 173 bath solution. In the presence of TEA, a variety of action potential waveforms were elicited by 174 current injection ( Figure 2D). In particular, the initial upstroke of the action potential was 175 unaffected, but there was a significant broadening of the repolarizing phase ( Figure 2D, upper 176 panel, dark blue) that often elicited a second calcium-dependent action potential ( Figure 2D, 177 upper panel, light blue). This variety of action potential waveforms were always associated with 178 postsynaptic responses that could be clustered into two distinct groups ( Figure 2D, lower panel). 179 The prolongation of the action potential was, as expected, followed by a matching increase in the Finally, we examined the occurrence of spikelets in neurons that exhibited subthreshold 186 oscillatory activity. Since these oscillations occurred simultaneously in several neurons (Lefler,187 Torben-Nielsen and Yarom, 2013), it was expected that spikelet occurrence will be correlated 188 with the oscillatory activity. About 50% of the olivary neurons showed spontaneous subthreshold 189 oscillations ( Figure 2E). Careful examination of the peaks of the oscillations ( Figure 2F) revealed 190 that they were crowned with spikelets. To quantify this observation, we calculated the 191 distribution of the inter-spikelet-interval (ISLI, Figure 2G, black bars), and found distinct groups 192 appearing at intervals of 200 ms. We then calculated the autocorrelation function of the 193 subthreshold oscillations ( Figure 2G, yellow line) and found that it matched the ISLI perfectly. It 194 is important to note that a similar fit was observed in 60% of the oscillating neurons (n=18) 195 whereas in non-oscillating neurons (n=70) the ISLI exhibited a Poisson-like distribution ( Figure  196 2H). The strong correlation between oscillatory behavior and the occurrence of spikelets further 197 supports the conclusion that these events represent activity in adjacent electrically coupled 198

neurons. 199
To summarize, the observations in Figures 1 and 2 strongly suggest that the small and slow 200 spontaneous events represent action potentials occurring in electrically coupled neurons and 201 therefore we will refer to them as 'spikelets'. 202 with relatively long durations (upper panel, blue traces) were elicited in cell 1 by 50 pA ,20 ms current pulse. Occasionally they were followed by a second response (cyan traces). These action potentials elicited post junctional responses in cell 2 with corresponding waveforms (lower traces). (E) Subthreshold events recorded in oscillating olivary neuron. (F) Superposition of the gray rectangles in E, at higher magnification. Note that spikelets were only present for 50 ms along the peak of the oscillations. (G) Inter-spikelet interval (ISLI) from the same neuron (using 4 ms bins), and an autocorrelation (yellow line) of the membrane potential. (H) The ISLI distribution in a non-oscillating neuron.

205
The large and fast events reflect localized regenerative responses. 206 The source of the fast events is a mystery. On the one hand, they seem to be independent of 207 chemical synaptic transmission, and on the other, they are not activated by pre-junctional action 208 potentials in coupled neurons. To resolve this mystery, we first examined the effect of the 209 membrane voltage on the occurrence and waveform of both types of unitary events. We next 210 used an optogenetic approach combined with pharmacological manipulations to gain insights 211 into the source of these events. Finally, we used a computational approach to examine possible 212 mechanisms to account for the experimental observations. 213 The effect of membrane potential was examined by DC current injection, which on average set 214 the membrane potential to a range of -33 to -90mV. Figure 3A-B shows the aligned superimposed 215 traces of spikelets ( Figure 3A) and fast events ( Figure 3B) from one neuron. Normalizing the event 216 amplitudes ( Figure 3A-B, right panels) shows that whereas the spikelet shape was unaffected by 217 the current injection (A), the fast events showed a slowdown of the late repolarizing phase with 218 hyperpolarization (B). Quantifying the effect of the injected current (see Methods) on the 219 amplitude ( Figure 3C-D) and duration at 20% of the amplitude ( Figure 3E-F) in 16 neurons 220 revealed no effect on the amplitude of either type of events (average slope R 2 = 0.0057 and 221 0.0035, respectively). The spikelet duration was slightly, but not significantly, affected ( Figure 3E; 222 R 2 = 0.44; one-sample t-test p=0.057). In contrast, it significantly increased the duration of the 223 fast events ( Figure 3F; R 2 = 0.712; one-sample t-test p=0.0005). Comparing the two sets of data 224 revealed a significant difference ( Figure 3E vs. F; paired t-test p=0.017). Finally, we measured the 225 effect of the DC current injection on the rate of occurrence of the subthreshold unitary events 226 ( Figure 3G-H). Whereas the frequency of the spikelets remained unaffected ( Figure 3G, R 2 = 227 0.012), the frequency of the fast events increased by a factor of up to 5 ( Figure 3H; R 2 = 0.836). 228 This difference between the occurrence of spikelets and the fast events was highly significant 229 ( Figure 3G vs. H, p=0.005, paired t-test). These differences further support our presumption that 230 two different mechanisms generate the two types of subthreshold unitary events. Since spikelets 231 are evoked by action potentials in paired cells, DC current injection to the post-junctional cell will 232 not have a significant effect on the number of action potentials in the pre-junctional cell. On the 233 other hand, the voltage dependency of the fast events strongly suggests that they reflect 234 regenerative responses initiated within the neuron. We therefore refer to them as Internal 235 Regenerative Events or IREs. 236 Next, we examined the ability of the IREs to trigger a full-blown action potential. Figure 4A shows 238 the superposition of spontaneously occurring action potentials (black and purple traces) and IREs 239 (gray traces) measured from the same neuron. It is clear that some of the action potentials 240 (purple) are preceded by a fast pre-potential that resembles the spontaneous IREs. This 241 possibility was further supported by phase plotting the action potential ( Figure 4B, purple and 242 black traces), which showed that 62% of the spontaneous action potentials were triggered by 243 depolarizing events that resembled the IREs in amplitude and waveform. On the population level 244 29.2±26.7% of the IREs triggered action potentials, whereas 41.6±21.2% of the action potentials 245 were triggered by IREs ( Figure 4C). It should be noted that this is an underestimation, since action 246 potentials that seem to arise from the baseline might also be evoked by an IRE that is masked by 247 the beginning of the action potential (see inset in Figure 4A). 248

249
We then examined the possibility of activating these events by synaptic inputs, using Thy1-ChR2-250 YFP transgenic mice in which ChR2 is expressed in axons that deliver excitatory inputs to the IO 251 but not in IO neurons (Lefler, Torben-Nielsen and Yarom, 2013). As shown in Figure 4D

264
These results strongly support the possibility that the fast events are generated within the neuron 265 either spontaneously or in response to synaptic input. Their fast kinetic suggests that sodium 266 current is involved; in fact, these fast events were not detected in 27 neurons that were injected 267 with QX-314 (2-5 mM). This contrasts with normal conditions in which only 28% of the cells were 268 devoid of fast events. Furthermore, in a complex experiment with two patch electrodes in the 269 same neuron, one filled with regular intracellular solution and the other with QX (Supplemental 270 information, Figure S1), we found that in the control condition both IREs and spikelets were 271 readily detected ( Figure  characteristic feature of the IREs is that for any given cell, their amplitude clearly segregates into 283 distinct groups, whereas their waveform is essentially identical ( Figure 1C-D). We used a 284 modeling approach to explore the conditions in which this feature could be reproduced, either 285 for dendritic 'hot spot' or for axonal spike failure. 286 To examine the possibility of "dendritic hotspots" we developed detailed compartmental models 287 of three 3D reconstructed olivary neurons; a sample modeled cell is shown in Figure 5A-C. To 288 simulate "hot spots", we injected, at every dendritic location, a simulated spike (Stuart and 289 Sakmann, 1994) with a 100 mV amplitude and a 2 ms half-width ( Figure 5D). We searched for the 290 dendritic locations where the resultant somatic voltage response resembled the experimental 291 IREs, with half-widths ranging from 3.25 to 3.5 ms (see Figure 1F). This yielded nine possible 292 confined dendritic locations in the modeled cell (colored dendritic segments in Figure 5C). The 293 respective simulated responses are shown in Figure 5E. Thus, although the amplitudes of these 294 responses varied from 2 mV to 6.5 mV ( Figure 5E), their waveforms were almost identical ( Figure  295 5F). Similar results were found in two other modeled cells, thus demonstrating that only well-296 placed dendritic "hot spots" could replicate the experimental results. 297 To examine the possibility of axonal spike failure we constructed a model of a myelinated axon 298 connected via an initial segment to an isopotential soma connected to several passive dendrites 299 ( Figure 5G; see Methods). We examined the somatic voltage response to an antidromic action 300 potential (activated by short current injection at node 10) during successive blockade of the 301 Nodes of Ranvier (starting with node 1 only, then node 1+2, etc.). As more nodes are blocked, a 302 smaller somatic response is expected ( Figure 5H, brown, red, yellow and green traces, 303 respectively). Surprisingly, normalizing the amplitude of these responses ( Figure 5I) revealed that 304 the waveform was almost completely unaffected by the location of the propagation block. This 305 was due to the small effective capacitance of the myelin. Indeed, increasing the myelin 306 capacitance enhanced the difference between the waveforms of the somatic responses ( Figure  307 5J-K). 308 Both models reproduced the characteristic features of the IREs; i.e., identical waveforms and 310 variable amplitudes. However, it is difficult to envisage a biological mechanism that either 311 specifically localizes channels in a restricted dendritic "hot spot" or that simultaneously blocks hillock of olivary neurons. We modeled these axonal spines while taking their different neck 316 lengths into account ( Figure 6A) and assumed that their head membrane is excitable (see 317 Methods). Figure 6B shows the somatic responses, ranging from 6 to 9 mV, that resulted from 318 action potentials generated in these axonal spines by a short current injection. The somatic 319 responses had identical waveforms ( Figure 6C). Importantly, these axonal spines are innervated depolarizing currents into the model soma altered the action potential waveform in the spine 325 head ( Figure 3D), particularly the after-hyperpolarizing phase, which, as expected, turned into a 326 depolarizing phase (black arrow in Figure 6E, left panel). As a result, the duration of the somatic 327 response increased (red arrow in Figure 6E, right panel). This effect was observed for three levels 328 of somatic current injections (-150 pA, -30 pA and +30 pA) for all four spines ( Figure 6F). The 329 normalized responses are displayed in Figure 6G, and illustrate the prolongation of the 330 repolarizing phase of the somatic responses with hyperpolarization, as found experimentally. 331

333
Estimating network architecture from dual cell recordings of simultaneously occurring spikelets. 334 Figure 7A depicts the spontaneous activity recorded simultaneously from two neurons. As 335 described above (Figure 2) action potentials (diagonal bars) that occurred irregularly in either of 336 the two neurons were always associated with spikelets in the paired neuron ( Figure 7B). The 337 subthreshold activity, which was dominated by spikelets as well as IREs, appeared randomly in 338 the two neurons. However, occasionally spikelets occurred simultaneously in both cells (marked 339 in Figure 7A and shown at high resolution in Figure 7C) which we refer to as 'common spikelets'. 340 Figure 7C, which occurred without measurable time 341 difference, have variable amplitudes. The first and the third spikelets had larger amplitudes in 342 the red neuron (cell 2) whereas the middle spikelet had a larger amplitude in the black neuron 343 (cell 1). Since action potentials in one neuron evoke very similar spikelets in the other (Figure 2A-344 B), the most likely explanation is that each of these common spikelets represents the action 345 potential in an additional neuron that is coupled to both of the recorded neurons (see 346 Discussion). On the population level, 18 out of 38 pairs had common spikelets (47%). Of these 347 pairs, the occurrence of common spikelets varied from 0.02 to 1.1 /sec, which is 3.5-66% of the 348 total number of measurable spikelets ( Figure 7D). 349

Each of the three examples shown in
The occurrence of common spikelets can be used to estimate the number of neurons that are 350 electrically coupled to each neuron in the olivary network. In this example of a paired recording, 351 4 different groups were identified by clustering the detected common spikelets by the amplitude 352 value and amplitude ratio between the two neurons ( Figure 7E Further insights into the organization of the network can be extracted from the distribution of 370 the number of groups of common spikelets. As shown in Figure 7G, the number of common 371 groups varied from 0 to 7 with a likely higher incidence at 2 -4 groups. As will be shown below, 372 this type of distribution cannot be the result of random connectivity where the probability of 373 connection depends solely on the distance between the recorded cells.

376
We re-examined the approach we used to estimate the number of connections per neuron by 377 reconstructing a realistic olivary network ( Figure 8A, see Methods). The firing rate of neurons in 378 the network was set to 0.058 Hz ± 0.04 Hz (as observed experimentally) and the number of 379 common spikelets in pairs of neurons occurring within 15 min of simulation was measured. 380 Recordings from a sample pair are shown Figure 8B. In this example, four groups of spikelets that 381 appear 26, 16, 65 and 32 times were detected ( Figure 8B). By applying the same calculation as 382 performed in the experimental observations (Figure 7), we concluded that the red neuron was 383 electrically connected to 24 neurons whereas the black neuron was connected 26 neurons. In 384 this model, the red and black neurons were actually connected to 17 and 20 neurons, respectively 385 (pink circles in Figure 8C). We performed the same calculation in 20 randomly selected pairs of 386 neurons and plotted the estimated versus the real number of connections per cell ( Figure 8C). 387 The results were distributed along the diagonal, (with the average marked by + sign), thus 388 demonstrating the validity of this approach in estimating the number of connections for each 389 neuron. However, the accuracy of the estimation depends strongly on the variability in firing rate 390 and the recording duration; the accuracy decreases the shorter the recording duration and the 391 higher the firing rate variability, but the error in estimating the average number of connections 392 for each cell was still small (Supplemental information, Figure S2, see discussion). 393 The experimental distribution of the number of groups of common spikelets ( Figure 7G) provides 394 additional information on the organization of the network. We examined two possible 395 connectivity patterns that might support this type of distribution (see Methods; network 396 connectivity matrices). The first is a network where the probability of a connection between two 397 neurons depend solely on their inter-somatic distance ( Figure 8D). The second assumes that the 398 network is organized into clusters of neurons ( Figure 8G), where the connection probability 399 within a cluster is larger than between clusters (both probabilities are distance-dependent). We 400 searched for the parameters that yielded the best fit to the experimental observations; namely, 401 the distance-dependent connection probability (see Devor and Yarom, 2002a and Figure 8E, 402 H) and the distribution of common neighbor (proxy for groups of common spikelets in Figure 8F, 403 I). Our search yielded an interesting result. Whereas the simple distance-dependent network 404 captured the distance-dependent probability of a connection ( Figure 8E), it failed to reproduce 405 the distribution of common groups as found experimentally ( Figure 8F). On the other hand, when 406 the modeled networks were organized in clusters, it replicated both the experimental 407 distribution of common groups ( Figure 8I) as well as the distance-dependent connection 408 probability ( Figure 8H). Note that in all modeled networks, each neuron was connected to about 409 11-21 neurons, which lay within the numbers estimated from the experimental observations 410 ( Figure 7F). 411 information, Figure S3). However, out of over 45 paired recordings we never encountered a pair 446 where an action potential in one neuron elicited an IRE in the coupled neuron. Furthermore, an 447 IRE in one neuron never coincided with a spikelet in the other neuron. Moreover, whereas the 448 probability of spikelets was independent of the membrane voltage, the IRE showed an increased 449 probability at depolarization (Figure 3). Finally, IREs were recorded from the olivary neuron in CX to account for the failure of axonal propagation, one has to assume that action potentials can be 474 generated down the axon away from the cell body (Pinault, 1995;Sasaki, 2013). Moreover, the 475 fact that the axon, or more likely the axon initial segment, is usually regarded as the site with the 476 lowest threshold (Palmer and Stuart, 2006), cannot be readily applied in our case, since 477 simultaneous failure at multiple nodes are needed to generate different amplitudes. 478 Finally, we were able to simulate the experimental observation assuming that the source of the 479 IREs are the spines that emerge from the axon initial segment. If these spines can support action 480 potentials, as suggested in other systems (Segev and Rall, 1998 where the specific activation of individual spines are bound to occur. Thus, it is more likely that 486 the source of the IREs are these axonal spines. 487 Spines emerging from an axon is not a common feature of central neurons. As we demonstrated, 488 IREs are strong enough to activate the neurons (Figure 4) and thus generate a unique situation 489 where the activation of the neuron is independent of the background synaptic activity that takes 490 place in the soma and the dendritic tree. It is tempting to speculate that this arrangement ensures 491 activation of the neurons in a state of "emergencies". Thus, when the system is in need of a fast 492 response, it activates these inputs that bypass the ongoing activity, thus insuring an immediate 493 response. 494

The electrically coupled network in the olivary nucleus 495
Spikelets and subthreshold oscillations: The formation of electrically coupled network within the 496 olivary nucleus is well-established. Here we suggested that this network, operating as a 497 synchronous rhythmic device, is capable of generating precise temporal patterns (Jacobson,498 Rokni and Yarom, 2008). Both synchronicity and rhythmicity are generated by the delicate 499 interplay between electrical coupling and ionic conductances. Thus, a single cell by itself cannot 500 oscillate, whereas in a network formation the neurons generate subthreshold oscillations (Manor 501 et al., 1997;Loewenstein, Yarom and Sompolinsky, 2001;Chorev, Yarom and Lampl, 2007). In 502 this work, we studied the relationship between spikelets and subthreshold oscillatory activity and 503 found that in oscillating cells the occurrence of spikelets coincided with the depolarizing phase 504 of the oscillation whereas in non-oscillating cells they seemed to be randomly distributed (Figure  505 2). This result strongly supports our previous hypothesis that the occurrence of subthreshold 506 oscillations is a network phenomenon. Therefore, when the recording cell is oscillating, the entire 507 network is synchronously oscillating, thus generating action potentials at the peak of the 508 oscillation that appear in the recorded cell as spikelets. Theoretically, by calculating the number 509 of spikelets at the peak of the oscillatory activity, one should be able to calculate the number of 510 coupled neurons in the oscillating network. Although tempting, this is practically impossible 511 because spikelets, given their small amplitude and noisy environment, cannot be classified into 512 groups. Therefore, we used the common spikelets to estimate the number of coupled neurons. 513 Common spikelets: Recording from two olivary neurons revealed spikelets that occurred 514 simultaneously in both recorded neurons. Given that the average rate of spikelets is 0.56±0.62Hz, 515 the probability that these common spikelets reflect random occurrence is extremely low. 516 Furthermore, the repeated appearance of common spikelets with the same relationships 517 (amplitude ratio and waveforms) further supports the non-random occurrence of these events. 518 Thus, the accurate timing of common spikelets can only be attributed to a common source; i.e., 519 a single pre-junctional neuron. The number of neurons that were coupled to the two recorded 520 neurons, which varied from 1 to 7, should be correlated with the size of the network; more 521 common spikelets are expected in larger interconnected networks. 522 Our simple method of calculating the size of a coupled network is based on data obtained during 523 simultaneous recordings from two neurons and on the assumption that the neurons display a 524 similar rate of spontaneous spiking activity. Our simulations show that the accuracy of this 525 method is mainly affected by variability in the firing rate (Supplemental information, Figure S2). 526 To minimize the error in estimating the firing rate of the neurons in the network, we used only 527 the spontaneous rate of the common spikelets, and not the firing rates in the recorded neurons In addition to calculating the number of connected cells, the distribution of common spikelets 536 enabled us to study the connectivity profile within the nucleus. Our data showed that each two 537 neurons had 1-7 common groups, and a normal distribution with an average of 3-4 groups. 538 However, there was a large number of paired recordings that failed to show common groups. 539 Using a theoretical approach, we demonstrated that this distribution should not be expected if 540 we assume that the probability of connection depends solely on the distance between the 541 neurons. On the other hand, if the nucleus is organized into clusters where the probability of 542 connection within the cluster is higher than between clusters, the observed distribution of 543 common groups can be reproduced. Although the size of the clusters, as well as the probability 544 of connection cannot be defined with the current data, this constitutes the first physiological 545 study that supports the assumption of clustered organization of the nucleus deduced mainly from 546 morphological studies. 547 In summary, we presented a comprehensive study that implemented a wide range of research 548 approaches to unravel the functional architecture of the inferior olivary network. We showed 549 that by analyzing spontaneous subthreshold events, new insights into the organization of the 550 network can be gained, thus paving the way for a novel experimental and theoretical approach 551 to the study of electrically-coupled networks. All experimental procedures were approved by the Hebrew University's Animal Care and Use 571 Committee. Brain stem slices were prepared from the following strains of mice: C57BL/6, B6.Cg-572 Tg (Thy1-COP4/EYFP)(Jackson Laboratory) and Gad2-tm2(cre)Zjh/J (Jackson Laboratory). 573

Electrophysiological recordings and ChR stimulation 581
The recording chamber was perfused with 95% O 2 and 5% CO 2 physiological solution at 24˚c-582 28˚c. Slices were visualized using a 40X water-immersion objective in an Olympus BX61WIF 583 microscope equipped with infrared differential interference contrast (DIC). In order to record 584 from intact olivary networks, recordings were targeted to the deepest neurons possible in the 585 slice. Whole-cell recordings were performed using 6-9 MΩ glass pipettes with intracellular 586 solution containing 4mM NaCl, 10 -3 mM CaCl 2 , 140mM K-gluconate, 10 -2 mM EGTA, 4 mM Mg-587 ATP, and 10 mM HEPES (pH 7.2). Signals were acquired at 10-20KHz using a Multiclamp 700B 588

Reconstructing neurons 595
In a few experiments (using C57BL/6 mice), Neurobiotin (0.5%; Vector Laboratories) was added 596 to the pipette solution to label the recorded neurons. Slices were then fixed in 4% 597 paraformaldehyde overnight, washed in PBS and stained with 1 µg/ml Streptavidin AlexaFluor 598

Data analysis and Statistics 600
Analysis was performed using MATLAB (R2014b and R2016a, MathWorks) for the experimental 601 data and Python 2.7 for the simulation data. The 71 neurons that were selected for detailed 602 analysis had a frequency of subthreshold events exceeding 0.02 Hz. The events were divided into 603 two different groups according to their amplitude and rise time. The event rise time was 604 calculated as 10%-90% of the amplitude. The IRE groups were clustered using the K-means 605 clustering method. The coupling coefficient (CC, Figure 2) was calculated as the ratio between 606 the change in the steady-state voltage of the post-junctional cell and that of the pre-junctional 607 cell in response to 250 ms current injection in the pre-junctional cell. The spike coupling 608 coefficient was measured as the amplitude of the post-junctional spikelet divided by the 609 amplitude of the pre-junctional action potential. A Pearson correlation was used to calculate the 610 p-value of the linear regression in Figure 2C. The Inter-spikelet-interval (ISLI) in oscillating 611 neurons ( Figure 2G) was only calculated in neurons that had more than 150 spikelets during the 612 session. The ISLI histogram was computed using 4 ms time bins, and the autocorrelation in 613 oscillating neurons was calculated using a lag of 1 ms. 614 The effect of the DC current injection in (Figure 3 C-H)  Common spikelets (Figure 7) were defined as spikelets that were detected in a paired recording 630 in both cells simultaneously. Groups of common spikelets were clustered using k-means analysis 631 on the peak voltage, followed by manual curation in which we verified that all common spikelets 632 in a group had a similar shape. To estimate the number of connections per cell, the total number 633 of spikelets (T spikelets ) was multiplied by the number of groups (N groups ) and divided by the number

636
If the two recorded cells were coupled (i.e., a spike in one cell gave rise to a spikelet in the other 637 cell), +1 was added to the estimation of connections for these two cells. 638

Neuron models 639
Using the Neurolucida software (MBF Bioscience), three olivary neurons were reconstructed 640 from fluorescence image stacks acquired using a Leica TCS SP5 confocal microscope (Leica 641 Microsystems). To compensate for tissue shrinkage, the z-axis of the reconstruction was 642 multiplied by a factor of 3. A compartmental model was generated from the morphological 643 reconstruction using NEURON (Carnevale and Hines, 2006). The axial resistance (R a ) was set to 644 100 Ωcm, the specific membrane capacitance (C m ) to 1 µF/cm 2 and the specific membrane 645 resistivity (R m ) for the three reconstructed cells were 4300, 4500, 3800 Ωcm 2 respectively. These 646 values were chosen to yield an input resistance (R in ) that was within the experimental range 647 (115±43 MΩ). 648

Axon model 649
The model was based on a previous axonal model (Chen et al., 2002) (ModelDB, accession 650 number 3793) with some modification. The model had the following compartments: soma, five 651 cylindrical dendrites, an axon hillock (AH), an axon initial segment (AIS), ten myelinated segments 652 and ten Nodes of Ranvier (Node). The dimensions and density of the channel conductances for 653 these compartments are summarized in Table 1 and are partly based on previous EM studies (De 654 Zeeuw et al., 1990). Blocking of the different nodes was implemented by removing all active 655 conductances at that node ( Figure 5). The axial resistance was 150 Ωcm; the specific membrane 656 capacitance (C m ) was 1 µF/cm 2 in all compartments except for the myelinated segments, where 657 the C m was 0.01 µF/cm 2 . The specific membrane resistivity (R m ) was 12,000 Ωcm 2 in all 658 compartments except for the myelinated segments, where R m was 1,200,000 Ωcm 2 so that the 659 myelin membrane time constant was 12 ms. When the myelin C m was modified, R m was set to a 660 value that kept the membrane time constant at 12 ms. 661

Model of axon with spines 664
Four spines were added to the axon model in the AH, 11.7, 15.3, 4.5 and 8.1 µm from the soma. 665 Each spine had a neck length of 8.5, 4.5, 7, 3.5 µm, respectively, and a diameter of 0.2 µm. The 666 spine head diameter and length were 2 µm each. The sodium and potassium conductances (gNa 667 and gK, respectively) in the AH and in the spine heads were set as in the AIS (Table 1). For these 668 simulations, the active conductances were removed from the Nodes of Ranvier. 669

Building the IO network connectivity matrices 670
We constructed a network of IO composed of 1134 neurons randomly distributed within a 671 volume of 250x500x200 µm, which resulted in 0.045 neurons per 10 µm 3 . We then clustered the 672 neurons by their location using k-means clustering, and varied the number of neurons in a cluster 673 by choosing k to be 1134 divided by the number of neurons in a cluster. The probability of a 674 connection between two neurons decays with distance according to a Gaussian profile: 675 Σ * 3 4 5 6 * 7 5 100 676 where S is the maximal probability for connection (when the distance between the neurons is 0, 677 x is the distance between neurons and σ sets the decay of connection probability with distance 678 (see Figure 8H inset for examples). Note that the shape profile of neuron connectivity within a 679 cluster could have a different S and σ than the connectivity profile of neurons belonging to 680 different clusters. The common neighbor distribution ( Figure 8F, I) was calculated on randomly 681 selected pairs of neurons within a distance of 40 µm. 682

Constructing the IO network model 683
To simulate a realistic network of IO neurons, we followed the steps described above but with a 684 few modifications. The network volume was 125x250x100 µm, and populated with 180 neurons 685 (0.057 neurons for 10 µm 3 ). These neurons were cloned from the three 3D-reconstructed olivary 686 neurons, one of which is shown in Figure 5. We set the cluster size, k, to 40. S and σ were 77 and 687 42 15 within a cluster and 45 and 20 between clusters, respectively ( Figure 8H, inset, gray). The 688 electrical connection between two neurons was mediated by 2 gap junctions (GJs). AGJ 689 conductance (GJc) of 0.3 nS resulted in a coupling coefficient of 0.03±0.019 as in the experimental 690 range (0.039±0.029). After adding GJc to the modeled cell, R m was modified to maintain the 691 experimental value of R in (See details in (Amsalem et al., 2016)). The spikes in the networks were 692 created by current injection to the soma (simulated spikes) following a Poisson process. 693

Supplementary Figures
Figure S1   Figure 8C) as a function of simulation duration for 6 different firing rate variabilities (std). The calculation was done only on neurons that had common neighbors (n=278 pairs). Mean firing rate was 0.058 Hz. (C) average normalized spikelets from B. This network is similar to the network in Figure 8A, with two modification -the number of GJs was changed to 6 and the GJc to 0.25 nS.