Abstract
Social communication is crucial for the survival of many species. In most vertebrates, a dedicated chemosensory system, the vomeronasal system (VNS), evolved to process ethologically relevant chemosensory cues. The first central processing stage of the VNS is the accessory olfactory bulb (AOB), which sends information to downstream brain regions via AOB mitral cells (AMCs). Recent studies provided important insights about the functional properties of AMCs, but little is known about the principles that govern their coordinated activity. Here, we recorded local field potentials (LFPs) and single-unit activity in the AOB of adult male and female mice during presentation of natural stimuli. Our recordings reveal prominent LFP theta-band oscillatory episodes with a characteristic spatial pattern across the AOB. Throughout an experiment, the AOB network shows varying degrees of similarity to this pattern, in a manner that depends on the sensory stimulus. Analysis of LFP signal polarity and single-unit activity indicates that oscillatory episodes are generated locally within the AOB, likely representing a reciprocal interaction between AMCs and granule cells. Notably, spike times of many AMCs are constrained to the negative LFP oscillation phase in a manner that can drastically affect integration by downstream processing stages. Based on these observations, we propose that LFP oscillations may gate, bind, and organize outgoing signals from individual AOB neurons to downstream processing stages. Our findings suggest that, as in other neuronal systems and brain regions, population-level oscillations play a key role in organizing and enhancing transmission of socially relevant chemosensory information.
SIGNIFICANCE STATEMENT The accessory olfactory bulb (AOB) is the first central stage of the vomeronasal system, a chemosensory system dedicated to processing cues from other organisms. Information from the AOB is conveyed to other brain regions via activity of its principal neurons, AOB mitral cells (AMCs). Here, we show that socially relevant sensory stimulation of the mouse vomeronasal system leads not only to changes in AMC activity, but also to distinct theta-band (∼5 Hz) oscillatory episodes in the local field potential. Notably AMCs favor the negative phase of these oscillatory events. Our findings suggest a novel mechanism for the temporal coordination of distributed patterns of neuronal activity, which can serve to efficiently activate downstream processing stages.
Introduction
Social, reproductive, and defensive behaviors are crucial for survival, and in many organisms their implementation depends on chemical signaling and processing. In most vertebrates, including mice, this function is primarily served by the vomeronasal system (VNS), a compact but essential sensory system dedicated to processing chemical signals from other organisms (Dulac and Wagner, 2006; Dhawale et al., 2010; Mohrhardt et al., 2018). Although the vomeronasal system resembles the main olfactory system, there are fundamental differences in physiology and function. Studies in behaving (Luo et al., 2003) and anesthetized (Hendrickson et al., 2008; Ben-Shaul et al., 2010) animals, and ex vivo preparations (Meeks et al., 2010) examined spiking responses of neurons in the mouse accessory olfactory bulb (AOB) during sensory stimulation. These studies showed that AOB neurons display low baseline rates, and typically respond to chemical stimuli from conspecifics and nonconspecifics with rate elevations (Tolokh et al., 2013; Kahan and Ben-Shaul, 2016). Anatomically, the first brain relay of the vomeronasal system, the AOB, appears like an integrative network, with its output neurons, AOB mitral cells (AMCs), sampling activity from multiple input channels (Wagner et al., 2006; Larriva-Sahd, 2008; Mohrhardt et al., 2018). While individual vomeronasal ligands can elicit physiological and behavioral responses (Kimoto et al., 2005; Ferrero et al., 2013; Li and Liberles, 2015; Ishii et al., 2017; Ishii and Touhara, 2019), most natural vomeronasal stimuli are complex bodily secretions which contain numerous chemical components. Correspondingly, information about natural stimuli is likely represented by AMCs in a distributed manner (Ben-Shaul et al., 2010; Tolokh et al., 2013; Hammen et al., 2014; Fu et al., 2015; Kahan and Ben-Shaul, 2016).
In the main olfactory system, information about odorants is also conveyed to downstream regions by the concerted activity of ensembles of output neurons (Malnic et al., 1999; Buck, 2004; Koulakov et al., 2007; Uchida et al., 2014; Chong and Rinberg, 2018), and several processes were shown to organize and format output patterns of main olfactory bulb neurons. These are often revealed by local field potential (LFP) signals, which typically reflect population-level activity (Buzsáki, 2006; Buzsáki et al., 2012; Einevoll et al., 2013). Specifically, a variety of olfactory bulb LFP oscillatory patterns, spanning the theta, beta, and gamma frequency ranges, were described previously (Kay, 2014, 2015; Gretenkord et al., 2019). Although the exact mechanisms underlying these oscillations are only partly understood (Bathellier et al., 2006; Fukunaga et al., 2014; Li and Cleland, 2017; Fourcaud-Trocmé et al., 2018; Burton and Urban, 2021; Kersen et al., 2022), it is known that they arise from both local (e.g., within the main olfactory bulb) and global interactions (e.g., involving multiple olfactory brain regions), and are directly influenced by breathing. Together, these oscillatory patterns reflect processes that may gate, bind, and convey parallel streams of information, as it is transmitted from the main olfactory bulb to downstream regions. More specifically, activity of neurons with respect to each other, or with respect to the breathing cycle, can be informative about the stimulus (Shusterman et al., 2011; Smear et al., 2011; Haddad et al., 2013; Smear et al., 2013; Cleland, 2014; Uchida et al., 2014; Wilson et al., 2017; Chong and Rinberg, 2018; Gill et al., 2020; Dalal and Haddad, 2022), suggesting that these relationships are important for chemosensory information processing. To date, similar phenomena were not reported in the vomeronasal system.
One of the most notable differences between the main olfactory system and the vomeronasal system involves stimulus uptake. In the main olfactory system, stimulus uptake is linked to breathing, and active sampling modifies sniffing behavior and hence the manner by which information is relayed to downstream regions (Kepecs et al., 2006; Verhagen et al., 2007; Wachowiak, 2011). In contrast, vomeronasal system stimulus uptake is largely uncoupled from breathing, and instead relies on active suction by the vomeronasal organ (Meredith and O'Connell, 1979; Wysocki et al., 1980; Ben-Shaul et al., 2010). Other differences between the systems include the intrinsic properties of olfactory bulb neurons (Zibman et al., 2011; Shpak et al., 2012; Gao et al., 2017), their interactions with each other (Takami and Graziadei, 1991; Dulac and Wagner, 2006; Larriva-Sahd, 2008), and the properties of downstream regions that receive the information (Davison and Ehlers, 2011; Haddad et al., 2013; Uchida et al., 2014; Dalal and Haddad, 2022). Together, these observations imply that coordinated patterns of population activity in the VNS, if found, may be different from those in the main olfactory system.
Previously, we and others (Zylbertal et al., 2015, 2017b; Gorin et al., 2016; Tsitoura et al., 2020) have identified infraslow oscillations in the AOB. In addition, a number of studies have shown a broad range of LFP AOB oscillations that can reflect experience, context, and potentially active arousal/sampling (Binns and Brennan, 2005; Leszkowicz et al., 2012; Tendler and Wagner, 2015; Pardo-Bellver et al., 2017; John et al., 2022). However, explicit relationships between these rhythmic population activity patterns and sensory stimulation, and, more importantly, their potential interaction with single-unit spiking activity are unknown.
Here, we describe prominent theta-band LFP oscillations in the AOB. We show that they are induced by stimulus presentation and that their magnitude is stimulus dependent. Notably, in addition to full-blown oscillations, the AOB network also displays a more subtle pattern of stimulus-dependent activity that closely resembles the oscillations in its spatial distribution. Our analysis suggests that, once sensory stimulation is sufficiently intense, these latent patterns develop into prominent oscillations. Based on their spatial and temporal features, we propose that these LFP patterns reflect a periodic exchange within the AOB network, likely involving AMCs and granule cells. Importantly, we show that these oscillations constrain spike times of individual AOB neurons and can therefore serve a thresholding, amplifying, and potentially also multiplexing function. Thus, despite prominent differences in their constituent elements and connectivity, it appears that, like the main olfactory system, the vomeronasal system includes circuit elements that produce population-level oscillations, which serve to orchestrate information transfer to downstream regions.
Materials and Methods
Subjects
For recordings, we used adult sexually unexperienced male and female mice from the BALB/C and ICR strains (Envigo Laboratories). All experiments complied with the Hebrew University Animal Care and Use Committee. Table 1 provides a detailed description of subjects used in each session.
Details of the sessions used in this article
Stimuli and dataset
The dataset contains recordings from 34 sessions. Of 41 sessions initially considered, 7 were discarded because they lacked oscillatory activity and clear modulations of neuronal activity. This can result from several factors including unsuccessful cuff electrode implantation or poor physiological condition of the subject mouse. Stimuli were collected from adult estrus and nonestrus sexually naive female mice of the BALB/C, C57BL/6, and ICR strains (Envigo Laboratories), and from castrated, naive, and sexually experienced male mice from the BALB/C, C57BL/6 and ICR strains. For collection of mouse urine, mice were gently held over a plastic sheet until they urinated. The urine was collected in plastic tubes, flash frozen in liquid nitrogen, and then stored at −80°C. Predator urine was purchased (PredatorPee, Maine Outdoor Solutions) and stored at −80°C. Three sessions include urine, saliva, and vaginal secretion stimuli. See Kahan and Ben-Shaul (2016) for details on stimulus collection for these stimuli. One session included urine stimuli with different pH values, previously described in the study by Cichy et al. (2015). Table 1 provides a detailed description of the stimulus sets used in each session.
Surgical procedures and electrode positioning
Experimental procedures were detailed in the study by Yoles-Frenkel et al. (2017). Briefly, mice were anesthetized with 100 mg/kg ketamine and 10 mg/kg xylazine, and tracheotomized, and a cuff electrode was implanted around the sympathetic nerve trunk. Incisions were closed with Vetbond (3M) glue, and the mouse was placed in a stereotaxic apparatus where anesthesia was maintained throughout the experiment (0.5–1% isoflurane in oxygen). In parallel penetrations, a craniotomy was made immediately rostral to the rhinal sinus, the dura was removed, and electrodes were advanced into the AOB at an angle of ∼30° with an electronic micromanipulator (model MP-285, Sutter Instruments). In perpendicular sessions, the craniotomy was caudal to the rhinal sulcus and electrodes were advanced at ∼60°. We used 32-channel probes (NeuroNexus) with either 8 channels on each of four shanks or 16 channels on each of two shanks (Table 1). Some sessions include LFP recordings from only 16 channels. Before recordings, electrodes were dipped in fluorescent dye (DiI, Thermo Fisher Scientific) to allow subsequent confirmation of placement within the AOB external cellular layer (ECL).
Stimulus presentation procedure
Each session contains several distinct stimuli (Table 1), and each stimulus was presented three to six times within a session (mean ±SD, 5.2 ± 0.78 repeats per stimulus; median and mode, 5 repeats per stimulus). Stimuli were presented in blocks, and within each block, all stimuli were presented in a random order. During each presentation, 2 µl of stimulus was applied directly into the nostril. After a delay of 20 s, a square-wave stimulation train (duration, 1.6 s; current, ±120 µA; frequency, 30 Hz) was delivered through the sympathetic nerve cuff electrode to induce vomeronasal organ suction. Following another 40 s delay, the nasal cavity was flushed with 1–2 ml of Ringer's solution, and drained with suction from the nasopalatine duct. The cleansing procedure lasted 50 s and included sympathetic trunk stimulation to facilitate stimulus elimination from the vomeronasal organ lumen.
Electrophysiology
Neuronal data were recorded using an RZ2 processor, PZ2 preamplifier, and two RA16CH head-stage amplifiers (all from TDT). LFP and single-unit data were saved as two parallel streams with different sampling rates and bandwidths. LFP data were collected at 508.63 Hz and bandpass filtered at 0.5–300 Hz. A notch filter was applied offline to remove 50 Hz line noise. Single-unit and multiunit activity was sampled at 24,414 Hz and bandpass filtered at 300–5000 Hz. Custom MATLAB codes were used to extract spike waveforms. Spikes were sorted automatically according to their projections on two principal components on eight channels of each shank using KlustaKwik (Csicsvari et al., 2003; Hazan et al., 2006) and then were manually verified and adjusted using Klusters (Hazan et al., 2006). Spike clusters were evaluated by their waveforms, projections on principal component space (calculated for each session individually), and autocorrelation functions. A cluster was defined as a single-unit if it displayed a distinct spike shape, was fully separable from both the origin (noise) and other clusters, and its autocorrelation function demonstrated a clear trough around time 0 of at least 10 ms. Clusters apparently comprising more than one single-unit were designated as multiunits. Under these definitions, multiunits could represent the activity of as few as two units. As a further criterion for single-unit classification, we examined the “shape symmetry” of all spike shapes. Shape symmetry was defined as the correlation coefficient of the two sides (margins) of the spike waveform around its peak (in absolute value, with one of the margins reversed). For example, if the waveform contains nine samples with a peak at sample 5, the correlation will be calculated between the vectors defined by samples 1–4 and a reversed version of the vector defined by samples 6–9. If the shape is symmetric, the symmetry will be 1. When the peak is not in the middle, as is usually the case, the correlation is calculated for vectors determined by the smaller margin. A high symmetry (>0.95) is helpful for identifying noise accidentally misclassified as spikes. Using a graphical representation of all spike shapes, we examined the entire dataset for such noise-like spikes and removed them from further analyses. Neuronal firing rates were calculated in 0.2 s windows and smoothed with the smoothdata function, (using a 1 s Gaussian window).
Experimental design and statistical analysis
Analyses were conducted in MATLAB, using built-in MATLAB functions when available, or custom code as needed. All italicized function names refer to MATLAB functions. Figures were generated in MATLAB, exported to Adobe Illustrator, and edited for clearer presentation.
The basic data-processing procedure, applied to each session separately, begins with manual definition of prominent oscillatory events. These events allow us to define a “template,” which is then used to provide an ongoing measure of network activity (the template match signature), and to derive a larger set of episodes with high resemblance to the manually defined oscillatory pattern. The spectral and temporal features of each of these episodes are then analyzed separately. The key stages in this procedure are listed as follows.
Manual definition of oscillatory events and derivation of the oscillatory template. For each session, oscillatory events were marked using an interactive custom-written data browser that allowed viewing LFP, unit activity, and trial events, and marking oscillatory epochs. These epochs were saved, and their average LFP envelope pattern was defined as the oscillatory template for that session. The LFP envelope was derived with the envelope function using the rms option, and a window of 0.1 s. The high stereotypy between all user-defined events (see Fig. 2C,H–J) implies that the oscillatory template is robust to small variations in manual event definition.
Derivation of the LFP oscillatory template match signature (TMS). The TMS is the linear correlation (corr function) between the oscillatory template and the LFP envelope, calculated at each time point (see Fig. 3). The result is smoothed (smoothdata function with a 0.5 s Gaussian window), rectified, and raised to the 10th power to enhance contrast.
Detection of candidate oscillatory episodes. All data segments in which the TMS was above its mean value (calculated for each session separately; see Fig. 3), continuously for at least 1 s, were considered for further analysis (selection procedure described in stage 5 below). The mean was used as a threshold since it was typically lower than thresholds obtained with other approaches, resulting in more detected episodes and hence larger samples for statistical analyses.
Spectral analysis of oscillatory events. Spectral analysis was applied to all candidate oscillatory episodes obtained in the previous stage. We used the wavelet synchro-squeezed transform (WSST), which includes a “reassignment” step that enhances temporal and spectral estimates by compensating for spreading effects of the base wavelet (Daubechies et al., 2011; Thakur et al., 2013). The WSST (wsst function) was applied to the LFP channel with the highest energy using the analytic Morlet wavelet family. Focusing on one channel is warranted given the high zero-lag correlation among all channels during individual episodes, as calculated by the cross-covariance (xcov function; Extended Data Fig. 5-1). Then, “ridges,” defined as continuous frequency sweeps in time, were extracted (wsstridge function with a penalty value of 5; higher penalty values favor ridges without abrupt changes in the frequency domain). Frequency contributions >40 Hz (which are well above our range of interest) were ignored. The first ridge extracted using this procedure contains the most energy, with subsequent ridges containing less energy. The inverse-WSST (iwsst function) reproduces the original time-domain signal from the output of the WSST or from a subset of ridges. The instantaneous phase is the angle (angle function) of the Hilbert transform (hilbert function) of the reconstruction using a single ridge. In our analysis, we focused on episodes that are well approximated by a single ridge, as determined with the reconstruction score, which is defined as the correlation coefficient between the original signal and its reconstruction (in the time domain). The mean frequency of a ridge is its weighted mean frequency across all bins, with weights given by the amplitude (absolute value of the WSST in the corresponding bin). Analysis based on ridges rather than rigid frequency bands is suitable for oscillatory events whose frequency varies in time. In addition, since a ridge contains a single frequency at each time point, each time point is associated with a well defined instantaneous phase. This stage yields a set of episodes along with their reconstruction scores, mean frequency, mean TMS values, and mean envelope. See Figure 4 for an illustration of this procedure and Extended Data Figure 4-1 for reconstruction of oscillatory episodes using one or three ridges.
Selection of episodes for further analyses. For further analysis, we selected episodes that fulfilled the following two criteria: a reconstruction score >0.5; and a mean frequency in the theta range (2–14 Hz). This range was selected after manual inspection of oscillatory events indicated that virtually all are within that range. We then manually excluded episodes whose high reconstruction score was not due to a dominant theta component (“false positives”) and identified “complete” episodes that display prominent noise-free oscillatory modes. Of 2609 candidate episodes, 266 were designated as false positives, and thus were excluded from further analysis, and 89 were designated as complete. The selection criteria ensure inclusion of episodes with a dominant theta-band component. Complete episodes were primarily used for current source density and independent component analyses (ICAs). When a distinction is made between all episodes and complete episodes, this is indicated in the text. The relationships between various features of individual oscillating episodes are shown in Extended Data Figure 4-2A–D. In several analyses, oscillatory episodes were associated with specific stimuli. An episode was associated with a stimulus if it was entirely included within a stimulus presentation trial (beginning at stimulus application and ending before the wash).
Current source density analysis.
Current source density is the second spatial derivative of the LFP signal. It was calculated for each complete episode along contact sites within a given shank (Mitzdorf, 1985). The time-evolving current source density was smoothed using the smoothdata function (20 ms, Gaussian window).
Analysis of the LFP and unit signals with respect to sensory events.
To study the relationship between the TMS and the sensory events, we considered stimulus presentation epochs beginning with stimulus application to the nostril and ending 40 s after electrical nerve stimulation (for a total period of 60 s). Some analyses and displays include margins of 5 s at the beginning and end of this epoch (total 70 s). For each trial, we derived the LFP template match signature and the spiking rates of each of the units (calculated in 0.2 s bins and smoothed with the smoothdata function with a 1 s Gaussian window). This allowed us to calculate the average TMS across all trials, and all trials for a given stimulus. In addition, we derived the LFP envelope associated with each stimulus, evaluated in time points during which the TMS was above threshold. Since TMS values vary across individual trials and stimuli, the reliability of stimulus-specific envelope estimates may vary across stimuli.
Correlation between unit firing rates and the TMS.
Spiking rates were resampled into bins matching the LFP sampling rate and smoothed using a 5 s Gaussian window (smoothdata function). The relatively long window was chosen to reduce noise in firing rate estimates, and thus in the correlation between the spiking and LFP signals. The correlations between spike and LFP envelope signals were calculated using the xcov function (normalized scaling).
Spike time phase analysis.
For each unit, spike times were assigned with the phase of the closest LFP sample within oscillatory episodes. Shuffled phase distributions were obtained by randomly reordering the phases associated with each episode (randperm function). The random permutation was applied independently to each unit within each of the episodes, thus abolishing phase dependencies of simultaneously recorded units. For both real and shuffled data, we obtain spike phase distributions of each unit for each of the stimuli, and a composite phase distribution is obtained by pooling across all stimuli. For statistical analyses of phases, we used the circ_stat MATLAB toolbox (Berens, 2009); specifically, circ_otest for uniformity of data, circ_mean for the circular mean, and circ_cmtest for a nonparametric comparison of medians of two phase distributions. Across all of these analyses, we only evaluated phase distributions (whether for individual neurons or for specific unit-stimulus combinations) if the number of spikes within a sample was at least 50.
Pairwise comparisons of phase distributions.
For each pair of (simultaneously recorded) units, we conducted a comparison if at least one of the two units under comparison (or their shuffled distributions) displayed a nonuniform phase distribution (p < 0.05 as determined using the circ_otest). The same procedure was applied for comparison of phase distributions associated with pairs of distinct stimuli for a given unit. Namely, we only compared distributions if at least one of the two distributions (in either the real or shuffled data) significantly deviated from uniformity (p < 0.05, using the circ_otest). The criteria were applied to avoid comparisons between two uniform distributions (which are meaningless), and to equalize the selection criteria for the real and shuffled distributions. Our results hold with unequivocal differences between real and shuffled data using other criteria for inclusion (e.g., when considering all samples regardless of deviation from uniformity, or conversely, when considering only pairs for which both samples significantly deviate from uniformity).
Assessing the effect of phase locking on downstream neurons.
If the firing rate is sufficiently small, a Poisson process with a baseline rate of F Hz has a firing probability of F/1000 within a 1 ms window. Thus, the probability to obtain x or more spikes from n neurons within a temporal window of s ms, is given by the upper tail of the binomial cumulative distribution function with parameters, x, n · s, and F/1000 [binocdf(x, n · s, F/1000, “upper”)].
Independent component analysis.
To evaluate “sources” (components) of individual oscillatory episodes, we used the MATLAB rica (reconstruction-independent component analysis) function. We specified eight components (a larger value resulted in additional noisy sources) and an iteration limit of 1 · 106 (all runs converged before this limit). The transform function was then used to reconstruct the components from the model structure returned by rica. Because ICA cannot determine source scales (magnitude), we normalized all to a maximum (absolute) value of 1. The contribution of each component to each channel is evaluated by projecting the original LFP signal (each channel separately) on each of the sources derived for that episode (see Figs. 9, 10, for components and their projections). Like the scale, the order and the sign of the components is arbitrary.
Phase-locking value.
The phase-locking value (PLV), either across all spikes from a neuron or of mean phases across neurons (Lachaux et al., 1999), was calculated as follows: PLV = |Σei·θj|/n, where the sum is across all (n) unit vectors with phases θj, i is the imaginary unit, and ‖ indicates vector magnitude. Intuitively, the PLV is the ratio between the length of the vector derived by adding all unit vectors with their actual orientations (phases) θ, and the length of the vector that would be obtained if all (n) unit vectors shared the same orientation. To calculate whether the phase-locking value of the neuronal population is larger than expected by chance, we randomly drew n phase values (number of single-units in the population) from a uniform distribution and calculated their phase-locking value. This procedure was repeated 1000 times, allowing us to assess the probability to obtain the observed phase-locking value by chance.
Synchronization among units.
To determine the degree of temporal synchronization within oscillatory episodes, we divided each episode to temporal windows of varying resolutions ranging from 1 ms to 2 s (see Fig. 11). Then, for each episode, we counted the total number of spikes per bin, combined across all single-units. This set of counts defines a binned spike vector, whose (squared) magnitude is given by R = ΣCj2, where the sum is over counts in each of the temporal bins. Note that the number of bins depends on the duration of each episode and the bin size. We only considered cases with at least two bins (a limiting factor with brief episodes and large bins), and in which the number of spikes was larger than the number of bins (a limiting factor with small bins). The magnitude (R) reflects synchronicity, since for a given number of spikes, it increases when spikes are concentrated in a few bins (the largest magnitude occurs when all spikes are concentrated in a single bin, and the smallest when spikes are distributed equally between bins). Then, for each episode, we repeat the same calculation 5000 times with shuffled spike times. Shuffling is achieved by a random (circular) shift of spike times of each unit within each episode (see Fig. 11). The circular shift conserves the firing statistics of each unit, but disrupts its relationship to other units (and to the oscillations). We then count the fraction of shuffled vectors that are larger than the observed vector (R) and use it as an estimate of the probability to obtain the observed degree of synchronization by chance (for that episode). Low probabilities thus support synchronicity. Finally, we count the number of episodes for which this probability is smaller than a certain threshold (0.05). Under the null hypothesis (i.e., spikes are not coordinated across neurons), this probability is given by the upper tail of the cumulative binomial distribution. Thus, for each temporal resolution, it can be considered as the probability to obtain this number of significant cases of synchronization by chance (i.e., significance). In practice, because of the above-mentioned inclusion criteria, the smallest bin size that yielded enough spikes for analysis was 10 ms.
Data availability
All the raw data, processed data, and codes will be provided by the authors on request.
Results
To investigate AOB LFPs during sensory processing, we presented vomeronasal stimuli to anesthetized mice by direct application to the nostril, and electrically activated the vomeronasal organ (Fig. 1A). During each session, we present several distinct stimuli. Each stimulus is presented three to six times in an interleaved order. During stimulus presentation, we record both LFP (0.5–300 Hz) and single-unit activity using planar multielectrode arrays, which sample a substantial extent of the AOB. Our dataset contains 34 recording sessions from adult male and female mice. Details about each session are given in Table 1. We first describe recordings made with 32-channel probes, advanced in parallel to the AOB ECL (Fig. 1A–C), thereby maximizing the number of contact sites within the ECL, which contains AMC somata (Larriva-Sahd, 2008).
Experimental setup and recording sites. A, Experimental setup showing electrode array. The lightning bolt indicates stimulation of the sympathetic nerve trunk, which innervates the vomeronasal organ. During each trial, a drop of stimulus (2 µl) is applied to the nostril (application event), and, after 20 s, the vomeronasal organ is activated via a cuff electrode (stimulation event) placed around the sympathetic nerve trunk. After 40 s, the nasal cavity is flushed and another stimulus is presented. B, Schematic of the 4 × 8 recording array. Intersite distances along a shank are 50 or 100 µm. C, Atlas images of the AOB in sagittal (left) and coronal (right) views. Images are adapted from Paxinos and Franklin (2001). The black line in the sagittal section shows the approximate orientation of the recording probe. The green line shows the approximate extent of the recording region along one shank (with 100 µm intersite spacing). Histologic image shows the path of a DiI-stained probe (red), with nuclear DAPI counterstaining (blue). The broken white lines delineate the AOB ECL. The square within the coronal section shows the approximate horizontal extent of the probes (600 µm). ECL, External cellular layer; EPL, External plexiform layer; MCL, mitral cell layer; GL, glomerular layer; GCL, granule cell layer; lot, lateral olfactory tract; VN, vomeronasal nerve layer; VNO, vomeronasal organ. Table 1 provides details of recording sessions.
The AOB displays a spatially stereotypical pattern of oscillations
To detect prominent LFP patterns, we first manually browsed the multichannel data. Occasionally, we observed distinct oscillatory events that appeared across a subset of channels at varying amplitudes (Fig. 2A,B). As detailed below, these oscillatory episodes vary in frequency (both within and across an episode), but broadly fall within the theta range (2–12 Hz). Notably, within a session, all manually detected oscillatory events shared a similar across-channel envelope pattern (Fig. 2C,D,H,I). Plotting LFP signals from each channel according to their relative spatial locations reveals that regions with high LFP activity correspond to a central active area, surrounded by a region with considerably lower magnitudes. This is shown for one oscillatory event in Figure 2E, and for the average envelope across all manually detected events (from the same session) in Figure 2F. As electrode positioning within the AOB is variable, these envelope distributions are distinct for each recording session. Yet, almost all reveal a central hotspot of activity. Examples from different sessions are shown in Figure 2G and Extended Data Figure 2-1. Considering AOB anatomy and probe dimensions (Fig. 1C), these patterns imply that signals are generated locally and do not represent volume conduction from other regions. Further evidence for this claim is provided below.
Basic features of LFP data. A, LFP signals from 32 channels during one trial (∼70 s). Stimulus application is indicated by the black vertical line and drop icon, vomeronasal organ stimulation is indicated by the red vertical lines and the lightning bolt icon, and the rightmost vertical blue line represents the beginning of the wash period between trials. Channels are arranged according to their serial numbers (not spatial positions). Channel colors indicate the shank of the recording site (as in Fig. 1B). LFP amplitudes are normalized to fit within the plot. A vertical distance of 1 (i.e., difference between channel baselines) corresponds to 357 µV. B, Magnified view of a section shown in A. C, LFP envelopes (across all channels) associated with all manually defined events for the same session from which A and B are taken. D, The mean LFP envelope across all events in C. E, A section containing the oscillatory episode in B, with channels plotted according to their physical arrangement. In this session, horizontal distances (between shanks, columns) are 200 µm, and vertical distances among sites are 100 µm. Channel numbers are indicated in each panel and correspond to the numbers in A and B and Figure 1B. F, Mean LFP envelope for the session that includes the oscillatory template in E. The mean is calculated across all manually defined oscillatory events and is shown as a heat map according to relative spatial positions of channels. G, Mean LFP envelopes for four other sessions (Extended Data Fig. 2-1, envelopes in all recording sessions). In F and G, numbers above each image indicate the low and high color bar values. H–J, Stereotypy of manually defined events. H, Distribution of cross-correlation coefficients among all manually defined events within a session, pooled over all sessions (n = 2017 pairs). I, Distribution of the mean cross-correlations among events, for each session (n = 34). J, Distribution of the number of manually defined events per session (n = 34).
Figure 2-1
Spatial LFP envelopes over channels (i.e., the template), for each of the sessions in the dataset. Each image is normalized separately from minimal to maximal values, spanning the color map shown in the bottom right. The ranges are indicated for each session in µV. For example, the top left panel ranges from 0 to 150 µV. Note that values of 0 correspond to absent channels (the first 11 recordings only include 16 channels). Note also that some recordings contain nonfunctional recording sites. (i.e., 14–16 in some of the recordings in the lowest row). These broken channels were excluded for oscillatory template calculation. Download Figure 2-1, TIF file.
Detection of oscillatory episodes using the TMS
To automatically and objectively detect oscillatory activity, we used the observation that oscillatory episodes within a session are characterized by a stereotypical energy distribution across recording sites. We use the term TMS to denote the similarity to the energy distribution during oscillatory episodes. To derive the TMS, we calculate the linear correlation coefficient between the oscillatory template and the LFP envelope across all channels at each point in time (Fig. 3). The result is smoothed, rectified, and sharpened to increase contrast, yielding a measure of similarity between ongoing network activity and the oscillatory template. As shown in Figure 3D, the TMS captures oscillatory network activity and avoids periods of high power not due to LFP (e.g., during electrical stimulation), which exhibit a distinct spread of energies across channels. The TMS can be viewed in two complementary ways. First, it can be used to detect discrete episodes of oscillatory activity. Second, it can be treated as a continuous physiological signal. We first use the TMS to study discrete oscillating episodes (Figs. 4–6) and then continue to analyze it as a continuous measure of network activity (Figs. 7, 8).
Derivation of the template match signature. A–C, Illustration with simulated data. A, Manually detect and mark all oscillatory periods (of which one is shown and marked by the rectangle). B, Average the LFP envelope associated with all detected periods. This yields the template, which is shown linearly (corresponding to channel numbers) and spatially (according to relative locations of recording channels). C, Calculate the template match signature by correlating the oscillatory template with LFP energies at each sampling point. Here, the template is shown at three time points. In practice, it is correlated with each and every time point yielding a continuous output. D, Application of the procedure to real data. The template is shown at the inset. The TMS is given by the blue trace at the bottom of the panel. Note the low TMS values during electrical stimulation of the sympathetic nerve trunk (marked by the red vertical lines). This is because during these epochs, the spread of energies (induced by stimulation artifacts or other experimental noise) is not correlated with that occurring during oscillatory episodes. E, Application of the template to an entire recording session. The red line indicates the threshold value during this session. The bottom trace shows the thresholded TMS. F, Expanded view of a section in E.
Analysis of oscillatory episodes and their varying spectral and temporal features. A, An oscillatory event (one channel). Original signal (black) and reconstruction from the first extracted ridge (green). B, Time-evolving frequency (blue) and amplitude (green) of the first ridge for the episode in A. C, Instantaneous phase of the first ridge. D, Examples of spectral analysis of 12 oscillatory episodes from two recordings. Here, white traces show the original signals, with amplitude scale indicated on the right side. The colored trace is a spectrogram showing the time-evolving frequency, with power indicated by color. The same vertical axes are used for all episodes, but horizontal (time) scales vary. Note that many episodes show an increase and then a decrease in frequency. All episodes on the top row and the left episode on the bottom row are from the same session. All the other episodes are from another session. The episode in A is shown here on the second row and second column. E, Distribution of mean frequencies and durations for all complete episodes (n = 89). F, Like E for all episodes (n = 2343). See Extended Data Figure 4-1 for more examples of oscillatory episodes, and Extended Data Figure 4-2 for additional features of these episodes.
Figure 4-1
Complete episodes and their reconstructions. Each episode is shown twice. The top representation shows the original episode (black), and its reconstruction using 1 (green) or all 3 ridges (red). The bottom representation shows the original episode (µV scale on the right) and the first ridge time–frequency decomposition as a heat map. Consecutive episodes with the same background shade were recorded in the same session (i.e., the top row and the first example on the second row are from the same session). The session from which the episodes are derived is shown above each panel. Download Figure 4-1, TIF file.
Figure 4-2
General episode features. A, Reconstruction score (correlation coefficient between original signal and reconstruction with the first ridge) versus mean TMS (across the entire episode). B, Reconstruction score versus the mean first ridge frequency. C, Mean TMS versus normalized envelope amplitude (envelope within each session was normalized by the envelope of the episode with the largest envelope). D, Reconstruction score as a function of normalized envelope. In all panels, black dots represent all episodes, and black dots with red circles correspond to complete episodes. Note the correlation between mean template match and an above-threshold reconstruction score and the concentration of first ridge frequencies in the ∼5 Hz for the complete episodes. Download Figure 4-2, TIF file.
Oscillatory episodes comprise frequency sweeps within the theta range
We began by analyzing the spectrotemporal characteristics of the oscillatory episodes. The procedure for identification of episodes is detailed in Materials and Methods. First, we consider data sections for which the TMS is above threshold for ≥1 s (Fig. 3E,F). Then, for each of these, we select the highest energy channel and apply a wavelet transform to obtain its frequency domain representation. We then extract the dominant “time–frequency ridge” from this representation. Focusing on a single ridge facilitates characterization of the frequency and phase of oscillations as a function of time. Application of this procedure to one oscillatory episode is demonstrated in Figure 4A–C. Only cases for which a single ridge provided a good approximation of the original signal were considered (for details, see Materials and Methods). We distinguish between complete episodes (n = 89), which contain a complete uninterrupted oscillation, and all other episodes (n = 2254), which contain a dominant theta frequency component, but do not contain the entire event. Figure 4D shows examples of 12 complete episodes from two recording sessions. These examples underscore the diversity in episode durations and frequencies. Similar plots for all complete episodes are shown in Extended Data Figure 4-1. While many episodes follow a smooth up and down sweep in the frequency domain, others assume a more complex structure. Although individual events vary in mean frequencies, all fall within the broadly defined theta range (Fig. 4E, left; average of mean frequencies, 5.47 Hz). As seen for these examples (Extended Data Fig. 4-1), episodes also vary in duration. While most episodes are <5 s, some last ≥15 s (mean duration, 5.2 s; Fig. 4E, right). Distributions of episode durations and frequencies for all episodes (not only the complete episode) are shown in Figure 4F (n = 2343; average of mean frequencies, 7.1 Hz; mean duration, 2.3 s). In conclusion, although oscillatory episodes fall within the theta range, they vary in both duration and spectral characteristics. Notably, episode features also vary within individual sessions, indicating that variability is not simply due to different experimental preparations or electrode positions.
Spiking activity is locked to the LFP phase
In the main olfactory bulb, spike times are constrained to certain stages of the breathing cycle (Cury and Uchida, 2010; Shusterman et al., 2011; Smear et al., 2011; Fukunaga et al., 2014; Uchida et al., 2014). To investigate whether a similar phenomenon occurs in the AOB, we studied the relationship between spike timing and LFP phase. Importantly, the phase is derived from the channel with the highest LFP energy. This is justified by the fact that all channels are well aligned during oscillatory episodes (Extended Data Fig. 5-1).
By determining spike phase distributions and comparing them to phase-shuffled data, we can address several hypotheses. First, we tested whether AOB spiking activity depends on the LFP phase. Specifically, we tested whether the population of single-units displays phase distributions that deviate from uniformity, using the omnibus test for circular data (Berens, 2009). Figure 5A shows examples of phase distributions of seven single-unit with significant deviations from uniformity. To quantify the degree of phase locking, we applied the phase-locking value, which ranges between 0 (uniform phase distribution) and 1 (all phases are identical). Phase-locking values for the examples in Figure 5A range between 0.12 and 0.21.
Distributions of single-unit spike times relative to the oscillation phase are nonuniform. A, Phase histograms of seven single-units. The cartoon on the top left illustrates spike times with respect to the LFP phase. The number of spikes, PLVs, and p-values for tests of uniformity are indicated for each unit. Broken lines represent the LFP cycle for reference. Green bars show the phase distribution of units. Lighter bars show the distributions of the LFP phases across all episodes during the corresponding recording session and are flat as expected. B, Cumulative distribution of p-values for the test of uniformity (of spike times relative to LFP phase) across all single-units. Green and gray traces correspond to real and shuffled phase data, respectively. C, Distribution of phase-locking values for all single-units (top) and all tuned single-units (bottom). In both histograms, one data value >0.5 is not shown. D, E, Preferred phase distributions across all single-units (D) and tuned single-units (E; p < 0.05), in both linear (top), and polar (bottom) histograms. The red bars in the polar plots show the mean phase. Broken lines in linear plots show the LFP cycle for comparison. F, Like E, but for phase-shuffled data. G, Population phase-locking values (for all and for tuned single-units), and distribution of population phase-locking values when preferred phases are drawn from a uniform distribution (1000 shuffles). Vertical lines indicate observed phase-locking values. See Extended Data Figure 5-1 for an analysis of phase relationships of channels within episodes.
Figure 5-1
Analysis of temporal lags across channels within an episode. For each episode, we used the cross-covariance function (xcov function in MATLAB) to find the best match (highest covariance) between the channel with the maximum energy (for that episode) and all other channels. A, Distribution of all cross-covariance values obtained with this procedure. B, Distributions of cross-covariance values for complete episodes. C, Cumulative distribution of lags (in absolute values) that yield the highest cross-covariance value. To exclude nonoscillating channels from this analysis, we only considered channels that had a cross-covariance of at least 0.5 (a low threshold considering the distribution shown in A). The rationale for setting a threshold on the cross-covariance is that nonoscillating channels will not yield high correlations with the highest-amplitude, reference-oscillating channel D. Cumulative lag distributions when considering channels with cross-covariance values of at least 0.8. Mean values are indicated in C and D and are very small relative to a theta cycle. This analysis revealed that the percentage of channels with lags of <10 ms (in absolute value) are 85% across all channels in all episodes (no cross-covariance threshold), 98% for all channels with a cross-covariance >0.5 (as in C), and 99.9% when considering channels with cross-covariance >0.8 (D). A period of 10 ms represents one-twentieth of a 5 Hz theta cycle and thus assigning the oscillation phase based on one channel is expected to induce a negligible jitter in our phase estimates. CC, Cross-covariance. Download Figure 5-1, TIF file.
Comparison of phase distributions among units and stimuli (for a given neuron). Cartoons at the top illustrate the comparisons shown below. A, Cumulative distribution of p-values for pairwise comparison of phase medians across all simultaneously recorded unit pairs. B, Cumulative distribution of p-values for each unit and stimulus combination. C, Cumulative distribution of p-values for pairwise comparison of phase medians across all stimulus pairs for each unit.
Stimulus dependence of the TMS and the LFP envelope. A, An example of the TMS during an entire recording session. Horizontal green lines correspond to stimulus presentation periods (60 s), starting with application to the nostril and ending 40 s after nerve stimulation. Note that the TMS increases during stimulus presentation periods. B, Stimulus-triggered TMS, around all stimulus presentation times, in one type of experiment (Table 1, D2). The TMS rises following stimulus application (time 0), and then again during electric stimulation. The black vertical line corresponds to stimulus application. The two red vertical bars correspond to the start and end of the 1.6 s nerve stimulation. The blue vertical line marks the beginning of the wash period, which is effectively the end of the stimulation period. During electrical stimulation and washing, the TMS drops since electrodes record stimulation/flushing artifacts, respectively. C, Average stimulus-triggered TMS for different stimuli. Numbers represent fold dilutions. D, LFP envelopes associated with each of the stimuli. E, Mean LFP envelope across all sessions in this type of experiment (Table 1, D2). Each panel shows values for one session. The leftmost panel corresponds to the data in D; 100×, 30×, and 10× refer to dilutions of urine stimuli. See Extended Data Figure 7-1 for a statistical comparison of envelope magnitudes in this stimulus set, and for the same analysis applied to a different stimulus set.
Figure 7-1
The mean TMS, around all stimulus presentation times, in one type of experiment. Like Figure 7, B and C, but for a different stimulus set (Table 1, exp G). A, Stimulus-locked TMS, around all stimulus presentation times, in one type of experiment. B, Average TMS as a function of time for each stimulus separately. C, D, Histograms and statistical comparison of TMS values associated with different stimuli. C corresponds to the stimulus set shown here, and D corresponds to the stimulus set analyzed in Figure 7. Each panel includes data from all TMS samples associated with the particular stimulus. C includes data from 7 different sessions (n = 249,235, 7 sessions × 35605 samples per stimulus). D includes data from 5 sessions (178,025 samples, 5 sessions × 35,605 samples per stimulus). Median and mean values are shown in black and red, respectively. For visibility, some histogram bins are not contained within the plot limits (especially for bins near zero in C, and for large TMS values). In C, comparison of distributions using nonparametric ANOVA indicates a main effect of stimulus (11 df corresponding to 12 stimuli; p = 0 using the MATLAB kruskalwallis function). All pairs of stimuli are associated with significantly different distributions (p < 0.05), except pairs indicated with colored squares: panels that contain squares with the same colors are not significantly different from each other. All other pairs are different. In D, the significance for a main effect of stimulus is also 0 (8 df corresponding to 9 stimuli; n = 5 sessions), and all pairwise comparisons show significant differences as well (p < 0.05). Download Figure 7-1, TIF file.
Relationship between TMS and unit activity. A, LFP template match signature (blue) and summed neuronal activity signal (green) during an entire session. A temporally expanded view is shown below. Although both signals are correlated by virtue of their stimulus dependence (correlation of 0.565 at a lag of 0.076 s), closer examination shows that they do not perfectly overlap. Both single-unit and multiunit data are included in this panel. B, Maximal cross-covariance magnitudes, and lags at which they were obtained. Each dot represents one single-unit (n = 520). Medians of the cross-covariance (red line) and lags (green line) are both positive. C, Stimulus-triggered averages of the TMS and the summed single-unit activity (PSTH). Averages of both measures were calculated across all stimuli in each session. Each of the six panels represents one session. See Extended Data Figure 8-1 for additional examples.
Figure 8-1
Relationship between the TMS and summed single-unit activity, across stimulus presentations. Same conventions as in Figure 8C, but with more sessions included. Stimulus-triggered averages of both values were calculated across all stimuli. Green and blue traces correspond to the spike signal and the TMS, respectively. Note that the TMS generally follows the stimulation time course more closely and is more localized in time than the spike signal. Download Figure 8-1, TIF file.
Evidence for localization of oscillations within the AOB ECL (parallel penetrations). A, Schematic of a parallel penetration into the AOB. B, Spatial spread of LFP energies in one of the sessions (also shown in Fig. 1). C, LFP oscillations during one oscillatory episode, across all channels, plotted according to spatial position. D, Magnified temporal view of C, highlighting the phase similarity across all channels. E, Spike waveforms for the sites shown in B. Each square shows the spike waveforms recorded in that channel. Both single-unit and multiunit activity waveforms are shown. Each trace in each panel corresponds to one single-unit or multiunit. Although a given neuron can be recorded by multiple channels, here we plot the waveform on the channel with the largest amplitude. Note the correspondence between spiking activity and LFP amplitudes. F–I, Like B–E for another oscillatory episode from a different session. J, K, Independent component analysis of the oscillatory episodes shown in C and G, respectively. In each panel, the image on the left (components) shows eight components returned by the rica algorithm. All components are scaled to the same magnitude. The component resembling the oscillations is highlighted in red. The image on the right in each panel (projections) shows the projections of each of the components on each of the LFP channels. Note the pattern of weights for the main oscillatory component.
Evidence for localization of oscillations within the AOB ECL (perpendicular penetrations). A, Schematic of a perpendicular penetration into the AOB, using probes with a 2 × 16 configuration. Also shown is a dye tract of the electrodes in the session shown in B. B, Left, LFP signals during one episode recorded with a 2 × 16 electrode array. Note the reversal of the LFP along the probe (shaded regions). Right, Spike waveforms recorded on the sites shown on the left (including both single-unit and multiunit activity waveforms). C, D, Like B but for episodes from two other sessions. Note that in B, C, and D spiking activity is confined to sites at or above the LFP reversal. E, Current source density analysis of the areas indicated by gray squares in B, C, and D. The relevant sections are indicated by the italicized letters. These plots illustrate periods with alternating sinks and sources around the reversal point. Temporal scales in all plots are indicated by the green horizontal bars. F, G, H, Independent component analysis of the oscillatory episodes shown in B, C, and D, respectively. In each panel, the image on the left (components) shows eight components returned by the rica algorithm, all of which are scaled to the same magnitude. The component resembling the oscillations is highlighted in red. The image on the right in each panel (projections) shows the projections of each of the components on each of the LFP channels. Note the pattern of weights for the main oscillatory component.
Analysis of coordinated firing during oscillatory episodes. A, Illustration of the binning approach, for real (above) and shuffled (below) data. Note that the shuffled binned spike counts are a “circular” shift of the original. Numbers are color coded to make this more apparent. The magnitude of the binned vector (across all neurons active in each episode) is calculated and compared with magnitudes of vectors obtained by shuffling. B, Fraction of significant p-values across all episodes. A total of 2343 episodes was analyzed, but not all of them could be used for all bin sizes. The number of episodes actually tested is indicated in parenthesis below the bars. Bin sizes are shown in milliseconds, except 1 and 2 s (on the right). Probabilities denote the chance to obtain the observed number of significant cases under the null hypothesis of no synchronization (5%, indicated by the red horizontal line). Although we also analyzed windows as small as 1 and 5 ms, there were no valid episodes for these small bin sizes (see Materials and Methods).
Under our null hypothesis, spike times are independent of the LFP phase, and the distribution of p-values (for deviations from uniformity) across the population should be uniform, resulting in a cumulative distribution that follows the diagonal. As shown in Figure 5B, this is not the case. There is a marked excess of small p-values for the real data, but not for the shuffled data. The probability of obtaining an excess of small p-values (e.g., <0.05) is given by the binomial distribution, and in this case is very low (single-units, n = 331; p = 3.9 ·10−36). The corresponding p-value distribution for the phase-shuffled data is uniform, without evidence for an excess of small p-values (binomial p = 0.90). Thus, we first conclude that, as a population, AMCs display a clear spike time preference with respect to the LFP phase. The distribution of phase-locking values across the entire population of single-units is shown in Figure 5C for all units (mean, 0.1), and specifically for “tuned units” (mean, 0.15; defined as units with p-value < 0.05 for nonuniformity).
We next asked whether AMCs tend to prefer similar LFP phases. The preferred-phase distribution (Fig. 5D) again shows a marked deviation from uniformity (N = 331; mean phase, 207°; omnibus test for uniformity, p = 4.8 · 10−27). Similar results are obtained when considering only tuned single-units (Fig. 5E; n = 84; p = 1.3 · 10−11; mean phase, 203°). Notably, there is a strong preference for the negative LFP phase, consistent with a current sink in the ECL. For phase-shuffled data, preferred phase distributions do not deviate from uniformity (Fig. 5F). To quantify the degree of preferred phase similarity across units, we again applied the phase-locking value. The population phase-locking values for all units and for tuned single-units are 0.47 and 0.67, respectively. The probability of obtaining these values by chance (assuming a uniform distribution of preferred phases) is extremely low (Fig. 5G), indicating that there is a high degree of similarity among preferred phases across the population.
Phase comparison between individual units and stimuli
While preferred phases of distinct units are similar, they are not identical (Fig. 5D). If different units have distinct preferred phases, this could play a role in temporal coordination of specific subsets of neurons that share similar phases. To test this possibility, we conducted pairwise comparisons of phase distributions between simultaneously recorded units, using a nonparametric multisample test for equal medians (Berens, 2009). This analysis revealed a marked excess of significant differences among pairs (Fig. 6A; n = 1560 unit pairs; probability to obtain such an excess of p-values < 0.05 under the binomial distribution: 1.18 · 10−85; for the shuffled data, p = 0.038). Although the probability for the shuffled data is <0.05, it is orders of magnitude lower for the real data. We therefore conclude that, along with the general preference for the negative phase of the LFP cycle, there are frequent differences between phase preferences of individual units.
Finally, we test a more elaborate hypothesis, according to which individual units show distinct phase distributions when different stimuli are presented. First, we consider the population of unit–stimulus pairs (i.e., phase distributions of individual units during presentation of specific stimuli) and ask whether there is an excess of significant p-values (under the null hypothesis of uniform distributions). Indeed, this is the case for real data (n = 743 unit–event pairs; probability under the binomial model, 3.04 · 10−15; Fig. 6B), but not for shuffled data (p = 0.99). Finally, we compared phase distributions of individual units during presentation of different stimulus pairs. As shown in Figure 6C, there is a marked excess of significant pairwise differences for the real data (1675 comparisons; under the binomial mode p = 9.36 · 10−23), but much less for the phase-shuffled data (p = 0.03). This analysis thus supports the idea that the phase distribution of a given unit may also be stimulus dependent.
Together, we have shown that individual neurons in the AOB display preferred firing phases within the LFP theta cycle, and that the preferred phase varies among neurons and is stimulus dependent. This temporal structure of firing can play important roles in information transfer to downstream regions.
Oscillatory periods represent an extreme expression of a continuously varying network state that reflects sensory stimulation
Recall that the TMS can be treated as a continuous physiological signal that reflects the similarity to the stereotypical energy distribution associated with oscillatory episodes (Fig. 3). We now begin to explore the expression of this energy distribution throughout the experiment and, specifically, whether it is corelated with stimulus presentation. Figure 7A shows an example of the TMS during an entire recording session, with epochs of stimulus presentation denoted by green horizontal bars. As illustrated in this representative example, the TMS increases during virtually every period of stimulus presentation. To analyze the average effect of stimulus presentation on the TMS, we derived the stimulus-triggered TMS. For this analysis, we focus on one stimulus set, which includes three stimuli, each at three different dilutions allowing us to relate the TMS to stimulus intensity (Table 1, set D2). The global stimulus-triggered TMS, based on all trials from all sessions is shown in Figure 7B. The stimulus-triggered TMS is tightly linked to stimulus presentation, first increasing during stimulus application, and further rising following electrical vomeronasal organ stimulation. We note that this temporal profile is consistent with the observation that AOB neurons sometimes respond after stimulus application, before electrical activation of the sympathetic nerve trunk (Bergan et al., 2014; Bansal et al., 2021). Note also that the TMS (Fig. 7B) is essentially zero during electric stimulation (Fig. 7B, two red bars). This is because stimulation artifacts assume a distinct energy distribution over channels compared with that associated with LFP oscillations. A very similar profile is observed in another set of experiments using different stimuli (Extended Data Fig. 7-1A). In summary, the TMS and hence the oscillating network state are closely correlated with chemosensory stimulus processing in the AOB.
Both the TMS and the LFP envelope depend on stimulus identity
We have shown above that the TMS rises during stimulus presentation. But is the TMS merely a reflection of stimulus uptake, or does it instead depend on sensory stimulus features? The fact that TMS magnitude varies between trials (Fig. 7A) already suggests that the network pattern reflected by it might depend on the sensory stimulus. To address this hypothesis, we derived the stimulus-triggered TMS separately for each stimulus (and averaged it across sessions with the same stimulus set). Comparison of stimulus-triggered TMS values reveals a stimulus-specific and dose-dependent response, with more concentrated stimuli eliciting higher stimulus-triggered TMS profiles (Fig. 7C). Larger responses to female over male stimuli are also consistent with previous reports by us and others (Hendrickson et al., 2008; Ben-Shaul et al., 2010; Bansal et al., 2021). Statistical comparison of TMS distributions associated with each of the stimuli reveals a strong effect of the stimulus with widespread pairwise differences among stimuli (Extended Data Fig. 7-1).
Finally, recall that the TMS is merely a measure of similarity to the oscillating state, and, thus, high TMS values do not necessarily reflect high LFP intensities. To test whether LFP intensity is also stimulus dependent, we derived the average LFP envelope. Here, we consider time points during which the TMS is above threshold (Fig. 3E). Mean LFP envelopes for each of the nine stimuli (during one session) are plotted according to their spatial location in Figure 7D. These plots indicate that, like the TMS, higher LFP magnitudes (envelope) are associated with presentation of specific (and often more concentrated) stimuli. The mean LFP envelopes (averaged across the 32 channels) for all five sessions from this dataset are shown in Figure 7E. These examples demonstrate that, in all sessions, LFP magnitude varies with stimulus identity and, in many sessions, also with stimulus dilution. Together, our results indicate that both the spatial pattern of the LFP (reflected by the TMS), and its intensity, are stimulus dependent.
Relationship between the template match signature and single-unit firing rates
Because AOB single-unit activity typically increases following stimulus presentation (Luo et al., 2003; Ben-Shaul et al., 2010; Arnson and Holy, 2013; Tolokh et al., 2013; Doyle et al., 2016; Bansal et al., 2021), it may be expected that the TMS and spiking rates would be correlated. To test whether this is the case, we plotted the TMS and the summed rate (combined across all recorded single-units and multiunits) during an entire session (Fig. 8A). Indeed, they appear correlated at this coarse timescale. Yet, finer temporal resolution reveals that the overlap is not precise (Fig. 8A, bottom).
To reveal the temporal relationship between the two signals, we calculated the cross-covariance between the firing rate of each single-unit (n = 520) and the TMS within the same session, and noted the lag at which the maximum cross-covariance is obtained. A positive lag indicates that the TMS precedes the single-unit rate vector. As shown in Figure 8B, most, but not all, of the units are positively correlated with the TMS. While lags associated with high cross-covariance values may be either positive or negative, the highest correlations occur for smaller lags. Across all units, the median cross-covariance is 0.067 (n = 520; signed-rank p-value = 1.03 · 10−42), while the median lag is 0.247 s (signed rank p-value = 0.002).
To explore the relationship between these signals in the context of sensory processing, we superimposed the stimulus-triggered TMS and the summed single-unit peristimulus time histogram (PSTH), for each session. Examples for six sessions are shown in Figure 8C (Extended Data Fig. 8-1, more examples). Note that in many cases, both signals peak together, suggesting that they are physiologically related. Notably, the average TMS generally provides a more reliable indication of sensory stimulation than the combined spike rate. This is likely because single-unit sampling is sparse and random (with marked differences in the yield of stimulus-responsive neurons across sessions), while the LFP is inherently a more robust population level signal. These observations not only highlight the relevance of the TMS for sensory processing, but also as a reliable (yet low-dimensional) indicator of stimulus-induced vomeronasal system activity.
Oscillations are most prominent in the ECL and likely reflect a reciprocal exchange between AMCs and granule cells
Having shown the relationship between oscillatory events and single-unit activity, and their dependence on stimulus presentation, we turn to investigate their network source. Most of our recordings were conducted with the planar electrode array advanced in parallel to the AOB ECL (Fig. 9A). The ECL contains somata of AMCs and displays characteristic signatures of spiking activity (Luo et al., 2003; Ben-Shaul et al., 2010). In a subset of experiments (Table 1), we advanced electrodes perpendicular to the AOB, in an attempt to span its layers (Fig. 10A). At a coarse scale, these recordings support the notion that the origin of the oscillation is indeed the AOB, rather than volume transmission from another region. Specifically, in all planes defined by the long axes of the electrode array (Figs. 9A, 10A, Extended Data Fig. 2-1) and along the mediolateral axes (seen in parallel penetrations), maximal LFP intensities typically appear within the borders of the electrode array. We next consider the layer-specific localization of these oscillations within the AOB. Although we dip electrodes in a fluorescent dye to routinely confirm within-ECL location post hoc, reconstructions are not accurate enough to assign individual contact sites to specific layers. Yet, we gain important information about electrode localization by comparing hotspots of single-unit neuronal activity with LFP energy distributions (Fig. 9, compare B, E, and compare F, I). Comparison of LFP signatures and single-unit activity in individual sessions reveals a clear correspondence between oscillatory LFP power and neuronal firing activity, which is often associated with stimulus-induced neuronal responses in the ECL (data not shown).
Next, we reasoned that if LFP waveforms measured during parallel penetrations originate from the same neuronal structure, the oscillations are likely to be in-phase across channels, without reversal of polarity. Figure 9, C and G, shows all LFP channels during two oscillatory episodes, with expanded temporal views in Figure 9, D and H. These representative examples indicate that all channels within a session are indeed in-phase, suggesting that they either represent a single common process or, alternatively, multiple coordinated processes.
Continuing this line of reasoning, if oscillations result from an exchange between neuronal elements in different AOB layers, their polarity should reverse at the boundary between these layers. Indeed, in most perpendicular penetrations (Fig. 10A), we observe a clear LFP reversal between deeper and superficial sites. Examples from three different sessions are shown in Figure 10B–D. The simplest interpretation of this reversal is that it represents the transition along the lateral olfactory tract, which separates the ECL from the AOB granule cell layer. Consistent with the observation that the AOB granule cell layer activity is seldom detected by extracellular recordings, we observe neuronal spiking activity above, but not below the reversal (Fig. 10, compare B, left, right, C, left, right, D, left, right).
Further support for the hypothesis that LFP reversals represent the transition between ECL and granule cell layer is gained from current source density analysis. Current source density is the second derivative of the electric potential with respect to position, and provides a measure of current sources and sinks (Mitzdorf, 1985). We calculated it along the length of the recording probes that contain equally spaced recordings sites (and ignore the spread between adjacent electrodes). Figure 10E shows the current source density (represented by colors) along the sites of LFP polarity reversals (Fig. 10B–D, shaded regions) superimposed on the original LFP traces. These plots reveal an exchange of sinks and sources within each cycle period. Importantly, the current source density itself is reversed, matching the LFP reversal.
As a further test that the oscillations along the extent of the AOB emerge from a single source, we applied independent component analysis to individual oscillatory episodes. This analysis indicates that, for both parallel (Fig. 9J,K) and perpendicular penetrations (Fig. 10F–H), oscillations along the extent of the AOB are associated with a single dominant factor.
Integrating our observations, we propose the following scenario: vomeronasal organ activation elicits spiking activity of single AOB neurons. With sufficient single-unit activation, the network begins to recruit granule cells, which supply inhibitory feedback to AMCs and entrain oscillations. These oscillations then persist as long as incoming inputs provide sufficient drive to AMCs. Once AMC firing rates drop (because of reduction of incoming inputs), reduced granule cell activity ceases to drive oscillations. Under this scenario, full-blown LFP oscillations represent the peak of stimulus-induced AOB activity. This is consistent with the observation that oscillatory episodes (Fig. 4E) are briefer than stimulus-evoked responses of many individual AOB single-units (Yoles-Frenkel et al., 2018).
Firing synchrony among simultaneously recorded AOB neurons
Finally, we sought to evaluate the temporal scales of coordinated activity in our dataset. To this end, we tested whether spikes of simultaneously recorded neurons tend to co-occur during oscillatory periods. More specifically, we divided each oscillatory episode into temporal bins of varying durations and asked whether spikes from different neurons “cluster” in specific windows more than expected by chance (Fig. 11A). Chance levels are obtained by randomly shifting the spike times of each individual neuron (see Materials and Methods). Among other factors, the sensitivity of this analysis varies as a function of episodes per session (in our data: median, 41.5; minimum, 3; maximum, 395) and the number of active single-units recorded during these episodes (range, 2–42 units). Furthermore, we note that a fixed bin size is not optimal for detecting synchronization because of frequency modulated oscillations. Keeping these considerations in mind, we observe that window sizes from 100 ms to 2 s are all associated with a modest, yet highly significant overabundance of temporal coordination (Fig. 11B). Note that while synchronization at larger windows (i.e., larger than a typical 200 ms theta cycle) likely reflects the slower temporal dynamics of the rate-modulated neuronal response (Yoles-Frenkel et al., 2018), the significant synchronization at 100 ms windows (p = 0.005) is consistent with temporal clustering within a theta cycle. As noted above, our analysis is limited by the number of simultaneously recorded neurons, and this factor is more critical with higher temporal resolution. Thus, we cannot rule out the existence of synchronization at finer temporal windows, <100 ms).
The effects of phase-dependent firing markedly increase with neuronal population size
A key hypothesis that emerges from our work is that neuronal oscillations play a role in information processing by modulating the timing of AOB neurons. For that to occur, the degree of phase locking observed here (Fig. 5) should suffice to induce a substantial effect on downstream neurons. We now present a simple calculation to demonstrate that even modest deviations from uniformity can add up and become highly meaningful when considering coactive neuronal populations. Adopting the perspective of a downstream neuron receiving a fixed number of inputs (Fig. 12A), we calculate the probability of obtaining a particular number of spikes within a specific integration window. We compare this probability between a uniform and a nonuniform distribution of firing within a cycle, keeping the overall firing rates constant. For example, we can compare neurons that fire at a constant rate of 5 Hz throughout the cycle, to neurons that fire at 4.5 Hz within one-half of a cycle, and at 5.5 Hz during the other half (i.e., during the negative LFP phase). We then calculate the probability of obtaining a particular number of spikes during the negative phase of the LFP. Examples of such phase dependencies and their respective phase-locking values are shown in Figure 12B.
Effects of phase-dependent firing on the probability of coincident spikes. A, Schematic of the analysis. B, Examples of nonuniform spike time distributions within a cycle with simulated phase-locking values. Compare with real data in Figure 5A. C, Analysis of nonuniformity corresponding to the cases shown in B. Each panel shows the probability to obtain a certain number of spikes (or more), within a single integration window of 100 ms (i.e., half a 5 Hz theta cycle). Each row represents a different number of neurons (from top to bottom: 10, 50, 100, and 200), and each column represents a different degree of spike time locking, quantified by the firing rate within that period. In all of these examples, the baseline rate (without spike time locking) is 5 Hz, which corresponds to a probability of firing of 0.005 within a 1 ms window. In all plots, dashed blue traces indicate the probability under uniform firing, while dashed green traces indicate the probability with phase-dependent firing. Red traces represent the difference between these cumulative probabilities. Note the changes in horizontal scales between the different rows.
Specifically, we ask how the probability of obtaining a given number of action potentials within a window depends on the firing rate of each neuron within the window and the population size. This probability can be calculated using the binomial cumulative density function. As shown in Figure 12C, the differences, which increase with population size and the extent of phase locking, can be significant. For example, consider a population of 50 neurons, an integration window of 100 ms (half of a 5 Hz theta cycle), and a baseline probability of 5 Hz (uniform firing), compared with a mildly locked nonuniform distribution (50% increased probability = 7.5 Hz within the window). Under these conditions, the probability of observing 30 spikes within that window is 0.14 for the uniform distribution, and 0.88 for the phase-locked distribution. Thus, in this example, a moderate degree of phase locking can increase the probability of observing 30 spikes within a 100 ms window sixfold. Stronger effects are attained with tighter phase locking and larger populations. In conclusion, a modest degree of phase locking, as observed here, can exert dramatic effects on downstream neurons with sufficiently large populations.
Discussion
In this manuscript, we studied neuronal mechanisms underlying processing of social information in mice. Specifically, we focused on oscillatory LFP patterns in the AOB. While spiking activity is arguably the most informative measure of neuronal communication, LFPs can reveal population-level processes not easily observed by the activity of (even many) single-units. We provide evidence that LFP oscillations are generated locally within the AOB, are associated with stimulus sampling, and occur even in the absence of arousal or active investigation. We have shown that they fall within the theta band (2–12 Hz), but their sporadic appearance indicates that they are distinct from breathing-associated theta-band main olfactory bulb oscillations. Most importantly, we demonstrated that these oscillations exert a marked influence on neuronal spike timing, which can profoundly influence downstream processing. We thus propose that these oscillations shape AOB outputs as they are relayed to pathways that control behavioral and physiological function.
The template-matching approach
We discovered sporadic theta oscillations by manually browsing LFP data. To systematically analyze their relationship to stimulus presentation and unit activity, we relied on their highly conserved envelope distributions. Using the mean energy distribution as a template, we scanned entire recording sessions, thereby generating a continuous TMS. We found that the TMS fluctuates as a function of sensory stimulation, with peaks representing full-blown oscillations following efficient recruitment of the AOB network. The LFP reflects the activity of many neurons and is considerably easier to record than single-unit activity, and, thus, a small number of electrodes sampling LFP activity can provide a reliable readout of sensory activation. Therefore, practically, the TMS provides a reliable way to monitor stimulus-induced activity in the AOB and, potentially, other regions as well.
Source of the oscillations
The fact that oscillatory amplitudes often attain maxima within the central sites of the electrode array (in both parallel and perpendicular penetrations) indicates that oscillations are generated within the AOB. The high temporal correlation between signals across channels is consistent either with a “point source” that spreads in space, or a more distributed spread of coordinated generators. Independent component analysis (Figs. 9, 10) also supports the presence of a single source that underlies the oscillations across channels. We cannot determine which of these scenarios is more accurate, but the nonisotropic spread of the LFP envelope (Extended Data Fig. 2-1) favors the latter possibility. Important insights can be gained from previous work in AOB slices. Specifically, Jia et al. (1999) recorded LFP signals in different AOB layers during vomeronasal nerve activation. Based on their analyses and our own observations, we propose that the oscillations observed here result from ongoing AOB activation, which, if sufficiently strong, leads to reciprocal exchange between AMCs and granule cells. We speculate that stronger sensory activation not only increases the rate of AMC activation, but also that of granule cells, resulting in increased inhibitory drive onto AMCs and hence higher oscillatory frequencies. As the incoming signal wanes, there is no sufficient drive to activate AMCs and granule cells in synchrony, and oscillations decelerate and eventually disappear. Thus, frequency changes within a single oscillatory episode may result from dynamic changes in AMC inputs. It is important to note that, while these AOB oscillations fall within the theta band, they are mechanistically distinct from main olfactory bulb theta oscillations (Kay, 2014). Since all our experiments were conducted in anesthetized mice, we conclude that the oscillations are triggered by sensory activation rather than active exploration or top-down feedback from other areas. While most of our experiments included tracheotomy, some did not (in these experiments the nasal cavity was not flushed). Notably, also in these experiments, we did not observe breathing-associated theta oscillations, further supporting the notion that breathing does not play a direct role in the oscillations studied here.
Previous studies on AOB theta oscillations
Our study follows previous reports of theta-band AOB LFP oscillations. In addition to the in vitro analysis mentioned above (Jia et al., 1999), work from the Brennan laboratory (Binns and Brennan, 2005; Leszkowicz et al., 2012) and Lanuza laboratory (Pardo-Bellver et al., 2017) revealed stimulus-induced AOB theta-band oscillations. The former study also showed that oscillatory LFP bands can change following mating. Recordings in rats (Tendler and Wagner, 2015) revealed that coherence among theta LFP in various limbic regions (including AOB) depends on behavioral context. The widespread distribution of theta oscillations throughout limbic regions suggests that the patterns observed here may induce or resonate with downstream vomeronasal system regions. While our recordings indicate that oscillations can arise in the absence of arousal, the work of other groups highlights the relevance of oscillations in the freely moving state. Another prominent class of AOB oscillations, described by us and others, are infraslow oscillations (Gorin et al., 2016; Zylbertal et al., 2017a,b; Tsitoura et al., 2020). These are much slower than those observed here, likely play a different role, and arise from distinct mechanisms. Nevertheless, there may be a direct relationship between neuronal ensembles such as identified in Tsitoura et al. (2020) and those proposed here (see below; Fig. 13).
Models of physiological roles of theta oscillations and the effect of phase concentration. Cartoons of three different scenarios for phase-dependent spike timing. Each scenario depicts spike trains from four neurons, and several LFP oscillation periods. Stimulus presentation is indicated by the colored horizontal bars. Shaded regions behind neuron icons represent coactive ensembles. A, Uniform locking of spikes from all neurons regardless of the stimulus. Here, oscillations can amplify and thus gate the influence on downstream processing stages, in a manner that is common for all units and sensory stimuli. B, Locking of spikes, so that each unit has a preferred phase within the broad region of preference. Here, there may be distinct ensembles that represent distinct stimuli. C, Stimulus- and unit-dependent phase locking. In this scheme, an elaboration of B, ensemble memberships may change as a function of the stimulus.
Mechanisms underlying generation of the oscillations
Mechanistically, we proposed that AOB theta oscillations result from an interaction between AMCs and granule cells. This conclusion is based on spatial analysis of LFP waveforms (Figs. 9, 10) and previous work in slices (Jia et al., 1999). Notably, however, in the main olfactory bulb these circuit elements were proposed to underlie much faster gamma rhythms (Fukunaga et al., 2014; David et al., 2015; Li and Cleland, 2017; Kersen et al., 2022). This apparent discrepancy may be reconciled by the fact that AOB and main olfactory bulb neurons differ in intrinsic properties (Zibman et al., 2011; Shpak et al., 2012), connectivity (Larriva-Sahd, 2008), and the inputs they receive. All of these features (and others) can markedly affect oscillation frequencies (Li and Cleland, 2017; Burton and Urban, 2021; Kersen et al., 2022) and potentially account for the phenomenology of the oscillations observed here. Further experimental and modeling studies of the AOB network are required to pinpoint the mechanisms that produce AOB oscillations.
The role of active sampling and its influence on sensory processing is widely discussed in the context of the main olfactory system, where sampling patterns can modulate the activity of sensory neurons, and hence of downstream targets (Kepecs et al., 2007; Verhagen et al., 2007; Wachowiak, 2011; Lefèvre et al., 2016; Jordan et al., 2018; Ackels et al., 2020). This suggests another source that could modulate the intensity or pattern of the oscillations, namely, the vomeronasal pump. While there is evidence for intrinsic rhythmicity of the pump (Meredith and O'Connell, 1979; Meredith, 1994), given its slow dynamics, it is hard to imagine that these AOB oscillations directly reflect pumping rhythmicity. Furthermore, AOB oscillations are stimulus dependent, and pumping, at least under anesthesia, is unlikely to be. Nevertheless, the intensity of investigation may influence pump dynamics, and hence alter sensory neuron and AMC activity. To reveal the influence of active sampling on these oscillations, it will be necessary to record pumping activity (and neuronal activity) in awake behaving animals, a goal that is presently not within reach. Alternatively, vomeronasal organ sampling may not be as amenable to controlled voluntary modulation as is sniffing, raising the idea that the AOB realizes an intrinsic oscillatory mechanism serving a similar functional role (i.e., gating, binding, and multiplexing). Along the same lines, inputs from other brain regions could modify AOB neuron excitability and thus modulate the presence and expression of oscillations. Investigating this possibility will require enhancement or suppression of feedback inputs from other brain regions (Fan and Luo, 2009; Boyd et al., 2012; Markopoulos et al., 2012; Otazu et al., 2015; Cádiz-Moretti et al., 2016).
The physiological relevance of AOB theta oscillations
The existence and prevalence of the stimulus-induced theta oscillations beg the question whether they play a role in information processing. We suggest that they are not an epiphenomenon, but rather an important element of AOB physiology. Below we propose several models, all of which are inspired by previous work in the main olfactory system (Chaput, 1986; Schoppa, 2006; Cenier et al., 2009; Cury and Uchida, 2010; Smear et al., 2011; Gschwend et al., 2012; Uchida et al., 2014; David et al., 2015; Kay, 2015; Frederick et al., 2016; Osinski and Kay, 2016; Zelano et al., 2016; Jiang et al., 2017; Fourcaud-Trocmé et al., 2018; Losacco et al., 2020; Burton and Urban, 2021) and other regions (Buzsáki, 2006; Colgin, 2013).
One prominent aspect of the vomeronasal system and the AOB is its integrative structure. Along with established patterns of labeled line coding in the AOB (Ferrero et al., 2013; Ishii et al., 2017; Osakada et al., 2018), many natural stimuli are complex blends of molecules (Ben-Shaul et al., 2010; Hammen et al., 2014; Kahan and Ben-Shaul, 2016), implying that stimulus recognition may involve integration of information across channels (Kaur et al., 2014). Indeed, AOB circuitry itself already realizes substantial integration, as AMCs often sample information from multiple glomeruli (Takami and Graziadei, 1991; Wagner et al., 2006; Larriva-Sahd, 2008; Mohrhardt et al., 2018). If distinct channels convey related information, there may be an advantage for their coordinated activation. Thus, one simple hypothesis (Fig. 13A) is that the oscillations serve to constrain neuronal firing to limited temporal windows, and thus activate downstream regions more efficiently (Fig. 12). In this scenario, oscillations serve a gating and amplification function: low neuronal activity levels will not give rise to oscillations, but once the activity is sufficiently high, oscillations will appear, constrain spike timing to smaller windows, and thereby facilitate and gate information transfer to downstream, regions.
A more elaborate hypothesis suggests the existence of multiple ensembles, so that members of each are activated together. Under this scenario (Fig. 13B), neurons within an ensemble share a similar phase preference and are thus more likely to exert a synergistic effect on downstream targets. The notion that different units follow distinct phase preferences is supported by the analysis shown in Figure 6. However, unlike the main olfactory system, where spike times of neurons “tile” the theta breathing cycle (Cury and Uchida, 2010; Shusterman et al., 2011), here most of the activity is concentrated within about a third of a cycle. Finally, the most elaborate scenario holds that phase distributions of units are stimulus dependent (Fig. 13C). According to this model, each neuron can participate in different ensembles depending on the stimulus. Support for this idea comes from the analyses shown in Figure 6, B and C. The two latter models may seem inconsistent with the idea, promoted also by us (Yoles-Frenkel et al., 2018) and others (Zibman et al., 2011), that the vomeronasal system is slow and temporally imprecise. However, there is no real discrepancy, since the synchronization proposed here involves joint activity of neurons, rather than locking to stimulus uptake.
Potential effects of coordinated activity on downstream processing stages
It should be noted that the term gating may refer to two distinct aspects. The first is that oscillations influence spike timing of AMCs, and thus, at least to a degree, gate their activity. The second interpretation is that since spike times of AMCs are clustered in time, they can more efficiently activate downstream neurons, and thus gate their activation. Definitive proof of either of these interpretations would require intracellular recordings of either AMCs or downstream regions. However, our emphasis is on the second interpretation, namely gating of activation of downstream neurons, and our analyses support the feasibility of this scenario. Specifically, as we have shown above (Fig. 12), even a modest degree of phase-dependent firing can exert a profound effect on joint neuronal firing probability within a limited temporal window. Furthermore, we have provided evidence for coordination at a temporal scale finer than a theta cycle (i.e., 100 ms; Fig. 11). It is notable that we were able to detect such coordination despite the modest number of simultaneously recorded neurons, and despite the fact that we have used fixed bin widths to detect synchronization associated with frequency-modulated oscillations. Parallel recordings of a larger number of units may reveal that timescales of coordination are even finer than observed here. Ultimately, the effects of timing on information processing depend on the integrative properties of downstream processing stages, whose details are not precisely known. This includes the degree of anatomic convergence, and the spatial and temporal integrative synaptic properties of neurons in these regions. Presently, our findings highlight a potential mechanism in which coordination of spike timing may not only improve the efficacy of activation, but may also organize activity from distinct neuronal subsets. Understanding which, if any, of these mechanisms play a role in AOB information readout, remains a central goal for future research.
Footnotes
This work was supported by Deutsche Forschungsgemeinschaft (German Research Foundation) Grant 378028035 (to Y.B.-S. and M.S.), German-Israeli Foundation for Scientific Research and Development Grant 1-1193-153.13/2012 (to Y.B.-S. and M.S.), and Israeli Science Foundation Grant 1703/16 (to Y.B.-S.). A.K. was supported by Lady Davis postdoctoral fellowship. S.T.M was supported by a Minerva Stiftung fellowship for this research.
The authors declare no competing financial interests.
- Correspondence should be addressed to Yoram Ben-Shaul at yoramb{at}ekmd.huji.ac.il
