Interspike intervals within retinal spike bursts combinatorially encode multiple stimulus features

Neurons in various regions of the brain generate spike bursts. While the number of spikes within a burst has been shown to carry information, information coding by interspike intervals (ISIs) is less well understood. In particular, a burst with k spikes has k−1 intraburst ISIs, and these k−1 ISIs could theoretically encode k−1 independent values. In this study, we demonstrate that such combinatorial coding occurs for retinal bursts. By recording ganglion cell spikes from isolated salamander retinae, we found that intraburst ISIs encode oscillatory light sequences that are much faster than the light intensity modulation encoded by the number of spikes. When a burst has three spikes, the two intraburst ISIs combinatorially encode the amplitude and phase of the oscillatory sequence. Analysis of trial-to-trial variability suggested that intraburst ISIs are regulated by two independent mechanisms responding to orthogonal oscillatory components, one of which is common to bursts with a different number of spikes. Therefore, the retina encodes multiple stimulus features by exploiting all degrees of freedom of burst spike patterns, i.e., the spike number and multiple intraburst ISIs.


Unfunded studies
Enter: The author(s) received no specific funding for this work. The authors have declared that no competing interests exist.
If the data are held or will be held in a public repository, include URLs, accession numbers or DOIs. If this information will only be available after acceptance, indicate this by ticking the box below. For example: All XXX files are available from the XXX database (accession number(s) XXX, XXX.).
• If the data are all contained within the manuscript and/or Supporting Information files, enter the following: All relevant data are within the manuscript and its Supporting Information files.
• If neither of these applies but you are able to provide details of access elsewhere, with or without limitations, please do so. For example: Data cannot be shared publicly because of [XXX]. Data are available from the XXX Institutional Data Access / Ethics Committee (contact via XXX) for researchers who meet the criteria for access to confidential data.
The data underlying the results presented in the study are available from (include the name of the third party and contact information or URL).
• This text is appropriate if the data are owned by a third party and authors do not have permission to share the data.

• * typeset
Additional data availability information: Neurons in various regions of the brain generate spike bursts. While the number of spikes within 14 a burst has been shown to carry information, information coding by interspike intervals (ISIs) is 15 less well understood. In particular, a burst with k spikes has k−1 intraburst ISIs, and these k−1 16 ISIs could theoretically encode k−1 independent values. In this study, we demonstrate that such 17 combinatorial coding occurs for retinal bursts. By recording ganglion cell spikes from isolated 18 salamander retinae, we found that intraburst ISIs encode oscillatory light sequences that are 19 much faster than the light intensity modulation encoded by the number of spikes. When a burst 20 has three spikes, the two intraburst ISIs combinatorially encode the amplitude and phase of the 21 oscillatory sequence. Analysis of trial-to-trial variability suggested that intraburst ISIs are 22 regulated by two independent mechanisms responding to orthogonal oscillatory components, one 23 of which is common to bursts with different number of spikes. Therefore, the retina encodes 24 multiple stimulus features by exploiting all degrees of freedom of burst spike patterns, i.e., the 25 spike number and multiple intraburst ISIs. 26

27
Neurons in various regions of the brain generate spike bursts. Bursts are typically composed of a 28 few spikes generated within dozens of milliseconds, and individual bursts are separated by much 29 longer periods of silence (~hundreds of milliseconds). Recent evidence indicates that the number 30 5 ganglion cells, generated spike bursts (Fig 1A). The majority of the spikes [82.0% ± 8.7%, mean 79 ± SD (standard deviation), n = 41 cells] were observed in bursts comprised of two or more 80 spikes, indicating that multi-spike bursts represent the major retinal code. indicate intervals shorter and longer than Tthresh, respectively. 6 During repeated presentation of the same stimulus, individual ganglion cells generated spike 92 bursts at similar time-points across repeats ( Fig 1A) [3,4]. This reproducibility enabled the 93 identification of corresponding bursts across repeats, which we termed "events" (Fig 1 and  94 Materials and Methods) [3,4]. Bursts generated in the same event had similar numbers of spikes 95 in different repeats of the stimulus, while those in different events often had different numbers of 96 spikes ( Fig 1A). Accordingly, the number of spikes within bursts carried information about the 97 stimulus (p < 0.01 for 41 of the 41 cells; the estimated mutual information was 0.80 ± 0.31 bits 98 per burst, mean ± SD, n = 41). We next calculated the burst-triggered averages (BTAs), which 99 represent the average stimulus sequence preceding isolated spikes and bursts with two, three, and 100 four spikes (1-, 2-, 3-, and 4-BTA, respectively). The BTAs were sequences of different 101 amplitudes (Fig 2A), and the difference between 1-BTA and 3-BTA had ON and OFF peaks 102 around −170 and −40 ms relative to the burst onset ( Fig 2B). This result indicates that the 103 number of spikes within a burst encodes the amplitude of ON-to-OFF light intensity modulation 104 within an interval of ~130 ms.

Burst ISIs encode oscillatory light intensity sequences 112
To investigate whether intraburst ISIs carry information, we first analyzed bursts composed of 113 two spikes (2-spike bursts). For each ganglion cell, 2-spike bursts in the same event tended to 114 have similar ISIs, whereas those in different events typically had different ISIs (Fig 3A). This  Two-spike burst ISIs were not correlated to the average number of spikes in events ( Fig 3B  134 and C; correlation = 0.0 ± 0.1, mean ± SD, n = 41), suggesting that these ISIs were modulated 135 according to stimulus features different from those modulating the spike number. To characterize 136 the stimulus features modulating 2-spike burst ISIs, we extracted the stimulus sequences 137 preceding 2-spike bursts with different ISIs (Fig 3D and E). The results show that the stimulus 138 sequences had systematic differences depending on the ISIs. We next subtracted the 2-BTA 139 (black in Fig 3E) from the average of sequences preceding 2-spike bursts with long ISIs. The 140 result was an oscillating sequence with two ON and two OFF peaks (yellow in Fig 3F and G). 141 The subtraction from the 2-spike bursts with short ISIs gave the same sequence, but with the 142 opposite sign (blue in Fig 3F and G). These results indicate that 2-spike burst ISIs encode an 143 oscillatory deviation from the 2-BTA. Intervals between the ON and OFF peaks were ~40 ms 144 (Fig 3F and G) and, therefore, were much shorter than those of the stimulus feature encoded by 145 the spike number (~130 ms; Fig 2B). Thus, two-spike burst ISIs encode oscillatory sequences 146 much faster than the intensity modulation encoded by the spike number. 9 The two ISIs of three-spike bursts encode the phase and amplitude of oscillatory 148 components 149 We next investigated the characteristics of 3-spike bursts. The first and second intraburst ISIs 150 (ISI1 and ISI2) tended to be different for different events (Fig 4A and B) and carried information 151 about the stimulus (p < 0.01 for 40 of the 41 cells; 0.46 ± 0.21 bits per burst, n = 41). In addition, 152 the data suggest that ISI1 and ISI2 were modulated differently. For example, in Fig 4A, ISI1 was 153 similar between the first (red) and second (green) events, but tended to be different for the third 154 event (blue). In contrast, ISI2 was similar between the second and third events, but different for 155 the first event. Accordingly, in the two-dimensional plot of ISI1 and ISI2, bursts of different 156 events occupied different locations, and events did not align one-dimensionally, but were 157 distributed two-dimensionally ( Fig 4B). 158 The above results suggest that the modulation of ISI1 and ISI2 has two degrees of freedom. 173 The distribution of ISI1 and ISI2 had features that complicate further analysis. First, events with 174 longer ISIs tended to have larger trial-to-trial ISI variability than events with shorter ISIs (Fig  175   4C). This inhomogeneous variability suggests that shorter ISIs represent information with a 176 resolution higher than that represented by longer ISIs. Second, although ISI1 and ISI2 were 177 modulated differently, they were not completely independent, but were correlated (Fig 4B;  178 correlation coefficient = 0.16 ± 0.14, mean ± SD, n = 41). We were able to correct the first point 179 by non-linearly scaling ISI1 and ISI2 so that different events had similar variability (Fig 4C and  180 D). To correct the second point, variables were further linearly scaled to have an approximately 181 circularly symmetric distribution so that the correlation was negligible (Fig 4E; correlation 182 coefficient = −0.01 ± 0.13, n = 41). The resultant variables, " and # , had a circularly 183 symmetric distribution and, therefore, were approximately independent, with different events 184 occupying a similar amount of area ( Fig 4E). 185 Using " and # , we defined the "burst phase" for each burst (Fig 4B and E). Stimulus 186 sequences preceding 3-spike bursts exhibited systematic differences according to the burst phase 187 ( Fig 5A). We subtracted the 3-BTA (Fig 5B) from the sequences preceding 3-spike bursts. The 188 resulting deviations were oscillatory components with two or three ON peaks separated by 70-80 189 ms, with OFF peaks among them (Fig 5C-E). The intervals between these peaks were almost 190 constant, but their timing relative to the onset of bursts shifted depending on the burst phase ( Fig  191   5C-E). When the burst phase increased from 0° to 360°, the peaks moved closer to the burst 192 onset (Fig 5C-E), and the timing of the major ON peaks showed approximately linear 193 dependence on the burst phase (Fig 5C-E). These results indicate that the phase of 3-spike bursts 194 encodes the temporal phase of an oscillatory component. In addition, we found that the distance 195 of the point ( " , # ) from the origin of the " -# plane encodes the amplitude of the 196 oscillatory component (Fig 5F). Similar coding was found for bursts elicited with natural scene 197 movies (Fig 6). indicate the average and SEM values among 8 cells that generated more than 1000 3-spike 228

bursts. 229
As shown in Fig 4B, bursts with the phase 0° and 180° had only ~5-and ~1-ms differences in 230 ISI1 and ISI2, respectively. Nevertheless, for bursts with the phase 0°, the ON and OFF peaks in 231 the encoded sequences were separated by ~80 ms, while the separation was only ~50 ms for 232 bursts with the phase 180° ( Fig 5A). These results suggest that single spikes in bursts do not 233 simply indicate the occurrence of a single characteristic in light intensity sequence. To further 234 confirm this point, we conducted a simple reconstruction analysis. We calculated the spike-235 triggered averages (STA) for three-spike bursts and then generated an estimated stimulus 236 sequence by adding the STA aligned according to the burst spikes (Fig 5G). This reconstruction 237 failed to reproduce the observed dependence of the ON and OFF peaks on the burst phases 238 (compare Fig 5A and H), and the deviation from the 3-BTA was much smaller than that of the 239 actual data (compare Fig 5E and I). To quantify the amplitude of the deviation, the root mean 240 square of the deviation from the 3-BTA was calculated for the period between −200 and +25 ms 241 and averaged for all burst phases. The values were significantly larger for the actual stimuli than 242 for the reconstructed stimuli (0.107 ± 0.061 for actual stimuli and 0.003 ± 0.002 for 243 reconstructed stimuli, n = 19, P = 1.5 × 10 −7 , two-tailed Mann-Whitney-Wilcoxon test). Thus, 244 the simple reconstruction model failed to explain the observed burst coding. 245

Two independent components of burst patterns 246
The oscillatory component encoded by 2-spike burst ISIs had peak-to-peak intervals that are 247 similar to those of the components encoded by 3-spike bursts (~80 ms; compare Fig 3G and 5E), 248 suggesting that 2-and 3-spike burst ISIs are modulated by related stimulus features. To further 249 characterize this similarity, we analyzed events in which both 2-and 3-spike bursts were 250 generated (e.g., Fig 4A, bottom). For each of these events we calculated the average values of " 251 and # for 3-spikes bursts and the average value of the 2-spike burst ISIs (Fig 7A). Plotting the 252 data on the " -# plane showed that 2-spike burst ISIs differ systematically depending on the 253 position of the events along the direction of ~30° (Fig 7B). A linear fitting indicated that the 254 optimum direction was 33.1° ± 15.6° (circular mean ± SD, n = 41; Fig 7C,  The above result raises the hypothesis that the retinal mechanism that modulates 2-spike ISIs 277 also modulates 3-spike burst patterns in the orientation of ~33.1° on the " -# plane. Since " 278 and # are suggested to have two degrees of freedom, one possibility is that another 279 independent mechanism modulates 3-spike burst patterns in the orthogonal orientation, i.e., 280 ~123.1°, on the " -# plane. If two such independent mechanisms were present, modulation of 281 " and # in the two orthogonal orientations would have independent trial-to-trial variations, 282 and we tested this prediction. Although " and # had an approximately circularly symmetric 283 distribution (Fig 4E), their trial-to-trial variations within each event had an asymmetric principal axes corresponding to the smaller and larger variances were in the orientations of 29.4° 286 ± 7.3° and 119.4° ± 7.3° (n = 41; orange and green in Fig 7C and D). This analysis indicates that 287 modulations of " and # in these two orientations have approximately independent trial-to-288 trial variation. Because these axes (~29.4° and ~119.4°) are close to those proposed for the 289 hypothesis (~33.1° and ~123.1°, see Fig 7C), the results are in accordance with the presence of 290 two independent mechanisms, one of which (~30°) is common to 2-and 3-spike bursts.  Fig 3G, yellow). Given that the total distribution of " and # is circularly 294 symmetric, the smaller trial-to-trial variability of the ~30° component as compared with the 295 ~120° component indicates that the former is more precise. Therefore, the component common 296 to 2-and 3-spike bursts is more informative than that specific to 3-spike bursts. 297 The two approximately independent components encode oscillatory components that are 298 approximately orthogonal to each other, i.e., components with ~1/4 cycles difference in phase 299 ( Fig 7F). This result suggests that the two orthogonal oscillatory stimulus components modulate 300 3-spike burst patterns in the orthogonal orientations on the " -# plane. This model is 301 consistent with the above result that 3-spike burst patterns encode the amplitude and phase of an 302 oscillatory component, since an oscillatory sequence with an arbitrary amplitude and phase can 303 be approximated by a sum of two orthogonal oscillatory components with fixed phases and the 304 same frequency. 305

306
Our results reveal that intraburst ISIs of retinal bursts encode oscillatory light intensity sequences 307 that are much faster than the sequence encoded by the spike number. When a burst has three 308 spikes, the two intraburst ISIs combinatorially encode the amplitude and phase of the oscillatory 309 component. These results therefore suggest that a k-spike burst (k = 2, 3) encodes k different 310 stimulus features by exploiting all the k degrees of freedom, i.e., the spike number and k−1 ISIs. 311 This simultaneous representation of multiple stimulus features enables multiplexed information 312 coding, a mechanism that greatly increases the information transmission capacity [19,27,28]. 313 Whether this combinatorial coding occurs for bursts with k ≥ 4 spikes remains unknown, as 314 these bursts were rarely observed under our experimental conditions (Fig 1D). 315 When a burst has three spikes, the two mechanisms combinatorially determine the two ISIs. 323

Mechanisms of the combinatorial ISI coding 316
Because the two mechanisms are driven by the two orthogonal oscillatory components, the two 324 few milliseconds are much more effective in eliciting LGN spikes than those with ISIs of > ~20 345 ms [39][40][41][42][43][44][45]. Consistently, LGN burst ISIs are sensitive to the millisecond-scale structure of 346 current input [19]. Similar to synaptic connections from the retina to the LGN, those from the 347 LGN to cortical neurons are more responsive to spikes with short ISIs than those with long ISIs 348 Therefore, it is possible that burst ISIs encode spatial information as well as temporal 355 information. 356

Conclusions 357
The present results suggest that the retina employs mechanisms to regulate multiple components 358 of intraburst ISIs, and thereby encodes multiple stimulus features by exploiting all degrees of 359 freedom of burst spike patterns, i.e., the spike number and multiple intraburst ISIs. This burst 360 coding is likely to affect visual information transmission, as synaptic transmission is sensitive to 361

Animals 365
All experiments were approved by the RIKEN Wako animal experiments committee and were 366 performed according to the guidelines of the animal facilities of the RIKEN Center for Brain 367 Science. Larval tiger salamanders were provided by Charles D. Sullivan Co. Inc., Nashville, 368 Tennessee, USA. 369

Recording and Stimulation 370
Retinal recording was performed as described previously [48]. Dark-adapted retinae from larval

Identification of Bursts and Events 386
Histograms of the ISIs were generated for each isolated ganglion cell. The histograms often had 387 two distinct peaks (Fig 1B) representing shorter and longer ISIs, corresponding to intra-and 388 interburst ISIs, respectively [3][4][5]. The threshold interval Tthresh was set at the trough between the 389 two peaks in the ISI histogram (Fig 1B). Tthresh was 38.6 ± 20.0 ms (mean ± SD) for 41 cells 390 stimulated with the spatially uniform stimulation, and 87.7 ± 18.5 ms for the 16 cells stimulated with 391 the natural scene movie. If two consecutive spikes occurred with an interval shorter than Tthresh, they 392 were incorporated into the same burst, while they were separated into two consecutive bursts if the 393 interval was longer than Tthresh ( Fig 1C). The robustness of this method was examined as follows. 394 Bursts were defined using various threshold intervals of ~10 ms to ~100 ms, and the rates of the 395 isolated spikes and bursts with 2-7 spikes were measured (Fig 1D). Events were determined as follows. The first spikes of bursts were extracted (Fig 1E, top) and 405 merged across different repeats of the stimulus presentation (black in Fig 1E, middle). The intervals were composed of short intervals representing inter-trial fluctuation (red in Fig 1E, middle, and red  409 in Fig 1F) and longer intervals separating the consecutive events (blue in Fig 1E, middle, and blue in 410 Fig 1F). Thus, if two consecutively merged first spikes were closer than Tthresh, they were 411 incorporated into the same event; otherwise, they were assigned into two consecutive events (Fig 1E,  412 bottom). The robustness of the method was evaluated as follows. If bursts occurred with a large 413 timing jitter in different repeats, two consecutive bursts in one repeat were incorporated into one 414 event. However, this occurred for only 2.7% ± 2.7% of bursts for the cells (mean ± SD, n = 41) 415 stimulated with the spatially uniform stimulation, and 3.5% ± 1.1% for the cells stimulated with the 416 natural scene movie (n = 16). The values were small, indicating that the definition of events was 417 robust. 418

Experimental Design and Statistical Analyses 419
The numbers of the analyzed cells and retinae were as follows. Salamanders: n = 41 cells in 15 420 retinae for spatially uniform stimulation, and n = 16 cells in 4 retinae for natural scene 421

stimulation. 422
Correlations between the spike number and the intraburst ISIs of the bursts were investigated 423 using data from ganglion cells stimulated with the spatially uniform stimulation as follows. For 424 each event j in which at least one 2-spike burst occurred, the average number of spikes (n (j) ) and 425 the average intraburst ISIs of the 2-spike bursts (m (j) ) were determined. The correlation between 426 n (j) and m (j) was then calculated across events. Similar correlations were calculated for ISI1 and 427

ISI2 of the 3-spike bursts. 428
Mutual information conveyed by the burst spike number was determined as follows. When 429 the stimulation was repeated, bursts occurred in a limited number of discrete events that were 430 defined at specific time points of the sequence (Fig 1A). Therefore, it was investigated how far 431 the receiver of bursts can specify the events by knowing the spike number, as compared to the 432 case where the receiver receives bursts without knowing the spike number. Rep  , and the total number of bursts is Burst = ∑ Burst (4) @ Ev 4?" . When a burst 438 was generated, the probability that the burst was in the j-th event is Burst Coordinate transformation of 3-spike burst ISIs was performed as follows (Fig 4B-D). For 3-460 spike bursts in event j, the average and standard deviation of ISI1 were designated as ISI RRRR " (4) and 461 SD " (4) , respectively. SD " (4) was linearly fitted with ISI RRRR " (4) ( = 1, ⋯ , Ev ) as SD " (4) ≅ 462 " UISI RRRR " (4) − " W, where " and " are constants (Fig 4C, top). For a variable " = 463 log "-(ISI " [ms] − " [ms]), \]\ ISI^= (ISI " − " ) ," (log 10) ," , and thus the standard deviation 464 of " in each event was similar among different events (see bursts shown in different colors in 465 Fig 4D). Bursts with ISI 1 ≤ " were removed from the analysis. # was similarly defined for 466 ISI # (Fig 4C, bottom). The principle axes of the distribution of " and # were determined 467 (Fig 4D), and new variables " and # were defined by scaling " and # along these axes 468 so that the standard deviations along these axes were equal to 1 (Fig 4E). " and # were 469 shifted so that their averages were zero. The burst phase was atan2( # , " ) ( Fig 4E). 470 To characterize the stimulus features encoded by bursts, the stimulus sequences preceding 471 bursts of a specific spike number, intraburst ISIs within a specific range, or burst phase within a 472 specific range, were collected. The average and SEM values of the collected sequences were then 473 used for the analyses. Neurons that generated only small numbers of 2-or 3-spike bursts were 474 removed from the analyses (see legends for Figs 3G, 5D, 6C, and 6D). For linear reconstruction, 475 stimulus sequences preceding the spikes in three-spike bursts were collected and averaged. The 476 average sequence was divided by three and then used as the STA for the reconstruction (Fig 5G-477   I). 478

479
The central data and computer codes used in this paper are available the open science framework 480 database at: https://osf.io/29ect/. Other data and codes are available upon request. 481