Time-course of population activity along the dorsoventral extent of the superior colliculus during delayed saccade tasks

The superior colliculus (SC) is an excellent substrate to study functional organization of sensorimotor transformations. We used linear multi-contact array recordings to analyze the spatial and temporal properties of population activity along the SC dorsoventral axis during delayed saccade tasks. During the visual epoch, information appeared first in dorsal layers and systematically later in ventral layers. In the ensuing delay period, the laminar organization of low-spiking rate activity matched that of the visual epoch. During the pre-saccadic epoch, spiking activity emerged first in a more ventral layer, ∼100ms before saccade onset. This buildup of activity appeared later on nearby neurons situated both dorsally and ventrally, culminating in a synchronous burst across the dorsoventral axis, ∼28ms before saccade onset. Stimulation of individual contacts on the laminar probe produced saccades of similar vectors. Collectively, the results reveal a principled spatiotemporal organization of SC population activity underlying sensorimotor transformation for the control of gaze.

and ventral channels. The most notable difference between VG and MG trials is the smaller, but not 135 significantly different (Wilcoxon rank-sum test, P>0.067, for all channels), peak amplitude of movement 136 activity for MG trials. Overall these results show that the amplitudes of two bursts are systematically 137 organized across depths. The peak visual activity is situated between 0.75mm and 1.05mm (between 138 channels 5 and 7) more dorsally than the peak movement activity, reflecting a visual preference for dorsal 139 channels and movement preference for ventral channels.  bin 5. Notably, these channels were the same that discharged maximally for the visual burst (thick black 184 trace). For the deepest channels, we noticed a difference between the two tasks. For MG tasks, the activity 185 was slightly higher than the peak activity of the visual burst while it remained at minimum level for VG tasks.

186
Note that the confidence intervals were very large and overlapping across bins which makes this increase 187 of activity not statistically significant in our data. Further investigation is needed to reveal the possible cause  along with the associated variability could be large enough to average out latency differences across 198 channels. To circumvent this concern, we aligned data on the onset of the visual burst. Briefly, for each trial 199 and for each channel, burst onset was estimated using the Poisson surprise detection method 24 . The  Figure 6b plots the visual latency estimate of each channel centered on the reference channel 206 7, i.e., the channel obtained from the CSD analysis for aligning data in depths across sessions. Note that 207 the reference and alignment channels need not be the same. A cubic polynomial fit of the visual latencies 208 reveals a general spatiotemporal trend for this session from dorsal to ventral depths (dashed trace, R 2 =0.83, 209 p=0.001). Figure 6c reports the session-averaged estimates of relative onset latencies in the visual epoch  in visual burst onset from dorsal to ventral depths within SC, spanning 7.3ms (between channel 5 and 219 channel -8) and 6.8ms (between channel -5 and channel -7) for VG and MG tasks, respectively. This is 220 consistent with single-synaptic transmission from one channel to another in depth, although other 221 mechanisms are also viable (see Discussion). A similar pattern in spike timing was recently reported across 222 cortical layers in rodents 25 .

223
Next, we analyzed the pre-saccadic epoch that leads to movement onset. Previous work has shown 224 that SC neurons display a buildup (or prelude) and/or burst of activity during that epoch 26,27 . In these 225 previous studies, buildup and burst activity were detected by threshold crossing of the averaged firing rate 226 activity computed over fixed temporal window defined for each event separately. This method was designed 227 with the only objective to detect the presence of these events. Here, we wanted to not only detect them but 228 also reliably estimate their onsets times and associated amplitudes.

229
We developed an algorithm that, first, estimates the onset latencies of well-defined events in the by a measure of reliability obtained through a bootstrapping procedure. Events that did not meet a reliability 237 criterion were discarded. As a consequence, anywhere between zero and three events were detected for 238 the average spiking activity of each channel. It is important to note that events E1, E2 and E3 capture 239 temporal characteristics of the activity and are interpretation-free in terms of neural circuit mechanism. A 240 subsequent classification procedure enabled the interpretation.
241 Figure 7a-i shows the detection of events E1, E2, E3 and P for the example dataset. The reliability 242 of the detection of each event was measured by the CIs obtained through bootstrapping (see Materials and 243 Methods). One can observe that for many channels CIs were well below the 0.6 threshold, indicating a 244 reliable estimation of all events. Panels f and h also report the slopes and CIs around the hinge points for 245 events E2 and E3. The slopes were not significantly different for most of E2 events, reflecting the slow 246 accumulation of the activity. In contrast, they were significantly different for E3, reflecting a sharp increase 247 of the activity. Figure 7i shows the final estimation of all four events, provided that they were statistically 248 significant (see Neuronal activity categorization, Materials and Methods). Even though the search window 249 for E2 was large (within 100ms before E1 and a total search window of 300ms), E2 events were 250 systematically detected within 50ms before E1, reflecting the onset of the accumulation that gives rise to 251 E1.

252
The next step was to classify these events into buildup or burst activity. When both E2 and E3 were 253 detected and were significantly different from each other (i.e., their CIs did not overlap), they were labeled 254 as buildup and burst events, respectively. If only E2 was reliably estimated, the activity displayed an early 255 buildup (<-50ms) but without a reliable hinge point for E3 (e.g., channel 3 of Figure 7i). When only E3 was 256 reliably estimated, the activity only displayed a burst just before saccade onset (>-50ms) but without an 257 initial buildup (e.g., channels 11 and 13 of Figure 7i). We plotted the distributions of events E2 and E3 258 pooled across all channels and bootstrap iterations, although cases when both or only one of the two events 259 were detected were considered separately (Figure 7k-m). It is important to note that the onset of each event 260 is only constrained temporally by the size of their respective search windows, which means that either onset 261 can occur at any time within 300ms before saccade onset with buildup preceding burst onset. The

262
distributions of event times are binomial and exhibit visual separation around −50 relative to saccade 263 onset. We therefore used this boundary to distinguish buildup (< −50 ) from burst (> −50 ) events.

264
Thus, many events detected as E2 (buildup) were reclassified as E3 (burst). Figure 7j replots the average 265 spike density functions of all channels with identified buildup and burst events superimposed. For this 266 session, peak activity (black tick marks) occurred around saccade onset for nearly all neurons along the 267 dorsoventral dimension, burst activity (purple tick marks) was also present in most neurons (but fewer than 268 those with peak activity), while buildup activity (cyan tick marks) was present primarily in neurons found 269 along the ventral half of the track.

270
The analysis shown in Figure 7 was performed for each laminar recording session. We then 276 Figure 8a show the results of the estimation of the events E1, E2, E3, and P. We found that 277 significant peaks of activity were detected across the whole track (from channel -8 to 7, Figure 8 a in Figure 8a,b indicate that events E1, E2 and E3 can be estimated reliably in the average waveform on 300 each channel and that they correspond to systematic events of the pre-saccadic activity.

301
Once the events (E2, E3) were detected, we classified them into buildup and burst categories.

320
The purple trace in Figure 8c shows the organization of burst onset across depths. Bursts were 321 found at all recorded depths, from channel -8 to channel 7. All bursts appeared on average -26.9ms (95%CI    345 Figure 8f plots the amplitude of spiking activity at the time of event P, E1, E2 and E3 across depths.

346
The activity at the onset of E1 was limited to low firing rate between 3.6spk/s ( Figure 8e,f indicate that events E1, E2 and E3 can be estimated reliably in the average 355 waveform.

356
The classification of events using classification labels (buildup, burst) are presented in Figure  to a true onset of activity relative to baseline.

366
The purple trace in Figure 8h shows the organization of burst onset across depths. Bursts were 367 detected on nearly every neuron encountered in the penetration, from channel -6 to channel 6.

374
Next, we analyzed the correspondence between the activity at burst onset and at the peak 375 response. To do so, we used the cubic fit computed on the distribution of the average firing rate at the onset 376 of the burst and at the peak P. A scaling factor was computed between the maximum firing rates of the two 377 fits. Figure 8d shows the rescaled fit of the bursts (green dashed line) for VG trials and shows the close 378 correspondence with the fit of the peak across all depths. Hence, the peak activity of the movement epoch 379 appears to be a scaled version of the activity at the burst onset across depths. Figure 8h shows the same 380 information for MG trials, and shows that the correspondence with the fit of the peak is close for ventral  Figure 9). Munoz and Wurtz found that burst-only neurons was the largest proportion (Table   395 1, 6%, 68% and 26%, respectively, from 26 ). Our results for both VG and MG trials show that buildup-burst Although the focus of this study was to examine SC activity patterns across depth, we did We used a multi-contact laminar probe to record simultaneously the activity of a population of 428 neurons along the dorsoventral axis of the SC of nonhuman primates performing delayed saccade tasks.

429
The categorization of the activity of each channel revealed, as summarized in Figure 11, a visual preference processes properties that were difficult to appreciate in previous studies that relied on single unit recordings.

446
These structural patterns of SC network architecture can inform the design of biologically-inspired models 447 that implement sensorimotor transformation 6,7 .

448
Past efforts to correlate activity features to neuron location in SC were limited by the uncertainties 449 of estimating the single electrode's depth. Even the most methodical approach e.g., 23 cannot account for 450 settling of neural tissue or of anisotropies in SC geometry (e.g., curvature). The constant intercontact 451 distance of the laminar probe combined with current-source density alignment circumvent many limitations  2015)) shows a linear trend of increasing motor dominance with depth, 456 but it isn't able to reveal the saturation of VMI at deeper locations that we report here ( Figure 3C). This

481
The transient burst of the visual response was followed by low-level activity in a subset of SC 482 neurons ( Figure 5). Given the rich balance of excitation and inhibition across all layers of the SC 35 , the 483 persistent activity can be readily generated through network dynamics e.g., 36 , although intrinsic, biophysical 484 features likely contribute as well 37,38 . That this low frequency activity is more prevalent in dorsal layers can 485 be inferred from data compiled with single electrode experiments see Figure 7C,D of 39 , although its 486 alignment, with rank order preserved, is best appreciated from the laminar probe data (compare Figure   487 3a,b with Figure 5). Notably, there was no transition from visually-dominant to motor-dominant layers during 488 the delay period, which we believe has major implications on reference-frame transformation research.

489
Previous studies have suggested that the transformation between reference frames occurs during the delay 490 period (under appropriate task design). Correlative evidence exists for craniocentric to oculocentric 491 representation 40,41 and for visual, target-centered to motor, gaze-centered coordinates 42-44 . One way to 492 reconcile these seemingly discrepant results, particularly with respect to the latter set of studies, is that the Once the animal receives permission to initiate a saccade, the population activity transitions to 499 more ventral layers (Figures 8 and 11), where neurons begin to accumulate activity ~100ms before saccade , but it has been debated whether the latter two are segregated along the dorsoventral axis. The ability 531 of a laminar probe to record simultaneously the activities of neurons along this axis revealed that there is a 532 gradual transition from visual to visuomotor to motor neurons with depth but the vast majority are visuomotor 533 (Figure 4). We also demonstrated previously 30 and discussed above that putative motor neurons can exhibit a visual response under certain conditions. Thus, we prefer to avoid making a distinction between 535 visuomotor and motor neurons.

536
Hypothesis-guided names have also been assigned to SC neurons, with different nomenclatures 537 being introduced over time. Saccade-related burst neurons discharge a high-frequency burst for optimal 538 vector saccades, and burst onset is tightly coupled to saccade onset, leading the movement by ~20ms 58 .

539
Neurons with low-frequency activity several hundred milliseconds before the burst were once called long-  were respectively 142±16ms and 207±20ms, relative to fixation offset). We sought to relate this finding to 555 the buildup or burst features of neural activity, and an initial glance suggests that accumulation onset is a 556 feature of buildup neurons. However, our analysis suggests that buildup onset occurs at least 30ms earlier, 557 ~100ms before saccade onset (Figure 8). We believe this discrepancy may be a result of the differences in 558 detection approaches. We developed a method to reliably detect and classify a buildup (accumulation) 559 and/or a burst (threshold) process. By contrast, the previous study, by using a two-piece linear regression, 560 limited the detection to only one neural event that spanned from 100 ms before fixation offset to the time of 561 peak in the saccade-related burst 65 . Their method was applied on individual trials and the amount of 562 stochastic noise in the spiking discharge may have prevented the distinction between two neural events.

563
Here we used trial-averaged waveform combined with bootstrapping of the trial sets, which allowed the 564 detection of at most two distinct neural events within a probabilistic framework. To the best of our 565 knowledge, this is the first time a method is developed to detect precisely and reliably the onset of buildup 566 and burst activity, beyond the use of predetermined temporal windows of analysis. This method revealed 567 that buildup activity was detected gradually on different channels spanning ~30ms since the initial buildup 568 around the center of the intermediate layers.

569
Burst onset occurred synchronously across all layers, ~28ms before saccade onset. This is

627
Neurophysiological signals were recorded using the Scout data acquisition system (Ripple, Salt 628 Lake City, USA). Data recording was synchronized with the beginning of the trial, and the timing of all trial 629 events were recorded simultaneously with raw neural activity. For each channel, the raw activity was parsed 630 into spiking activity (high pass filter at 250Hz and threshold at 3.5 times the RMS) and local field potential 631 (low pass filter at 250Hz). As the isolation of single neural activity was not achievable simultaneously on 632 every channel, the recorded spiking activity was always considered as being from multiple units in the 633 vicinity of each contact. A standard threshold crossing was used to determine spike times.

636
After recovery from surgery, each animal was trained to perform standard eye-movements tasks.

637
Eye movements were recorded using an infrared eye-tracker (EyeLink 1000 from SR Research, Ottawa,

668
Visual activity: The study of neural response to a stimulus is traditionally performed by aligning the 669 data on target presentation. This approach effectively ignores the trial-to-trial variability in both the actual 670 time of target presentation (due to frame rate limit) and neuronal stochasticity. We were concerned that this 671 variance may be greater than the time spanned for the visual response to emerge across layers. Therefore,

681
We next determined the relative burst onset times across the channels (Supplementary Figure 2).

682
First, each trial-averaged spike density waveform was baseline-corrected by the mean activity computed    Movement-related activity: Trial-averaged spike density waveforms were aligned on saccade 706 onset, which was detected using a velocity criterion (30deg/s) applied after 'go' cue. Each signal was also 707 corrected for baseline activity, defined as the average activity in the epoch [-100 0]ms relative to 'go' cue.

708
The movement-related activity was separated into two periods. The peri-saccadic or movement-related  Figure 3). In terms of stochastic accumulator 714 framework e.g., 65 , this is equivalent to detecting when the activity begins to accumulate, while accounting 715 for trends induced by baseline activity, and when it transitions into a burst. The objective was to detect up 716 to three events in this period: E1, corresponding to the time of significant change of activity compared to 717 baseline (which may include a linear trend, see below); E2, denoting when activity begins to accumulate 718 and corresponding to a 'hinge' point prior to E1; and E3, marking when the activity starts to burst and 719 corresponding to a 'hinge' point occurring between E2 and P, the time of peak activity around saccade 720 onset. To be able to flexibly detect one and/or two separate events (and subsequently onsets of buildup 721 and/or burst activity), the Poisson surprise method used before for the estimation of the visual onset latency 722 was discarded in favor of a 2-piecewise linear regression-based approach. We sought to limit the number 723 of ad-hoc parameters while relying on statistical measures of significance through data bootstrapping. Also, 724 this analysis was limited to saccades produced in the standard latency range of 200 to 400ms.

725
As an initial step in our analysis, we applied a detrending procedure to remove any potential bias 726 contributed by the low-frequency discharge from the delay period, well before the buildup and/or burst

815
Depth alignment of multiple sessions.

816
The above analyses focus on population neural activity across SC layers within a session. To 817 assess reliability, data must be averaged across sessions, which required appropriate alignment of data 818 collected in each penetration. To address this issue, we designed a method based on current-source 819 density analysis (CSD) 77 , which computes the second spatial derivative of LFPs and provides an estimate 820 of the distribution of the current sinks and sources as a function of space and time in a volume of tissue.

821
We estimated the CSD using the csdplotter toolbox that contains the implementation of the iCSD method

828
This feature was present in the recordings at both locations but translated in depth. The lower bound of the 829 sink pattern was at 1.28mm and 1.87mm in panels b and e, respectively. Such a strong sink appeared 830 systematically after target onset across all recording sites (data not shown here, but subject to a future 831 manuscript). We exploited this feature to align the data across sessions. The lower bound limit of the sink 832 pattern was used as a reference for estimating the relative depth of the probe. Supplementary Figure 1c  them. Thus, if any of these non-overlapping channels contains significant activity, then the total number of 850 channels across which data can be analyzed may be larger than the total number of channels of the probe.

851
In practice, only few channels were added and all data analysis that required depth alignment will be 852 presented between channels -8 and 8 (17 channels).

855
Trial-averaged activity was computed using large numbers of trials. Confidence intervals were 856 measured using bootstrapping (1000 bootstraps, bootci in MATLAB). Normality assumption was 857 systematically assessed using a Kolomogorov-Smirnov test (kstest in MATLAB). When the hypothesis of 858 normality was not rejected, a parametric test was applied (t-test, ttest in MATLAB) otherwise a non-859 parametric test was used (Wilcoxon rank sum test, ranksum in MATLAB). All tests were two-tailed.

860
To measure trends of latencies and activity amplitude across depths, we fitted a cubic function to 861 the data of the form