Predictive Ensemble Decoding of Acoustical Features Explains Context-Dependent Receptive Fields

A primary goal of auditory neuroscience is to identify the sound features extracted and represented by auditory neurons. Linear encoding models, which describe neural responses as a function of the stimulus, have been primarily used for this purpose. Here, we provide theoretical arguments and experimental evidence in support of an alternative approach, based on decoding the stimulus from the neural response. We used a Bayesian normative approach to predict the responses of neurons detecting relevant auditory features, despite ambiguities and noise. We compared the model predictions to recordings from the primary auditory cortex of ferrets and found that: (1) the decoding filters of auditory neurons resemble the filters learned from the statistics of speech sounds; (2) the decoding model captures the dynamics of responses better than a linear encoding model of similar complexity; and (3) the decoding model accounts for the accuracy with which the stimulus is represented in neural activity, whereas linear encoding model performs very poorly. Most importantly, our model predicts that neuronal responses are fundamentally shaped by “explaining away,” a divisive competition between alternative interpretations of the auditory scene. SIGNIFICANCE STATEMENT Neural responses in the auditory cortex are dynamic, nonlinear, and hard to predict. Traditionally, encoding models have been used to describe neural responses as a function of the stimulus. However, in addition to external stimulation, neural activity is strongly modulated by the responses of other neurons in the network. We hypothesized that auditory neurons aim to collectively decode their stimulus. In particular, a stimulus feature that is decoded (or explained away) by one neuron is not explained by another. We demonstrated that this novel Bayesian decoding model is better at capturing the dynamic responses of cortical neurons in ferrets. Whereas the linear encoding model poorly reflects selectivity of neurons, the decoding model can account for the strong nonlinearities observed in neural data.


Introduction
The response properties of cortical neurons are typically characterized by their receptive fields (RFs) (Hubel and Wiesel, 1968;Theunissen et al., 2000). The underlying assumption is that each neuron acts as a feedforward filter of its stimulus, somewhat independently from other neurons. However, in addition to thalamic inputs, cortical neurons receive both excitatory and inhibitory intracortical inputs, strongly shaping their spontaneous and sensory-evoked activities (Lamme et al., 1998;Braitenberg and Schü z, 1998;Lee and Winer, 2008).
Several studies have shown that auditory spectrotemporal RFs (STRFs) can dynamically change depending on the stimulus characteristics (Theunissen et al., 2000;Valentine and Eggermont, 2004;Rabinowitz et al., 2011;Schneider and Woolley, 2011). This has led to more complex encoding models that incorporates contextual effects and nonlinearities into the RF analysis (Ahrens et al., 2008;Rabinowitz et al., 2012;David and Shamma, 2013;Willmore et al., 2016). Here, we hypothesize that the nonlinearities and context-dependency of neuronal responses can be described by a Bayesian network model with static decoding filters where neurons collectively and interactively aim to reconstruct their stimulus. Interactions between neurons allow for the recognition of sensory features that represents the input. Central to this inference mechanism is a basic reasoning argument (Wellman and Henrion, 1993) stating that a highly probable justification can "explain away" the stimulus, suppressing other possible interpretations. In particular, a sensory feature that is decoded (or explained away) by one neuron is not explained by another. Explaining away implies that selectivity of one neuron cannot be understood in isolation, outside of the context of other neurons representing overlapping stimuli. If one naively looks at the linear relation between the stimulus and the neural response (i.e., RF), this representation systematically changes and dynamically reshaped by surround stimuli. Instead, recovering the underlying fixed decoding filters can give a more stable description for the selectivity of each neuron. It can also provide conceptual motivation for why context-dependent, adaptive, gain-control models, such as (Ahrens et al., 2008;Rabinowitz et al., 2011Rabinowitz et al., , 2012, improve prediction performance over linear encoding models. The main goal of this work is to introduce new methods to investigate the characteristics of this stable selectivity in auditory neuronal data. We start with a simple normative approach to model the structure and dynamics of the auditory system. In this architecture, auditory neurons are considered as predictors rather than filters of their input, as they collectively aim at inferring the underlying causes of the sensory stimuli (Smith and Lewicki, 2006;Lochmann et al., 2012). We approximate the decoding filters of auditory neurons from the primary auditory cortex of ferrets and use these filters in our Bayesian model to explain and predict the dynamics of neural responses better than an equivalent encoding model. The stimulus can be reconstructed from model responses as accurately as it can be reconstructed from the neural responses, and much more accurately than it can be reconstructed from the responses of the encoding model. This suggests that auditory neurons collectively represent the auditory scene (thus optimizing the capacity of subsequent layers to linearly reconstruct the stimulus from the neural responses), rather than responding to the auditory signal with fixed, independent RFs (which would require highly nonlinear decoding). Finally, we show that the predictive fields (PFs) learned from the statistics of speech are highly similar to neurons' actual decoding filters and that explaining away plays a crucial role in recovering a diverse set of features that fully represent the sensory environment.

A generative model of sounds
To be able to perceive a complex sensory scene, we need a schema that explains the observed stimulus by describing its underlying causes (Friston, 2005); for instance, hearing a knock on the door indicates a visitor hitting the surface of the door with the knuckles of his fist, which is not directly observed. This schema, so-called generative model, is a crucial ingredient for inference as it describes the internal hypothesis about the origin of the stimulus (Dayan et al., 1995;Friston, 2010). The presumed goal of the sensory system is to optimize this internal generative model until it can accurately explain the sensory data.
The generative model we consider (Deneve, 2008) (see Fig. 1, bottom) assumes that auditory events correspond to the coactivation of particular patterns of frequency. These events start and stop randomly, and combine linearly to generate the mean auditory signal. This auditory signal is corrupted by signal-dependent noise, as commonly observed in sensory neurons.
In more technical terms, elementary features correspond to variables X ϭ ͕ x j ͖ that randomly switch on (x j (t Ϫ 1) ϭ 0, x j (t) ϭ 1) or off (x j (t Ϫ 1) ϭ 1, x j (t) ϭ 0) at each time step t, with constant probabilities r on and r off . Each feature contributes to the resulting stimulus S by activating the sensory receptor responses s i . More precisely, the base firing rate q i0 of a receptor s i increases by q ij Ն 0 as a result of the presence of feature x j . We call the effect of a feature x j to the set of all receptors the PF of that feature, that is, q ij Ն 0, i ϭ 1 … M where M is the number of receptors. Sensory receptors produce noisy (Poisson) spike trains forming the likelihood as follows: Without loss of generality, we remove the "baseline" firing rates q i0 from the equations below since it is equivalent to having an additional hidden feature that is always "on"; that is, let: where x 0 ͑t͒ ϭ 1. This likelihood together with the prior probabilities r on and r off forms the generative model. Throughout the manuscript, the stimulus will be in the form of frequency-time representation of sounds (spectrogram) where each frequency channel covers a small portion of the whole range from ϳ100 Hz to ϳ8 kHz. s i (t) represents a sensory receptor receiving input from the i-th frequency channel.

Bayesian inference and neural implementation of explaining away
The inference model computes the posterior probability of hidden features, P͑x ͉ s͒, by using the Bayes rule: P͑x ͉ s͒ ϭ P͑s ͉ x͒ P͑x͒/P͑s͒. For conceptualizing the inference model, we introduce the "feature detector units," which keep track of the probability of x j being present at time t denoted by p j ͑t͒ ϭ P͑x j ͑t͒ ϭ 1 ͉ s͑0 ¡ t͒͒, where s͑0 ¡ t͒ represents the input received from all receptors from time 0 up to t. Feature detectors can be thought as the auditory cortical neurons, and sensory inputs are either the stimulus (in its spectrographic form) or neural responses from a lower-level computational step.
We focus on the log-odds ratio of the hidden features as follows: Using the Bayes rule, the inference can be performed in an iterative process as follows: The terms in A j represent the ratio of the likelihood of observations, and the terms in B j represent the ratio of priors. Because of the Markov property, the state of the object j at t ϩ 1 depends only on its current state. The probability of x j being in an on-state at time t ϩ 1 is possible either by being on at time t and staying on or switching from off to on (similarly for the probability of being off at time t ϩ 1). Therefore, we can write the following: Since each receptor unit, s i (t), is assumed to be independent from other receptors and its own history, the ratio of likelihoods can be written as follows: where ij ϩ ϭ Α kϭ0,k j N q ik x k ͑t ϩ 1͒ ϩ q ij , that is, it is the i given in Equation  1 with the condition that x j (t ϩ 1) ϭ 1 and ij The exact inference would require the computation of 2 N feature configurations, and it is computationally expensive even for small models. Therefore, during the computation of likelihood of x j ͑t ϩ 1͒, we use a mean-field approximation and replace all other hidden features with their expected values p k ͑t͒ from the previous time step: ij ϩ ϭ Α kϭ0,k j N q ik p k ͑t͒ ϩ q ij , and ij Ϫ ϭ Α kϭ0,k j N q ik p k ͑t͒. Taking the logarithm of Equation 4 and simplifying, we obtain the following: Combining Equations 3 and 5, we conclude the iterative inference scheme for L j ͑t ϩ 1͒. Once L j ͑t ϩ 1͒ is known, we can use Equation 2 to compute p j ͑t ϩ 1͒ ϭ e Lj͑tϩ1͒ 1 ϩ e Lj͑tϩ1͒ . We assume that the detector neuron's activity is proportional to the rectified log likelihood (log A j ) rather than the posterior probability L j (t). The rationale behind this hypothesis is that sensory processing is inherently hierarchical. Thus, the output of the feature detector units will in turn be used as "sensory inputs" by similar units in later sensory stages. These higher-order detectors will themselves be integrators cumulating evidence for more complex (or more stable) features. To ensure selfconsistency, sensory neurons should thus transmit new sensory evidence in favor of their preferred features A j , rather than sensory evidence cumulated over time L j . Therefore, we take (using log͑x ϩ 1͒ Ϸ x for simplicity) the following: Each detector unit can thus be interpreted as performing a weighted sum of sensory detector responses (see Fig. 1, black plain connections). However, the strength of connection from a sensory receptor to a feature detector is not fixed but divided by a prediction of this input by other detector units, that is, ŝ i j ͑t͒ ϭ Α k j q i0 ϩ q ik p k ͑t͒ (see Fig. 1, magenta lateral connections). We can also interpret the feature detector as implementing divisive predictive coding as follows: where ͑s i /ŝ i j ͑t͒ Ϫ 1͒ is the prediction error from the i-th sensory receptor. In particular, when other detectors perfectly predict the receptor response, that is, s i ͑t͒ ϭ ŝ i j ͑t͒, the contribution of this synapse is completely cancelled. Explaining away in a neural layer thus takes the form of a divisive inhibition targeting the neuron's inputs.

Model with no inhibition
Throughout the manuscript, we compare the responses of the intact model with the responses of the model without inhibition. In this version, feature detectors do not inhibit each other and the stimulus is processed in a feedforward fashion. The likelihood portion of the inference (Eq. 5) takes the form log A j ϭ Α iϭ1 M s i ͑t͒log ͩ 1 ϩ q ij q i0 ͪ Ϫ q ij , and the rest of the inference remains the same.

Learning the PFs
Thanks to the normative approach, the PFs (i.e., the elementary, positive, and independent features composing the auditory scene) can be learned from the statistics of natural sounds (Lewicki, 2002;Smith and Lewicki, 2006). Once the online estimate of the p j ͑t͒,j ϭ 1 … N has been inferred on a temporal window t ʦ ͓0,T͔ (Expectation step), the PFs q ij can be learned by maximizing the log likelihoods of observations given the expected values of the features as follows: The derivative at the maximum of the log likelihoods is zero as follows: We solve this equation using stochastic gradient ascent in the simulations. We take T ϭ 1 for online learning (at the same time-scale as inference). Since we defined q ij to be positive, if q ij becomes negative at a gradient ascent step, it is assigned to zero.

Neuronal recordings and stimulus
A detailed explanation of the experimental procedures was given by Mesgarani et al. (2009). Four adult, female, awake, and head-restrained ferrets were used in the neurophysiological recordings. The protocol for all surgical and experimental procedures was approved by the Institutional Animal Care and Use Committee at the University of Maryland and was consistent with National Institutes of Health Guidelines. Speech stimuli were taken from the Texas Instruments/Massachusetts Institute of Technology (TIMIT) database (Garofolo et al., 1993). Thirty different sentences (3 s, 16 kHz sampling) spoken by different speakers (15 men and 15 women) were used to sample a variety of speakers and contexts. Each sentence was presented five times during recordings. Speech spectrograms were binned at 10 ms in time and in 30 logarithmically spaced spectral bins between 125 and 8000 Hz. The total dataset consisted of 128 single-neuron recordings. We used N ϭ 60 neurons throughout the paper satisfying a certain reconstruction accuracy threshold as described in more detail below.

Encoding and decoding filters
Encoding. Encoding is the mapping from stimulus (spectrogram) to neuronal responses. We compute the encoding filters (STRFs; see Fig. 3A left) using normalized reverse correlation (Theunissen et al., 2001). Let r͑t, n͒ ϭ Α i Α h i ͑, n͒s i ͑t Ϫ ͒ where r͑t, n͒ is the estimated response of neuron n, h i ͑, n͒ is the i-th frequency channel of its RF, is the delay, and s i (t) is the i-th frequency channel of the stimulus. h͑, n͒ is approximated by minimizing e n ϭ Α t ͑r͑t, n͒ Ϫ r͑t, n͒͒ 2 where r͑t, n͒ is the average firing rate of neuron n over trials. It can be shown (Theunissen et al., 2001) that h͑, n͒ ϭ C SS Ϫ1 C Srn where C SS is autocorrelation of the stimulus and C Srn is the cross-correlation between stimulus and neuronal response (for more explicit definitions, see Mesgarani et al., 2009). In practice, some form of regularization is required to prevent overfitting to noise. We used ridge regression where the regularization parameters are optimized using fivefold cross-validation in each training set.
Spectral RF (SRF) filters. SRF (see Fig. 3A, right) is a reduced form of the full encoding filters (no time dimension). To compute the SRF for each neuron, we computed the normalized reverse correlation while we fix the delay. We used fivefold cross-validation in each training set to optimize the regularization parameters. To find the right delay between recordings and SRF predictions (see Fig. 4A), we shifted the recordings by 10 ms from 0 to 200 ms and reported the maximum correlation coefficient. This is done for all competing models.
Decoding. Decoding is the mapping from neuronal responses to stimulus (spectrogram) (Bialek et al., 1991;Stanley et al., 1999;Mesgarani et al., 2009). These filters (see Fig. 3B, left) can be computed in a similar manner: Let ŝ i ͑t͒ ϭ Α n Α g i ͑, n͒r͑t Ϫ , n͒ where g i ͑, n͒ is the contribution of neuron n with delay to the reconstruction of the i-th frequency channel of the stimulus, ŝ i ͑t͒. For each i, g i ͑, n͒ is approximated by minimizing e i ϭ Α t ͑s i ͑t͒ Ϫ ŝ i ͑t͒͒ 2 . Then g i ͑, n͒ can be collected into a single structure g͑, n͒ that forms the decoding filter for neuron n. It can be shown that g i ͑, n͒ ϭ C RR Ϫ1 C RSi where C RR is the autocorrelation of neuronal responses and C RS s the cross-correlation between responses and the stimulus (Mesgarani et al., 2009). Similarly, we used ridge regression where the regularization parameters are optimized using fivefold cross-validation in each training set.
Extracted PFs from decoding filters. Similar to the procedure with encoding filters, we can also extract simplified decoding filters that resemble the theoretically defined PFs. For each neuron n, we restricted the above analysis to a single time delay to recover frequency-only decoding filters. We used fivefold cross-validation in each training set to optimize the regularization parameters. To find the right delay between real stimulus and its model reconstruction (see Fig. 5A), we shifted the stimulus by 10 ms from 0 to 200 ms and reported the maximum correlation coefficient. This is done for all competing models.
Selection of neurons for analysis. The dataset consists of 128 individual recordings from four animals. We combined the data from all animals into a single model assuming that similar stimulus selectivities would exist across animals. We were interested in a subset of neurons that carry relevant information about the stimulus. For this, we reconstructed the stimulus using each neuron's firing rates individually. If the reconstruction and actual stimulus had cross-correlation Ͼ0.2, the neuron was included for further analysis. This resulted in 62 neurons. We excluded two more neurons since they had completely negative decoding filters, which resulted in flat PFs.

Predicting neural activity
To use the model to get predictions for the neural activity (see Fig. 4), we obtained frequency-only decoding filters from the neural data as described in the previous section. We used these decoding filters as the PFs in the model. There were three free parameters that were adjusted globally to optimize the fit to neural data (r on ϭ 1 and r off ϭ 20 for each neuron and a global gain parameter for the output of the spectrogram). In all models (i.e., intact model) without inhibition and SRF, we found an optimum delay for each model prediction that maximizes the correlation coefficient to the neural responses on the training set. These delays were not changed during predictions over the test set. The dataset was divided into 10 equal portions over time, and all model prediction accuracies were averaged over 10-fold cross-validation process.

Transient and sustained responses
To identify responses after a sudden change of input to a neuron (see Fig.  6), we first identified the best frequency of each neuron from the maximum of its one-dimensional, frequency-only encoding filter. Then for each neuron, we searched for sudden increases in the spectrogram of the stimulus restricted to Ϯ10 frequency bands around the best frequency and found the timings of the largest 150 local maxima in the derivative of the stimulus averaged over frequencies. For local maxima, we used findpeaks.m function in MATLAB (The MathWorks) with the condition that two consecutive local maxima are at least 200 ms apart. This gave us 150 events for each neuron. Transient responses are defined as the first 100 ms after such an event. Sustained responses are the following 100 ms after a transient response.

Results
In the next few sections, we underline the main features of our normative model and provide evidence that auditory cortical neurons carry similar computations. We start with a simulation which shows that the RF (encoding filter) of a specific neuron can drastically change depending on the number of neurons in a network as a result of explaining away. Therefore, we expect that it would be challenging to predict the activity of a neuron without knowing the activity of other neurons in the network. Despite this difficulty, we wanted to test whether decoding filters could give a more reliable representation of neural activity. We approximated the decoding filters of several neurons in the primary auditory cortex of ferrets using linear regression and plugged these filters into our model as PFs. This allowed us to simulate neural activity that can accurately replicate the dynamics of actual neurons. We found that this simulated activity can reconstruct the stimulus much better than actual neural activity and the activity obtained from an encoding model. Finally, we show that PFs can be directly learned from the statistics of stimulus and that explaining away is crucial in obtaining variety of filters.

Network effect on a single neuron
Despite the simplicity of the generative model (Fig. 1), the predicted unit responses (i.e., firing rates as given in Eq. 6) are highly nontrivial, and, in particular, cannot be understood in isolation, independently of other units. If, for example, we remove some of the PFs from the pool of available detectors (analogous to inactivating the corresponding sensory neurons), we find that the STRFs change (Fig. 2). This predicts that identifying the relevant features represented by sensory neurons and predicting their responses will be extremely challenging. In general, reaching perfect accuracy would require the simultaneous recordings of all neurons involved. Only a small subsets of all auditory neurons are recorded and included in this study, and they were not recorded simultaneously. Our primary goal in the next sections will thus be to clarify the functional and mechanistic implications of explaining away and accounting for particular aspects of the responses rather than focusing on precisely reproducing spike trains (as, e.g., in a GLM approach) (Paninski, 2004).
This normative approach illustrates three important predictions for auditory responses. First, if neurons are indeed involved in the optimal detection of overlapping auditory sounds, their responses will most likely be very poorly predicted by a linear encoding model. This is because the predicted responses are rendered highly interdependent by explaining away. Second, despite their apparent variability and absence of a fixed encoding model, the sound stimulus should be well reconstructed by a linear decoder applied to auditory responses. This linear decoder would indeed provide a better picture of neural selectivity than a linear encoding model. And third, even with the "right" decoding model, the normative model is unlikely to predict the neural responses with very high accuracy. This is because only a small subset of neurons are recorded and the rest are acting as "hidden variables" that affect the recorded neurons in unpredictable ways.
Encoding versus Bayesian decoding models applied to neural data In the auditory system, the response properties of neurons are usually characterized by STRF, which is the linear time-frequency filter best fitting the neuronal response when applied to stimulus (Aertsen and Johannesma, 1981). STRF describes how  stimulus is encoded in the neuronal activity; therefore, it is also called the encoding filter (for examples, see Fig. 3A, left). Another way to look at the relationship between neuronal data and stimulus is to ask how stimulus can be decoded from the neuronal activity (Bialek et al., 1991). This method finds the linear mapping from the response of a population of neurons to the stimulus (see Materials and Methods), and the resulting frequency-time filter for each neuron is called the decoding filter (Fig. 3B, left). Our Bayesian normative model suggests that context-dependent encoding filters could be replaced by the static decoding filters approximated from neuronal data, and then be used to predict the neuronal responses. To test this idea, we combined neuronal recordings from the primary auditory cortices of 4 ferrets listening to a 90 s speech stimulus from the TIMIT database (Garofolo et al., 1993) (see Materials and Methods). For each neuron, we computed a one-dimensional (frequencyonly) approximation of the decoding filter that summarizes the full decoding filter (Fig. 3B, right). We call this filter the PF of the corresponding auditory neuron as it represents the contribution of the neuron in predicting the stimulus. Then, we plugged these PFs into the model as q ij s (no learning is involved). r on and r off rates were assumed equal among the neurons and chosen to optimize the fit between model and real neural responses. We interpreted the firing rate of each model neuron (Eq. 6) as a prediction for the firing rate of the corresponding real auditory neuron. As a baseline model, we also generated predictions from a model where inhibitory connections (i.e., explaining away) are removed and call it the model Figure 4. Prediction of neuronal activity by different models. A, Comparison of average correlation coefficients (10-fold cross-validated) between the activity of 60 auditory neurons and their predicted responses by the intact Bayesian decoding model, SRF model, and Bayesian model without inhibition. Bayesian decoding model and SRF predictions are comparable (left; cc ϭ 0.31 for both) and provide a better fit to real responses compared with the model without inhibition (cc ϭ 0.26, p Ͻ 10 Ϫ10 , paired t test). For comparison, full STRFs reach a performance of cc ϭ 0.38, significantly higher than the simplified encoding and decoding models. However, the first and second derivatives of the intact Bayesian decoding model responses (middle and right panels) fit the neural responses' first and second derivatives significantly better than the SRF predictions ( p Ͻ 0.05 and p Ͻ 10 Ϫ6 , respectively). The Bayesian model better captures the detailed dynamics of the responses. Error bars indicate SEM. B, Comparison of measured neural responses for two example neurons (gray in both panels) with their inference model predictions with or without inhibition (blue and red, respectively), SRF model prediction (green), and STRF prediction (yellow). C, Cross-covariance between inference model versus real responses and between SRF versus real responses shows that model responses are temporally more precise than SRF responses. Shaded areas represent the 95% confidence intervals. *p Ͻ 0.05; ****p Ͻ 10 Ϫ4 . without inhibition. This altered model simulates a feedforward-only version of the intact model.
For comparison, we also computed a simplified encoding model for the neural data with the same basic structure. We reduced the full encoding filter of each neuron to one-dimensional (frequency-only) filter (Fig. 3A, right), which we call the SRF, convolved it with the stimulus, and rectified it to predict the neuronal activity.

Performances of the encoding and Bayesian decoding models
The average cross-validated performances of the encoding and Bayesian decoding models in accounting for the neural firing rates are shown in Figure 4A. The model predictions and SRF predictions are comparable and better than the predictions of model without inhibition (Fig. 4A, left; p Ͻ 10 Ϫ10 , paired t test, N ϭ 60 neurons). However, the ability of the Bayesian decoding model to capture the details of dynamic responses becomes more obvious when higher derivatives (first and second) of the neural responses and their predictions are compared (Fig. 4A, middle and right). The dynamics of the Bayesian model fits significantly better to the first ( p Ͻ 0.05) and second ( p Ͻ 10 Ϫ6 ) derivatives of the responses compared with the encoding (SRF) model. This is the case even though the Bayesian decoding model is indirect, as it optimizes the reconstruction of the stimulus, not the neural responses as the encoding model does. For comparison, the full STRFs that take the temporal dimension into account (and have 750 free parameters each, even though these parameters are effectively decreased to ϳ155 using ridge regression and only counting the parameters whose magnitude are Ͼ1 SD of the mean in each STRF) reach a performance of correlation coefficient ϭ 0.38, significantly higher than the simplified encoding and decoding models with 31 (frequency channels and optimum delay parameters) and 34 (additional r on , r off and a global gain parameter for the PFs) parameters, respectively. In Figure 4B, we show the fits of the model (blue), model without inhibition (red), SRF (green), and STRF (yellow) predictions for the activity of two auditory neurons (gray). The model predictions are better at capturing the transient dynamics of neuronal responses than SRF and model without inhibition. Temporal precision of the model is higher in the entire dataset, as illustrated by the mean cross-covariance of the model versus real responses (Fig. 4C, blue) and SRF versus real responses (Fig.  4C, red).
We also compared how well the stimulus can be reconstructed from the model neural responses using the full decoding filters (as in Fig. 3B, left) obtained from the real neural data. We took the responses of the Bayesian decoding model, SRF model, and model without inhibition and convolved them with the full decoding filters of the real neurons to reconstruct the stimulus. As summarized in Figure 5, we can observe that the Bayesian decoding model massively outperforms the encoding model. The stimulus can be reconstructed very accurately from the responses predicted by the Bayesian model (cc ϭ 0.81). It is as good as the stimulus reconstructed directly from the neural data (cc ϭ 0.79) or from the STRF predictions (cc ϭ 0.80). If, on the other hand, we predict neural responses from their SRFs, and then try to reconstruct the stimulus, the performance is poorer (cc ϭ 0.61) and clearly inferior to a reconstruction directly based on the neural data ( p Ͻ 10 Ϫ10 , paired t test, N ϭ 30 frequency channels).
In summary, the Bayesian decoding model outperforms the encoding model in how accurately it can represent the stimu-lus. Moreover, it is as good, if not better, than an encoding model of similar complexity when it comes to predicting neural responses, even when it is not directly fitted to those responses.

Response dynamics are captured by explaining away
We also observed other signatures of online inference, both in the neural data and in the Bayesian decoding model. We compared the transient and sustained responses of each neuron after the occurrence of a sudden transient in the stimulus (a sharp increase in the mean power averaged over frequency), such as the beginning of a new word or sentence (see Materials Comparison of reconstruction accuracies averaged over frequency channels (10-fold crossvalidated) for the neural responses, SRFs, STRFs, and inference models. Intact model reconstruction (cc ϭ 0.81) and STRF reconstructions (cc ϭ 0.80) are as good as the reconstruction by the real responses (cc ϭ 0.79), which is superior ( p Ͻ 10 Ϫ10 ) to the SRF and inference model without inhibition reconstruction (cc ϭ 0.61). Error bars indicate SEM. N ϭ 30 frequency channels, paired t test. B, Short samples of reconstructed stimuli from model, SRFs, and neural responses as well as the original stimulus. ****p Ͻ 10 Ϫ4 and Methods). Transient response is defined as the initial 100 ms after such an event, and sustained response is the next 100 ms after a transient response (for the average time course of responses after events, see Fig. 6C). Such dynamics are also consistent with the full STRF model. Immediately following these events, we observe a gradual sharpening of the RFs. Comparison of average RFs during transient and sustained responses for real neurons (Fig. 6A) shows the increased frequency selectivity as time progresses. A similar phenomenon is also replicated by the Bayesian model (Fig. 6B) due to the divisive competition between feature detectors. Following the occurrence of an event, many neurons are activated as they pool inputs from highly overlapping PFs. Transient response is strong but weakly selective. However, as more information is accumulated about the sound, some neurons become more confident than others about the presence of their preferred sound. This rise in divisive inhibition leads to decay in firing rates after transient response and sharpens the selectivity.

PFs, not RFs, reflect neural selectivity
Bayesian model neurons receive feedforward inputs pooled from their PFs (corresponding to the real decoding filters, Fig.  3B, right). However, their encoding filters (Fig. 3C) are reshaped by the competing activity of other model neurons. They are therefore narrower than the original PFs, with negative parts that are inexistent in the purely positive PFs (Fig. 3B, right vs Fig. 3C). The decoding filters of model neurons are wider and smoother than their encoding filters (Fig. 3C vs Fig. 3D). This phenomenon is also observed in the actual neural data. In particular, the measured encoding filters are markedly narrower than the corresponding decoding filters (Fig. 3A vs Fig. 3B), with additional negative subfields. To quantify this observation and present population trends, we plotted (Fig. 7A) the average encoding and decoding filters centered around the best frequency of each neuron (as long as the best frequencies of encoding and decoding filters are within Ϯ3 frequency channels of each other, which was the case for 41 of 60 neurons). Model neurons reproduce this qualitative difference on average (Fig. 7C), that is, wide, positive decoding filters and narrow, positive/negative encoding filters.
If inhibitory effects between neurons are removed (model without inhibition), the encoding filters became wider (Fig. 7B, middle) and resembled the underlying feedforward connections (Fig. 7A, middle). Therefore, these effects are entirely due to explaining away. Observed differences between encoding and decoding filters of auditory neurons support the hypothesis of broad feedforward tuning (through PFs) that is narrowed by lateral inhibition (  t test). These results show that the competition between feature detectors is crucial for shaping their frequency selectivity and efficiently representing natural stimuli.

Discussion
We presented a normative model where neurons interactively and collectively predict their input to optimally represent the underlying causes of their observations. A signature of this efficient representation is the competition between neurons representing similar stimulus features. The model predicts that feedforward inputs to a neuron (through PFs) are partially explained away by competing neurons which results in more narrow encoding filters (STRFs) compared with decoding filters. We also observed this signature in the analyzed auditory neurons (Figs. 3, 7). The dynamics of model responses are similar to neural data (Fig. 6) and can capture the activity better than a comparable feedforward model (Figs. 4,5). In addition, we show that the underlying selectivity of some auditory neurons (PFs) can be accurately captured from the statistics of natural stimuli using the model's inference and learning algorithms (Figs. 8,9). Despite originating from a functional hypothesis, the model is able to shed light on some fundamental mechanisms of sensory processing.
Biological plausibility of the model Identifying specific neural mechanism for inference and learning is not the main focus of this article. The model is normative, constrained only by optimality principles, and thus largely independent of its specific neural implementation. As an example, the algorithm presented here could be implemented in the auditory cortex by "projection" neurons (excitatory) representing the fractional prediction errors (s i /ŝ i j Ϫ 1) and inhibitory interneurons pooling the activity of projection neurons and inhibiting them divisively. In turn, our model neuron responses (the feature detector) could be obtained by summing the responses of these "projection" neurons. In any case, we note that the prediction errors (s i /ŝ i j Ϫ 1) used to update the activity of model units is zero on average, and is automatically kept as close to zero as possible by the inference algorithm. Thus, the experimentally observed balance between excitatory and inhibitory currents into cortical neurons (Wehr and Zador, 2003;Zhang et al., 2003;Okun and Lampl, 2008) could be interpreted as an imprint of an efficient network representing its input (Bourdoukan et al., 2012;Boerlin et al., 2013).
Moreover, explaining away provides a functional rationale for divisive normalization in which a neuron's response to a stimulus is given by its driving input divided by the summed activity of a pool of nearby neurons (Heeger, 1992;Heeger et al., 1996;Rabinowitz et al., 2011). This nonlinear phenomenon was proposed as a canonical computation for the brain because of the role it plays in the retina, visual cortex, olfaction, auditory system, and multisensory integration (Carandini and Heeger, 2011). Here, we argued that divisive normalization can naturally emerge in a neuronal system as a result of internal predictions through inputtargeted inhibition. Normalization occurs between neurons that share similar feature selectivity, and it acts as a mechanism to reduce redundancy (Schwartz and Simoncelli, 2001;Lyu and Simoncelli, 2009).
In the present model, the PF of feature detectors can be learned by maximizing the log-likelihood, which requires only the input from presynaptic receptor cells and postsynaptic activity of feature detectors. Furthermore, these PFs can be learned in an online fashion. Learning could be implemented locally by Hebbian-like plasticity (Caporale and Dan, 2008) (black) can be learned by model neurons (red) with an Expectation-Maximization algorithm starting from an initial condition (light blue) with a small bump at the peak of auditory neurons' PFs. B, Inhibition in the model is crucial for correctly approximating real PFs. The cc between real versus learned PFs in the intact model (horizontal axis, mean cc ϭ 0.68) is significantly higher than the cc between real versus learned PFs in the model without inhibition (vertical axis, mean cc ϭ 0.17). Inset, In addition to the inference model and inference model without inhibition (no inh) performances, we provide the mean cc between initial PFs (A, light blue) versus real PFs (initial, mean Ϸ 0.43) and the mean cc between two arbitrary real PFs (random, mean Ϸ 0). **p Ͻ 10 Ϫ2 ; ****p Ͻ 10 Ϫ4 .

Encoding, decoding, and correlations
What accounts for the experimentally observed differences between encoding and decoding filters? As detailed in Materials and Methods, both encoding and decoding filters are based on the cross-correlation between the stimulus and neuronal activity. For encoding filters, this cross-correlation is normalized by the autocorrelation of stimulus, whereas for decoding filters it is normalized by the autocorrelation of neuronal responses (Mesgarani et al., 2009). Therefore, narrow encoding filters can be explained by the high autocorrelation of the natural stimuli (specifically speech) compared with the relatively low autocorrelation of neuronal responses. This relatively low autocorrelation cannot be explained entirely by neural noise because we used the average firing rate over five trials. We proposed explaining away as a possible reason for decorrelated neuronal responses. We showed that, if this inhibitory effect is removed from model neurons, their encoding filters become wider and resemble their decoding filters more closely (Fig. 7). Such extension of frequency tuning has been shown previously by blocking cortical inhibition pharmacologically (Chen and Jen, 2000;Wang et al., 2002) and supports the mechanism proposed in the model. This prediction could be further tested using existing optogenetic techniques: We would expect that silencing the inhibitory interneurons (such as parvalbumin-positive neurons) (Moore and Wehr, 2013) in the auditory cortex would result in broadening of the frequency tuning of the remaining excitatory pyramidal neurons. We would further predict that encoding and decoding filters of these excitatory cells would become more similar and their spike trains would get more correlated after suppression of inhibitory neurons.

Relation to previous models
Predictive coding in the framework of Bayesian inference has been used successfully to explain classical and nonclassical RF properties in the visual domain such as end-stopping behavior (Rao and Ballard, 1999), dynamic adaptation to image statistics (Hosoya et al., 2005), and basic tuning properties and contextual effects (Spratling, 2010;Lochmann et al., 2012). However, it has been rarely used to explain neuronal phenomena in the auditory domain (Turner and Sahani, 2007;Ramirez et al., 2011;Yildiz and Kiebel, 2011;Wacongne et al., 2012;Yildiz et al., 2013) and, as far as we know, has never been used for predicting activity at a single-neuron level as achieved here.
Our functional approach assumes that single neurons represent specific events in the external world, and we deduce the neuronal mechanism from basic probabilistic principles. In contrast to previous functional approaches, such as efficient coding hypothesis (Olshausen and Field, 1996;Smith and Lewicki, 2006), we can learn and recognize stimuli in an online fashion, mimicking the biological auditory processing more closely, and we propose a plausible neural architecture with local, online learning rules for the parameters.
Our model is limited by its simple Markovian generative model. Markovian dynamics were chosen to introduce basic integrative properties in the model neurons while keeping the number of free parameters at a minimum. These dynamics could be improved with a more detailed generative model specific to the auditory system (Yildiz et al., 2013) and with inclusion of hierarchies representing different time scales and more complex features (Kiebel et al., 2008;Yildiz and Kiebel, 2011).
Top-down, normative modeling approach followed here is complementary to bottom-up, descriptive approaches, such as Linear-Nonlinear cascade (Chichilnisky, 2001), GLMs (Paninski, 2004), contextual RFs (Ahrens et al., 2008), gain control and adaptation (Rabinowitz et al., 2011;Rabinowitz et al., 2012;David and Shamma, 2013;Willmore et al., 2016), stimulus surprise (Gill et al., 2008), and maximally informative dimensions (Sharpee et al., 2004;Atencio et al., 2012). In particular, once these models are tuned for a specific type of stimulus, such as artificial or noise stimuli, they usually have to be retuned to explain neural responses with other types of stimuli. Instead, our model aims to identify the fixed selectivity of each neuron, independently of stimulus context and therefore presumably generalizable to different types of stimuli. Using protocols specifically designed to estimate these PFs rather than the stimulus-specific RFs of sensory neurons could unmask the underlying richness of these representations.