Natural Statistics as Inference Principles of Auditory Tuning in Biological and Artificial Midbrain Networks

Animals

Big brown bats (Eptesicus fuscus) were collected from an exclusion site under a state permit. All experimental procedures were conducted in accordance with a protocol approved by an Institutional Animal Care and Use Committee. A total of ∼100 bats were housed in our laboratory and used for vocal data recordings, and four (two male, two female) bats were used for neurophysiological data collection.

Audio recordings for training the biomimetic network

A bat call library was built from audio recordings of bats housed in a vivarium room where the temperature is kept at 70–80°F, and humidity is kept at 30–70%. This room holds ∼100 bats in groups of one to six separated in mesh cages. The recordings were made for 2 d using an Avisoft CM16/CMPA ultrasonic microphone and the Avisoft-RECORDER software. Mono audio was recorded at a sampling rate of 300 kHz.

Natural call recordings from big brown bats were processed to extract meaningful segments. An energy-based signal activity detection was performed on the entire database to remove the silences between calls and to split the recordings into segments containing bat calls (Park et al., 2014). As a result, we constructed species-specific databases containing 17,713 calls (∼10 min) for big brown bats. This call database was used for training artificial networks. The data were divided into a training set (15,000 randomly selected calls) to learn network parameters and test set (remaining 2713 calls) for verifying the network.

Social calls for natural sound representation

To investigate discriminability in the artificial network, we used a social call database that includes 26 audio clips for eight different types of bat calls (Fig. 2). These types include six calls, as defined in (Wright et al., 2013), specifically, Echo, frequency-modulated bout (FMB), upward frequency modulated (UFM), long frequency modulated (LFM), short frequency modulated (SFM), and chevron-shaped (CS); in addition to two additional calls types, long-wave and hook, which resemble a hook in time-frequency space. All audio clips were up-sampled from 250 to 300 kHz.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Example spectrograms of eight bat calls in social call database. A, Echo. B, FMB. C, UFM. D, LFM. E, SFM. F, CS. G, Hook. H, Long-wave.

Receptive field recordings

A head-post was adhered to the skull of bats for head fixation as described previously (Macías et al., 2018). The IC was located using skull and brain landmarks and a surgical drill was used to make a ≤1-mm diameter craniotomy preserving dura. The neurophysiological recordings were performed in a sound-attenuating and electrically shielded chamber (Industrial Acoustics Company). Each bat was restrained individually in a custom-made foam mold and the head was fixed by the head-post. Recording sessions were conducted over three to five consecutive days, each one lasting no more than 4 h. Water was offered to the bats every 2 h. No drugs were administered during recordings. During recordings a silver wire for grounding was placed in between muscle and skull ∼5 mm rostral to the craniotomy site. The 16-channel recording probe (Neuronexus A1x16-5 mm-50-177-A16) was inserted into the brain using a micromanipulator. The surface of the brain was registered as 0 μm for depth reference and the probe was advanced in 10 μm steps using a hydraulic microdrive (Stoelting Co). Recordings were taken at least 100 μm apart. An OmniPlex D. Neural Data Acquisition System recording system (Plexon) was used to obtain neural responses with 16-bit precision and 40-kHz sampling rate. A transistor-transistor-logic (TTL) pulse for each stimulus presentation was generated with the National Instrument card used for stimulus presentation and was recorded on channel 17 of the analog channels of the acquisition system for synchronization of acoustic stimuli and neural recordings. The stimuli were recorded on channel 18 of the acquisition system to corroborate synchronization.

Moving ripple stimuli

A set of ripple stimuli was generated to estimate STRFs of IC neurons (Kowalski et al., 1996; Depireux et al., 2001; Andoni et al., 2007). Ripples are modulated noise stimuli that are dynamic both in time and frequency. Each ripple can be described as S(t,x)=1 + ΔA×sin(2π(ωt + Ωx) + ϕ), (1)where t and x are indices for time and octave scaled frequency. ΔA and ϕ are amplitude and a phase, respectively. And ω and Ω represent modulation rates along temporal (Hz) and spectral (cyc/oct) axes. The temporal and spectral modulation parameters were varied from −176 to 176 Hz in steps of 32 Hz and 0.0–1.5 cyc/oct in steps of 0.15 cyc/oct spectrally (Fig. 3A). Each ripple spanned 6.66 octaves from 1.2 to 121 kHz and was 300 ms in duration.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

RTF extraction. A, A subset of ripple stimuli. B, Input patch sequence configuration from a ripple stimulus to a code vector for characterizing the network. C, Examples of responses on each node (each element of the code vector), and definition of magnitude m and phase ϕ in a ripple response. D, Magnitude and phase plots on one of nodes.

Audio playbacks for neural recordings

Extracellular recordings from the IC of awake animals were taken while they passively listened to broadcast of either ripple stimuli, or pure tones at 70 dB. All stimuli were generated at a sampling rate of 250 kHz using a National Instruments card (PXIe 6358) and transmitted with a calibrated custom-made electrostatic ultrasonic loudspeaker connected to an audio amplifier (Krohn-Hite 7500). The loudspeaker was placed at 60 cm (for all ripple and pure tones stimuli) from the bat’s ear. The frequency response of the loudspeaker was compensated by digitally filtering the playback stimuli with the inverse impulse response of the system as described previouosly (Luo and Moss, 2017).

Frequency tuning curves were built by recording neural responses to pure tones of 5-ms duration (with 0.5-ms ramping rise and fall). The tones ranged between 20 and 90 kHz (in 5 kHz steps) and the sound pressure levels ranged from 20 to 70 SPL (10 dB steps). At each recording site first, we played 20 repetitions of the randomized ripple stimulus and then 15 repetitions of each of the randomized pure tones at a different SPL.

Analysis of neuronal responses

For the analysis of auditory tuning in response to ripple and pure tone stimuli, responses were sorted offline, then single units were detected using the program Wave_clus (Quiroga et al., 2004). Each individual waveform was inspected and the acceptance threshold for clusters was <10% of spikes with <3-ms interspike interval, consistent with the neuronal refractory period. Any sites that showed no response to ripple stimuli were excluded from the spike sorting and further analysis in line with procedures used in other studies (Poon and Yu, 2000; Escabi and Schreiner, 2002; Andoni et al., 2007). After spike sorting, the Euclidian distance error between the mean and variance of number of spikes across trials was computed. Units whose error is <1.0 were selected for further analysis, following a Poisson model of spike representation (Corrado et al., 2005; Schwartz et al., 2006). This analysis resulted in 108 single units used for the current study.

Neural discriminability of conspecific calls

In order to examine selectivity of IC neurons to calls from the bat’s natural repertoire, we re-used neural data previously collected in an earlier study (Salles et al., 2020), where we collected neuronal responses to Echo calls versus FMB social calls. The study followed the same methodology for data collection as described here. Wave_clus was used to detect and classify single units from the recordings. The spikes responding to either FMB or Echo were counted in windows of 25 ms in duration, starting 5 ms after stimulus onset. Some units with an average of less than five spikes over 20 times recordings were excluded because they were considered as a non-responsive unit to the stimulus. Multiunit activity was determined from interspike intervals with <3 ms that were inconsistent with neuronal refractory period; and units with >10% of spikes with <3-ms interspike interval were excluded from analysis. As a result, total 575 units were finally obtained and their responses are used in the present work to contrast neural discriminability between Echo and FMB calls with artificial neurons.

Artificial network front-end processing

To develop a biomimetic architecture, a biologically-inspired auditory spectrogram is used as input for the network (Shamma, 1985a,b; Yang et al., 1992; Wang and Shamma, 1994). The auditory spectrogram incorporates four processing stages that emulate peripheral processing in the mammalian system: cochlear filtering, auditory-nerve transduction, hair cell responses, and lateral inhibition (Chi et al., 2005). Briefly, an incoming acoustic waveform is analyzed along a bank of constant-Q filters spanning a logarithmic scale. Then, each frequency channel undergoes a high-pass, nonlinear compression and low-pass filtering followed by lateral inhibition across frequency, following the implementation available in the NSL toolbox (Chi and Shamma, 2005) with the following settings: the frame length was set to 0.2 ms without overlap, and each octave was represented with 24 channels (i.e., 128 channels over 5.33 octaves). Octave-scaled center frequencies were represented as fc=440×2((c−32)/24+γ) where fc is a center frequency of the cth channel, and γ is a constant factor of octave shift (γ=4.38 ). Inputs to the artificial network were sampled as square patches of the spectrogram spanning 128 frequency channels (i.e., 5.33 octaves) and 160 time samples (i.e., 32 ms).

Structure of artificial network

An artificial neuron, i.e., node mimicking a biological neuron, is mathematically modelled by a linear combination of prenode outputs and a nonlinear activation function. An artificial network is constructed by connecting a large number of nodes to each other. Using nonlinear activation functions enables the network to perform nonlinear computations on feedforward propagation. For this study, we favored a generative architecture using an autoencoder composed of an encoder, which compresses original data into a compact code, and a decoder, which reconstructs the original signal from that code (Baldi, 2012; Doersch, 2016). The intuition is to directly test our hypothesis that the network would infer a statistical model of the training dataset of natural calls, and if successful should allow a faithful reconstruction of the inputs.

The proposed architecture is shown in Figure 4. First an encoder stage E is composed of convolutional layers, pooling layers, and a fully connected layer. A latent vector represents compressed features learned from the input data. A decoder stage D composed of reverse operations using transposed convolutions, reconstructs the input features from a latent vector. A sampling stage, interposed between the encoder and decoder, emulates neural activity yielding sparse binary activations.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Convolutional layered autoencoder structure for biomimetic network. The flow denoted in black shows a double stacking structure as a standard example. Based on this structure, a deeper structure can be constructed by stacking more modules, on the other hand, a shallow structure is created by removing a module on the top of the standard example.

Using the same general building block composed of convolution and pooling layers, this study investigates various configurations of the network by varying (1) depth, which is the number of blocks (in Fig. 4, the black-flow shows a double stacking structure as an example); a deeper network can be constructed by stacking more blocks, on the other hand, a shallow network can be created by removing a block; (2) nonlinearity, by varying the slope of nonlinear activation function employed; and (3) sparsity, by controlling the density of sampling in the latent space.

The encoder architecture E follows a convolutional neural network (CNN) framework to reduce the number of trainable parameters, hence controlling for overfitting issues and generalizability to unseen data (Dietterich, 1995). The convolutional layers compute output feature maps using 2D convolutions between input feature maps and several filters as follows: Iol[f,t,k]=∑ξ,τ,mIil[ξ,τ,m]fl[ξ−f,τ−t,m,k], (4)where f,t,l,k and m are indices for spectral, temporal, layer, channel of output feature map, and channel of input feature map, respectively. Ii , Io , and fl are feature maps for input and output, and convolutional filter applied in the lth layer, respectively. Multiscale filters are employed in each convolutional layer to balance broad span (in time and frequency) versus localized analyses. Then, output feature maps concatenate filter outputs using multisized filters (Fig. 5A; Szegedy et al., 2015). Specifics of both filter composition and dimensions of intermediate feature maps are summarized in Table 1. Neural activation by an acoustic feature is emulated by applying a nonlinear function after convolution as follows: Ial[f,t,k]=max(Iol[f,t,k],α×Iol[f,t,k]), (5)where α is a constant within an interval [0, 1] (Maas et al., 2013). Next, pooling layers compress the output from the previous convolutional layer by extracting a maximum among some values enclosed by a non-overlapping window (i.e., max pooling) Il (Scherer et al., 2010). As a result, the width and height of the output are reduced by half. At the top of the encoder, a fully connected layer is applied for mapping into a latent space, which involves natural statistics requiring to reconstruct original input, as vc=WlT×flatten(IL) where IL is a feature map in the last pooling layer, W is weight matrix in the fully connected layer, and flatten(.) is a reshape function from a 3D tensor to a vector.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Operations using multiscale filters. A, Convolution using multiscale filters. B, Transposed convolution using multiscale filters.

In the middle stage, a binary code vector vb is generated by performing a Bernoulli sampling process. A sigmoid function is applied to the latent vector to calculate prior probabilities. Thus, the output of the middle stage is represented as vb=Bernoulli(σ(vc)) where σ(.) is a sigmoid function.

The decoder D is composed of a fully connected layer and transposed convolution layers. In the fully connected layer, a latent vector is expanded into an initial space as v̂=Wl×vb , and the vector v̂ is reshaped to a 3D tensor as a set of initial feature maps as Îl=reshape(v̂) . From initial feature maps, a transposed convolution using multiscale filters is sequentially performed until the output has the same dimensions as the input patch (Radford et al., 2015; Shelhamer et al., 2017). Convolutional filters used in the encoder are applied for transposed convolution after transposing input channel from output channel dimension as f˙l[f,t,k,m] . A transposed convolution using multiscale filters is performed in three steps (Fig. 5B). First, the input feature map Îl is split into submaps, [Î1l,Î2l,...,ÎNl] , as many as the number of filters. Second, transposed convolution is individually performed for each pair of submap and filter. Finally, a set of output feature maps is obtained by averaging the results of the second step.

Training artificial network

The network was trained using the cost function: L=12∑n[(xn−D(E(xn)))2 + λ(ρ−∑iσ(vci))2], (6)where xn is an input patch with respect to the nth index, E(.) represents an encoder function while D(.) is for a decoder, and ρ means the average number of active nodes. The first term represents the mean square error between an input patch and its reconstruction by the autoencoder. The sparse constraint prevents overfitting as well as emulates sparsity of active neurons in the brain. Let Y be a random variable representing the number of active nodes by the Bernoulli process. Then, the distribution known as the Poisson binomial distribution is denoted as Pr(Y=ρ)=∑A[∏i∈Aσ(vci)∏j∈Ac(1−σ(vcj))], (7)where A is a set whose elements are possible combination for choosing ρ nodes from N nodes. This distribution can be approximated by Binomial(N,μ/N) where μ=∑iσ(vci) (Choi and Xia, 2002). The network training was implemented using TensorFlow (Abadi et al., 2016). AdamOptimization was applied for an optimizer with 1.0e−4 learning rate. And, λ was set to 1.0e−4 . For more details, readers can find the implementation on http://www.github.com/JHU-LCAP/BioSonar-IC-model/.

Comparisons between the statistic of bat calls and artificial neurons were performed to infer the network configuration that best matches the characteristics of IC neurons (as explained next). The best configuration composed of a triple stacking network, a parameter of nonlinearity α=0.2 in Equation 5 and 10% sparsity constraint in Equation 6.

Biomimetic STRFs

Once trained, the network was interrogated following the same procedure as biological neurons. The same ripple stimuli were given as input to the network and activity of the nodes before applying the sigmoid activation and the Bernoulli sampling, vc in Figure 4 was characterized. Each ripple was transformed into an auditory spectrogram (as described earlier). A sequence of input patches for each ripple were then composed by applying a sliding window (window length: 160 frames) in every 2 ms (sliding step: 10 frames; Fig. 3B). Input patches in the sequence were consecutively fed into the pretrained encoder, then a latent vector vc was obtained every 2 ms. The same procedure for extracting biological STRFs was followed (see above, Analysis of neuronal responses). To find the magnitude m and phase ϕ of the responses, we performed a 32-point FFT and derived the magnitude and unwrapped phase of the fundamental component (Fig. 3C). By repeating this procedure for all ripples, the magnitude and phase were collected in a matrix M(Ω,ω) and a Φ(Ω,ω) , respectively (Fig. 3D). These modulation responses were then converted into time-frequency STRF profiles by performing a 2D inverse FFT on the RTF (Fig. 6). Note that, in this study, all network architectures employed a total 100 artificial neurons (spanning a 100D latent space) so that 100-biomimetic STRFs were used for analysis.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

STRF calculation. A, expanded magnitude and phase matrices which are matching to Figure 2D. B, 2D (above) and 3D (bottom) representation of the STRF that is obtained by performing 64 × 64 interpolation and Gaussian smoothing sequentially. In 2D representation, red area represents excitation regions while blue represents inhibition regions.

Natural statistics and Auditory characteristics

Best velocity (BV)

We defined a BV as the center of mass with respect to response power in a magnitude plot. To estimate the center of mass, we performed a Gaussian surface fitting on the first quadrant of magnitude plot. After normalization as M¯e=Me/[∑Ω,ωMeΔΩΔω] where Δω and ΔΩ are, respectively, step size of temporal and spectral modulation rate, the fitting was performed to estimate mean vector and covariance matrix, by minimizing a square mean error function as Err=12∑Ω,ω(ln(M¯e)−ln(Gμ,Σ))2 , where Gμ,Σ is a Gaussian distribution with mean vector μ and covariance matrix Σ. By performing the least square error (LSE) estimator iteratively (Kay, 1993), we derived the Gaussian mean vector and covariance matrix. BV was defined as the slope of the mean vector (Fig. 7).

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

Descriptions for BV and orientation. A, Magnitude plot. B, The result of Gaussian surface fitting. The red ellipse represents Gaussian mean vector μ (center) and covariance matrix Σ (rotation), the BV is defined by a slope of Gaussian mean vector and the orientation error is defined by an angle difference between mean vector and covariance rotation.

A bootstrap for statistical comparison

We performed a bootstrap analysis to evaluate similarity between distributions of characteristics (e.g., FM velocity, BV) comparing natural calls, IC neurons, and artificial neurons. The procedure selects random 30 samples from natural calls in each iteration with replacement. For IC neurons, random samples from each of the four bats are used in each iteration to maintain a balanced representation across bats. In case of artificial neurons, we trained 10 independent networks (using different initialization procedures) and combined the neurons from each network into a complete set that was then sampled during the bootstrap procedure. For each comparison and each bootstrap repetition, the distance between means was noted. A total of 1000 repetitions were used to generate a distribution of mean distances d(μ,σ) where μ and σ are the mean and standard deviation. The p value for accepting null hypothesis was calculated as p=1−2∫0|ε|d(0,σ)(x)dx where d(0,σ) is a zero-mean Gaussian distribution with same variance σ, and ε was a real number satisfying d(0,σ)(ε)=d(μ,σ)(ε) .

Analysis of response selectivity in artificial neurons

We explored response selectivity to bat calls in biomimetic neurons. To replicate the study performed on IC neurons (Salles et al., 2020), FMB and Echo calls in the sound database were used to measure responses on artificial neurons. Each audio clip was fed into the network after converting to auditory spectrogram (see above, Artificial network front-end processing), then we obtained activation probabilities for 100-nodes as σ(vc) in Figure 4. We averaged the activation probabilities over the same type of calls, and placed the results for 100-nodes onto a 2D scatter plot. Since the IC neurons are categorized into three groups, FMB selective, Echo selective, and non-selective (Salles et al., 2020), we performed k-means clustering (k = 3) on the principal axis by the principal component analysis (PCA).

Social call representation with artificial neurons

We explored bat’s call representation with the biomimetic network. In order to perform stochastic analysis, we made 10-copies for each audio clip in the natural sound database (see above, Social calls for natural sound representation) by a data augmentation based on temporal shift so that 260 audio clips were ready for the response analysis on artificial network. After converting the audio clips for 8 types of bat call to auditory spectrogram (see above, Artificial network front-end processing), the spectrograms were fed into the network to obtain the network’s responses, the vc in Figure 4. Then, we estimated the Gaussian distributions for the responses to each call type, and measured a distance between two distributions by using the Jensen–Shannon divergence (JSD) as JSD(P,Q)=(KLD(P,M)+KLD(Q,M))/2 , where P and Q represent two target distributions, KLD is the Kullback–Leibler divergence (KLD), and M=(P+Q)/2 (Endres and Schindelin, 2003). Unlike the KLD, the JSD is bounded within the interval [0,1] where 0 means that two distributions are equal. Finally, we quantified a discriminability across the classes by averaging JSDs of all cases choosing 2 of 8. In evaluation, we calculated the averaging JSDs with 10-models for each configuration that were trained on different initial values and summarized the mean and standard deviation of the 10 results. Additionally, we explored the noise effect on the sound representation with simulated audios produced by adding Gaussian random noise to each of the 260 audio clips depending on signal-to-noise ratio (SNR).

To compare with neural data, we performed this analysis between FMB and Echo responses. In the same manner, we calculated JSD based on the networks. We adopted neural data used in the previous study (Salles et al., 2020). Among 575 neurons, we chose 351 neurons which were recorded with same version of stimuli, and constructed 351D vector to represent response pattern across the neurons by concatenating the number of spikes on each neuron. Once the vector is projected onto 100D space based on PCA, we estimated Gaussian distributions for FMB and Echo responses. Then, we calculated JSD between the distributions.