## Abstract

Optimal use of sensory information requires that the brain estimates the reliability of sensory cues, but the neural correlate of cue reliability relevant for behavior is not well defined. Here, we addressed this issue by examining how the reliability of spatial cue influences neuronal responses and behavior in the owl's auditory system. We show that the firing rate and spatial selectivity changed with cue reliability due to the mechanisms generating the tuning to the sound localization cue. We found that the correlated variability among neurons strongly depended on the shape of the tuning curves. Finally, we demonstrated that the change in the neurons' selectivity was necessary and sufficient for a network of stochastic neurons to predict behavior when sensory cues were corrupted with noise. This study demonstrates that the shape of tuning curves can stand alone as a coding dimension of environmental statistics.

**SIGNIFICANCE STATEMENT** In natural environments, sensory cues are often corrupted by noise and are therefore unreliable. To make the best decisions, the brain must estimate the degree to which a cue can be trusted. The behaviorally relevant neural correlates of cue reliability are debated. In this study, we used the barn owl's sound localization system to address this question. We demonstrated that the mechanisms that account for spatial selectivity also explained how neural responses changed with degraded signals. This allowed for the neurons' selectivity to capture cue reliability, influencing the population readout commanding the owl's sound-orienting behavior.

## Introduction

Uncertainty in sensory perception depends on the reliability of the sensory cues used to infer properties of the environment. The reliability of a sensory cue can be viewed as its precision across different environmental states (Fetsch et al., 2012; Ma and Jazayeri, 2014); that is, how consistent the cue is across trials and contexts. It is important for the brain to estimate the reliability of sensory cues so that decisions can be based on the most reliable signals and prior knowledge can be relied upon more heavily when sensory signals are degraded. Change in the reliability of sensory signals is inherent to natural environments. For example, reverberation (Devore et al., 2009) or extraneous noise sources (Keller and Takahashi, 2005) changes the signal-to-noise ratio of the auditory input (Saberi et al., 1998). In fact, high signal-to-noise ratio is unlikely in natural auditory scenes (Młynarski and Jost, 2014). The environmental noise, shared among neurons, can dramatically limit the information about the sensory cue encoded by neural populations (Moreno-Bote et al., 2014; Pitkow et al., 2015). Many studies have demonstrated that sensory reliability influences behavior (Ma and Jazayeri, 2014). For example, decreasing the reliability of spatial cues for signaling sound source direction increases the owl's sound localization bias (Saberi et al., 1998; Fischer and Peña, 2011). Although it is clear that sensory reliability influences behavior, how it is encoded by neurons is a matter of debate.

Several models of probabilistic computation have proposed how cue reliability could be represented in neuronal responses. These models make different assumptions about the role of tuning properties and spiking variability in coding sensory reliability. A population vector (PV) implementation of Bayesian inference with heterogeneous populations (Fischer and Peña, 2011; Girshick et al., 2011) hypothesized that sensory reliability is encoded in the neurons' spatial selectivity; that is, the sharpness of the tuning curves, and does not require any particular distribution of spiking variability. Conversely, in probabilistic population codes and sampling-based theories, sensory reliability affects the population activity through the entire shape of the tuning curves, along with their gain and correlated variability (Pouget et al., 2003; Jazayeri and Movshon, 2006; Ma et al., 2006; Fiser et al., 2010).

Here, we investigate the neural correlates of cue reliability in the space-specific neurons of the owl's midbrain. To localize sounds in the horizontal space, owls and other animals rely on the difference in the arrival time of a sound at the ears, called the interaural time difference (ITD) (Strutt, 1907; Blauert, 1997). ITD is initially processed in the owl's brain by neurons that compute the interaural phase difference (IPD) in narrow frequency bands (Carr and Konishi, 1990). Estimating sound direction from IPD carries a great deal of uncertainty due to corruption of the IPD by concurrent sounds (Keller and Takahashi, 2005). The owl's sound localization system is well suited to study how sensory reliability affects neural responses and behavior for several reasons. First, it is well established that the IPD is the relevant cue for the owl's localization in horizontal space (Moiseff and Konishi, 1983). Second, the spatial dependence of the cue variability in the presence of concurrent sounds has been estimated (Cazettes et al., 2014). Finally, noise in IPD can be quantitatively manipulated (Albeck and Konishi, 1995; Saberi et al., 1998).

To determine how cue reliability is represented in neural responses, we first examined the response of midbrain space-specific neurons as cue reliability varies. We found that changing cue reliability affects the response gain, the selectivity, and the correlated variability of neurons. We demonstrated that this coding emerges from the IPD detection mechanism and direction-dependent frequency tuning. Then, we showed that, when cue reliability varies, the broadening of the tuning curves was necessary and sufficient for a network of stochastic neurons to estimate sound direction near optimally and predict behavior.

## Materials and Methods

#### Electrophysiology

Surgery was performed as described previously (Wang et al., 2012). Briefly, adult barn owls (2 females) were anesthetized with intramuscular injections of ketamine hydrochloride (20 mg/kg; Ketaset) and xylazine (4 mg/kg; Anased). It has been shown that the responses of midbrain neurons are remarkably stable under ketamine anesthesia (Ter-Mikaelian et al., 2007; Schumacher et al., 2011). Prophylactic antibiotics (ampicillin; 20 mg/kg, IM) and lactated Ringer's solution (10 ml, SC) were injected at the beginning of each session. Carprofen (3 mg/kg, Rimadyl, IM) was administered to prevent inflammation and pain. These procedures complied with National Institutes of Health and guidelines of the Albert Einstein College of Medicine's Institute of Animal Studies.

The inferior nucleus of the external colliculus (ICx) was located stereotaxically and by the characteristic responses to ITD and interaural level difference (ILD) (Peña and Konishi, 2001). Responses were recorded using 1 MΩ tungsten electrodes (A-M Systems). Tucker Davis Technologies (TDT) System 3 and MATLAB software (The MathWorks) were used to record neural data.

#### Sound stimulation

All experiments were performed inside a double-walled sound-attenuating chamber. Custom-made earphones consisting of a calibrated speaker (Knowles model 1914) and a microphone (Knowles model EK-23024) were inserted in the owl's ear canal (Wang et al., 2012). Auditory stimuli delivered through the earphones consisted of 100 ms signals with a 5 ms rise–fall time at 10–20 dB above threshold. For each ICx neuron, we first measured the ITD, ILD, and rate-level response with broadband noise (0.5–10 kHz) and frequency tuning with tones. ITD was initially varied from ±300 μs with 30 μs steps over 5 trials and sound level varied from 0 to 70 dB SPL with 5 dB steps. Frequency tuning was estimated with tones ranging from 600 to 9000 Hz varied in 200 Hz steps over 15–20 trials. ITD tuning within the main peak of the rate–ITD curve was recorded at a finer resolution (10 μs steps; 20 repetitions) at the best ILD to measure the ITD-tuning width. Finally, neural responses to ITD with various tonal stimulations within the excitatory frequency range (±300 μs; 30 μs steps; 20 repetitions) were measured. Stimuli within all tested ranges were randomized during data collection.

To manipulate the binaural correlation (BC), three random noises (N1, N2, and N3) were generated on the computer. N1 was delivered to one ear and its copy with a time shift was delivered to the other ear. N2 was added to N1 in one ear and N3 was added to N1 in the other ear. These additional noise signals reduced the correlation between the signals in the two ears depending on the relative amplitudes of the uncorrelated and correlated noises. BC was calculated from BC = 1/(1 + *k*^{2}), where *k* is the ratio between the root-mean-square amplitudes of the uncorrelated and correlated noises (Licklider and Dzendolet, 1948; Jeffress and Robinson, 1962; Saberi et al., 1998). On a subset of neurons, we obtained ITD curves at the finer resolution (10 μs steps; 10 repetitions) using noise bursts at different BC levels (from 0 to 1 in 0.1 steps).

#### Data analysis

Wave_clus was used for spike sorting (Quiroga et al., 2004). Briefly, spikes were detected using a voltage threshold set at five times the estimated SD of the signal. To avoid double detection, spikes were separated by at least 1 ms. Neurons were considered isolated based on the presence of a refractory period of >1 ms in the interspike interval histogram and the similarity of spike shape. A complementary quality metric was the non-overlap of wavelet coefficients. In addition, the results of the sorting algorithm were visually inspected to confirm the quality of the sorting. Consistent with previous reports (Winkowski and Knudsen, 2006), no significant differences were found between the results of sorted and nonsorted traces.

For each stimulus parameter, a rate curve was computed by averaging the firing rate during the stimulation over trials. The spontaneous firing rate, obtained when no sound stimulation was presented, was subtracted from the response to the stimulus. The widths of the ITD-tuning curves were measured at half the distance between the minimum and the maximum response. Because the ITD tuning in owl's neurons are nearly symmetrical (Viete et al., 1997; Pérez and Peña, 2006), the center of the range used for the width measurement was used to estimate the preferred ITD. This method was used to avoid sampling bias in peak responses and was not significantly different from the best ITD measured at the maximum response (Wilcoxon rank-sum test, *p* < 10^{−10}). We used the absolute value of the best ITDs to combine data from contralateral and ipsilateral sides as a function of the eccentricity of the receptive field. We have previously tested different methods to quantify frequency selectivity, which yielded similar results (Cazettes et al., 2014). Therefore, here, we defined the best frequency as the center of the frequency range that elicited >30% of the maximum response (Cazettes et al., 2014).

In the owl, ITD increases monotonically with eccentricity (Moiseff, 1989), from small values in the front ITDs (∼0 μs) to the largest values in the periphery (∼200 μs). We estimated azimuth from ITD for the owl using published measurements (Moiseff, 1989) to compare the width and maximum firing rate of ICx neurons as a function of ITD with cue variability as a function of azimuth.

Variability of IPD at the input to the ears is dominated by the presence of concurrent sounds, which can shift the IPD of the target source. The filtering effect of the head causes IPD variability to be inherently direction and frequency dependent in the presence of concurrent sounds coming from different directions (Cazettes et al., 2014). Furthermore, the IPD variability may change depending on the amount of noise in the environment; that is, the signal-to-noise ratio. We quantified IPD variability across space for concurrent sounds as described previously (Cazettes et al., 2014). Briefly, to compute the IPD variability for concurrent sound sources, broadband noise signals with flat spectra between 0.5 and 9 kHz and equal level (50 dB) were each convolved with head-related impulse-responses of 10 barn owls at the appropriate source direction. The outputs were passed through a gamma-tone filter bank with center frequencies ranging from 1 to 8 kHz in 0.2 kHz steps and equal gains across frequency. The bandwidths of the gamma-tone filters were estimated from (Köppl, 1997) to match the bandwidths of the owl's cochlear filters, as described in (Fischer et al., 2009). The IPD in narrow frequency channels was calculated as the phase delay with the highest value of the cross-correlation between the left and right outputs of the gamma-tone filter bank. For each target direction, we obtained 37 estimates of IPD, each with a second source located at one of the directions covering the frontal hemisphere (between −90 and 90 degrees in steps of 5 degrees at elevation zero). The circular SD of IPD over directions of the second sound source was used as the estimate of IPD variability at each direction (Batschelet, 1981). To quantify the IPD variability at each direction, we averaged the SD of IPD across frequency.

#### Model of space map population

To model the population response in the owl's midbrain (Fig. 1), single broadband noise signals between 0.2 and 12 kHz were first passed through a gamma-tone filter bank with center frequencies ranging from 0.25 to 8 kHz in 0.25 kHz steps and equal gains across frequency (Köppl, 1997). The outputs were cross-correlated and transformed by an exponential input–output function to produce a set of frequency-dependent ITD curves (from 0.5 to 8 kHz). This step simulates the response of neurons directly upstream of ICx (Fischer et al., 2009). The input to ICx was then modeled as a weighted sum of the frequency-dependent ITD curves. The weights were direction dependent to model the direction-dependent frequency tuning in ICx (Cazettes et al., 2014). The weights for each neuron have a Gaussian shape with a center and width selected to match the observed frequency tuning at the preferred direction of the neuron. The weighted sum of the frequency-dependent ITD curves was passed through an exponential nonlinearity to produce the firing rates of Poisson neurons that modeled ICx responses (Fischer et al., 2009). The responses of space-specific neurons in the owl's optic tectum, which receive point-to-point projections from ICx (Knudsen and Knudsen, 1983), are influenced by feedback inhibition that serves to select the dominant stimulus when multiple peaks of activity are present in the map (Mysore and Knudsen, 2012, 2014). To model the output of the midbrain, we used a simplified version of the model of Mysore and Knudsen (2012, 2014) by multiplying the ICx activity pattern by a rectangular window function to select the dominant peak. The window is equal to one on an interval that is 28 degrees wide and centered on the location of the peak in a smoothed version of the ICx activity pattern and is zero outside this interval. The decoding analyses are applied to the windowed response pattern representing the population activity of the optic tectum.

#### Predictions for coding reliability in optimal and suboptimal decoding approaches

We considered the properties of optimal and suboptimal decoders in the context of a perceptual inference task consisting of inferring an unknown property of the environment, denoted θ, from a sensory cue *S* extracted from the sensory input. For the sound localization task, θ is the direction of a target sound source and *S* is the set of IPD cues across frequency channels computed in the sound localization pathway. The sensory cues determine the responses of a population of neurons, denoted *r*. In our study, *r* represents the responses of ICx neurons. For a given θ, the sensory cues *S* may vary due to environmental factors such as the presence of background noise (Cazettes et al., 2014). The variance of the sensory cues *S*, given θ, Var(*S*|θ), is used to quantify this variability. Similarly, the inverse of the variance is used to quantify the cue reliability (Landy et al., 2011). The cue reliability will have an effect on inferences made about θ by allowing the most reliable cues to dominate the estimate and by producing estimates that depend more on prior knowledge when sensory cues are unreliable. In the context of Bayesian inference, we expect to see an increase in behavioral bias when cue reliability decreases. This occurs because the likelihood function gets wider when cue reliability decreases and the prior has more influence on the posterior and thus on the estimate.

##### Population vector.

The suboptimal decoder considered was the PV, which has been used previously to describe the owl's sound localization behavior (Fischer and Peña, 2011). The PV is the average preferred direction of a population of *N* neurons, weighted by their respective spike counts as follows:
where the response of the *n*^{th} neuron is denoted *r _{n}* and

*u*(θ

*) is a unit vector pointing in the*

_{n}*n*

^{th}neuron's preferred direction θ

*. The model of Fischer and Peña (2011) proposes that the neural activities satisfy specific properties that allow the PV to accurately approximate a Bayesian estimate from the posterior*

_{n}*p*(θ|

*S*). It is assumed that the preferred directions of neurons in the population θ

_{n}are drawn from the prior distribution

*p*(θ). It is also assumed that the mean firing rates of the neurons are proportional to the likelihood function

*p*(

*S*|θ). That is, the mean firing rate of the neuron with preferred direction θ

_{n}can be written

*f*

_{n}(

*S*) ∝

*p*(

*S*|θ

_{n}). The response

*r*

_{n}(

*S*) on any trial may be modeled as a Poisson variable with mean

*f*

_{n}(

*S*). If these assumptions are satisfied, then the PV will point in the same direction as the Bayesian estimate from the posterior

*p*(θ|

*S*) (Shi and Griffiths, 2009; Fischer and Peña, 2011).

The hypotheses for how the population is constructed to allow for the PV to accurately approximate the Bayesian estimate determine how sensory cue reliability is encoded in the neural responses. Sensory cue reliability is captured by the likelihood function *p*(*S*|θ), which here is assumed to determine the response of the population to the stimulus. If the cue reliability decreases, then the likelihood function *p*(*S*|θ) will be wider and thus the population activity will be spread over more neurons. To produce this spread of population activity for each stimulus, the individual neurons must respond to more stimulus values, so the tuning curves will be wider. Conversely, the overall gain of the population response does not represent any property of the likelihood function *p*(*S*|θ) in this model. Therefore, the overall gain is not assumed to encode stimulus reliability.

The model also yields predictions for how tuning width and gain are incorporated into the neural decoding of θ. Increasing the widths of the neural tuning curves throughout the population results from the activity being distributed over more neurons. This allows the nonuniform distribution of preferred directions to increase the bias in the PV toward the region of highest concentration of preferred directions. Alternatively, changing the overall gain of the tuning curves will not influence the average behavior of the PV. Specifically, if the mean firing rates of all neurons are scaled by a constant *G*, then the average value of the PV over the Poisson neural variability would only be scaled by *G* as follows:
which would not influence the direction of the PV.

##### Maximum a posteriori.

A general approach to decoding neural population activity optimally is to compute a maximum a posteriori (MAP) estimate of the environmental property θ using a statistical model of the neural responses (Seung and Sompolinsky, 1993; Salinas and Abbott, 1995; Ma et al., 2006). Typically, neurons are assumed to have independent Poisson responses with mean values *f _{n}*(θ). This model does not explicitly take into account the sensory cue

*S*that generates neural activities in response to the environmental property θ. Instead, the mapping from θ to

*S*is contained in the activity function

*f*(θ). In this framework, sensory reliability may be encoded by any response feature that changes the information in the population response. For example, sensory reliability has previously been hypothesized to be represented in the gain of the activity function

_{n}*f*(θ) (Pouget et al., 2003). For neurons with Poisson responses, the SD of the spike count is equal to the square root of the mean. Therefore, reliability measured as the signal-to-noise ratio of the neural response, defined as the ratio of the mean to the SD, increases with the mean response. Therefore, higher response gain corresponds to higher stimulus reliability. Specifically, for a large population of independent neurons with uniform widths, the likelihood

_{n}*p*(

**r**|θ) has a Gaussian shape the variance of which is inversely proportional to the gain of the neural responses over the population (Pouget et al., 2003) (Ma et al., 2006). Notably, sensory reliability may also be represented in the widths of the tuning functions

*f*(θ) in MAP estimates. For example, if the gain of the tuning functions is held constant, the SD of the Gaussian-shaped likelihood becomes proportional to the width of the tuning functions. Therefore, when decoding using a MAP estimate, both the gain and widths of the neural tuning functions may encode sensory reliability (Ma et al., 2006).

_{n}To decode an estimate of sound source direction θ from the model population using a MAP estimate, we compute the stimulus direction θ at binaural correlation level *BC* that maximizes the posterior probability, or equivalently the log posterior probability, of θ given the neural activities as follows:
We computed the MAP estimate of stimulus direction at different BC values from the model population response assuming Gaussian variability with mean population response **f**(θ, *BC*). The covariance matrix Σ(θ, *BC*) characterized the correlations among neurons. The MAP estimate is computed assuming that the determinant of the covariance matrix is constant across stimulus direction as follows:

#### Effect of gain and width of tuning curves on the MAP and PV decoders

We tested the effect of the gain and tuning width on the MAP and PV by separately manipulating the gain and width in the ICx model responses. The gain was held constant by keeping the maximum response of the cross-correlation model fixed at the average gain at BC = 1 or BC = 0.1 before Poisson noise was added to the response. Alternatively, the shape of the tuning curves was held constant by keeping the curves of the cross-correlation model at BC = 1 before scaling them with the appropriate gain for each BC and applying the Poisson noise to the response.

#### Measure of population information

We computed a measure of the information about the stimulus in the neural population response, known as the linear Fisher information, from the model population response to assess the quality of the code. Here, the linear Fisher information estimated the amount of information about stimulus direction available in a locally optimal linear estimator (Beck et al., 2012; Moreno-Bote et al., 2014). The linear Fisher information is inversely proportional to the square of the discrimination threshold, which is the smallest difference between two stimulus directions that can be correctly determined by an optimal decoder of the neural activity. Specifically, the linear Fisher information quantifies the information contained in the neural responses that can be extracted without further nonlinear processing.

The linear Fisher information is a function of the population response at each stimulus direction θ and is determined by the tuning curves and the covariance of the neural responses. It is given by *I*(θ) = **f′*** ^{T}*(θ)∑

^{−1}(θ)

**f′**(θ), where

**f′**(θ) is a vector of tuning curve derivatives for all neurons in the population and ∑(θ) is the covariance matrix of the population response. Because some correlations among neurons can be extremely small, small errors in estimating these correlations can result in large biases in the information estimate when inverting the covariance matrix. It has been demonstrated recently that it is possible to compute an unbiased estimate of the linear Fisher information that is given by:

*Î*(θ) =

_{bc}**f̂′**

*(θ)*

^{T}*S*

^{−1}(θ)

**f̂′**(θ)

*N*is the number of neurons,

*M*is the number of trials, and

**f̂′**(θ) and

*S*(θ) are estimates of the tuning curve derivatives and covariance matrix, respectively, at direction θ (Kanitscheider et al., 2015). The tuning curve derivatives are estimated from the responses at stimulus directions θ

^{+}= θ +

*d*θ and θ

^{−}= θ −

*d*θ as

**f̂′**(θ) =

**f̂**

^{+/−}=

**r**

*(θ) is the population response on trial*

_{i}*i*. Similarly, the covariance matrix is estimated as the average of the sample covariance matrices at the two test directions

*S*(θ) =

*S*(θ

^{+/−}) =

*N*= 500 neurons and computed responses of

*M*= 4000 trials at directions θ = 0, 20, 40, 60, and 80 degrees with dθ = 1 degree.

#### Information-limiting noise

The Fisher information in the population response can be limited by noise that cannot be distinguished from a change in the stimulus (Moreno-Bote et al., 2014; Pitkow et al., 2015). Information-limiting noise affects the linear Fisher information by causing the covariance matrix to have a component that is proportional to the product of the derivatives of the tuning curves **f′**(θ)**f′*** ^{T}*(θ), called differential correlations (Moreno-Bote et al., 2014). The covariance matrix Σ(θ) can be written as Σ(θ) = Σ

_{0}(θ) + ε

**f′**(θ)

**f′**

*(θ) where Σ*

^{T}_{0}(θ) is the covariance matrix of the noise that does not limit information and ε is the variance of the information-limiting noise (Moreno-Bote et al., 2014; Pitkow et al., 2015). The quantity 1/ε is an upper bound on the linear Fisher information created by the presence of differential correlations. This means that, no matter how the population is decoded, the information-limiting noise prevents the linear Fisher information from exceeding 1/ε. The variance of the information-limiting noise can be computed from the covariance matrix Σ(θ) and the differential correlations

**f′**(θ)

**f′**

*(θ) as ε =*

^{T}*S*(θ) and the estimate of the differential correlations

**f̂′**(θ)

**f̂′**

*(θ) as described above.*

^{T}## Results

Data from 184 neurons recorded from the ICx of two adult barn owls were used to determine the neural correlates of IPD reliability in the owl's midbrain. In this study, we treated IPD reliability as inversely proportional to the IPD variance (Landy et al., 2011). Below, we examine the relationship between the gain and width of ITD tuning, correlated variability of responses, and the IPD variability in different contexts.

### Neuronal correlates of IPD variability

For each direction of a target source, we estimated the mean IPD variability by averaging the IPD SD across frequency and directions of concurrent sources. This is a way to assess how variable the sensory cue corresponding to a given source is in the presence of simultaneous sounds from other sources. This analysis showed that the variability due to concurrent sources was lower in frontal directions than in the periphery (Fig. 2*A*). We examined how the change in IPD variability with eccentric directions is correlated with the different parameters of the population responses. We found that ITD-tuning curves broadened with the eccentricity of preferred direction (Fig. 2*B*, top, *r*^{2} = 0.54, *p* < 0.001), as reported previously (Knudsen, 1982). Remarkably, the change in ITD-tuning curve widths was highly correlated with the increase in SD of IPD (Fig. 2*C*; *r*^{2} = 0.88, *p* < 0.001). In contrast, the mean firing rate at the peak of the tuning curves measured at a constant sound level at the ears (10–20 dB above each neuron's threshold) did not vary across space (Fig. 2*B*, bottom; *r*^{2} = 0.03, *p* = 0.1).

The reliability of IPD can also be quantitatively manipulated by changing BC (Saberi et al., 1998). Changing BC varies the similarity of the sounds reaching the left and right ears by adding independent noise to each side, thus instantaneously changing the signal-to-noise ratio (Saberi et al., 1998). IPD is most reliable when the sounds are perfectly correlated (BC = 1) and least reliable when the sounds are uncorrelated (BC = 0). We found that the relationship between BC and the resulting variability of IPD was nonlinear. The SD of IPD, computed from the cross-correlation of the left and right ear inputs, was well fitted by an exponentially decreasing function of BC (Fischer and Peña, 2011). Specifically, the IPD variability was low when the BC of the signals at the two ears is between 1 and 0.5 and steeply increased for BCs <0.5 (Fig. 3*A*). Therefore, a neural response parameter encoding IPD variability should change greatly at low BC values and less at higher BC values.

To measure the effect of changing stimulus reliability on the width and the gain of ITD-tuning curves, we presented stimuli with BCs between 0.1 and 1 while recording responses of ICx neurons (*n* = 46). We found that a high level of noise in the stimulus broadened the ITD-tuning curves (Fig. 3*B*) and decreased their gain (Fig. 3*C*). Although lowering BC affected both the width and the gain, they each varied with BC in different ways (Fig. 3*D*,*E*). The width of the ITD-tuning curves was not significantly affected by low levels of noise (BC >0.4; Wilcoxon rank-sum test, *p* > 0.05), but markedly increased at higher noise levels (BC <0.4; Wilcoxon rank-sum test, *p* < 0.002). In contrast, decreasing BC from 1 to 0.9 already reduced the firing rate significantly (Wilcoxon rank-sum test, *p* < 0.05).

If the IPD SD, which varies exponentially with BC (Fig. 3*A*), is linearly encoded in the width of the tuning curves, then the width of the curves should have an exponential dependence on BC as well (Fig. 3*A*). In fact, the broadening of the curves with BC was better fit by an exponential than by a linear fit (Fig. 3*D*); the root mean square error (RMSE) of the exponential fit (RMSE = 2.8%) was considerably smaller than for the linear fit (RMSE = 11.6%). If the gain is the response parameter that captures the IPD variability, the IPD SD should be proportional to inverse square root of the gain (see Materials and Methods). However, the relationship between the inverse square root of the gain and BC was better fit by a linear model than by an exponential (Fig. 3*E*; exponential fit RMSE = 10.3%, linear fit RMSE = 5.5%). Therefore, the change in IPD variability due to external noise was better matched by the change in the tuning width than by the change in gain. We next investigate the mechanism by which this happens.

### Mechanism for the change in neural responses with IPD reliability

We first considered the mechanism controlling the change in the width of tuning curves as a function of preferred direction. Spatial tuning in ICx neurons emerges by integrating inputs from ITD-sensitive neurons across frequency channels (Takahashi and Konishi, 1986). Before frequency channels converge in ICx, upstream ITD-sensitive neurons are narrowly tuned to frequency. In these upstream neurons, the shape of the ITD-tuning curves is determined by their preferred frequency (Wagner et al., 2002). Specifically, the width of the peaks of the ITD-tuning curves of upstream neurons depends on the period of the stimulating frequency, becoming narrower as frequency increases. We tested whether, in ICx, the width of the tuning-curve peaks from stimulation with tones also decreased as the frequency increased. We observed that this phenomenon still stands in ICx (Fig. 4*A*), contrasting with a previous report that omitted tones <3 kHz (Pérez and Peña, 2006). We have shown previously that preferred direction and preferred frequency are correlated in ICx such that neurons selective for directions in the frontal space are tuned to high frequencies and neurons selective for directions in the peripheral space are tuned to low frequencies (Cazettes et al., 2014). Therefore, we hypothesized that the change in width of the spatial tuning across the ICx population is determined by the integration across frequency (Fig. 4*B*).

We used a model similar to Fischer et al. (2009) to determine whether integration over different frequency ranges could explain the broadening of ITD-tuning width in the periphery. This model described frequency integration in ICx, starting by cross-correlating time-delayed sine waves to produce simulated output of coincidence detector neurons (Carr and Konishi, 1990; Fischer et al., 2009). The output of the cross-correlation was then transformed to a spiking response by a nonlinear input–output function, which yielded a bank of modeled tonal ITD curves corresponding to neurons that are presynaptic to ICx. The modeled tonal ITD curves were then linearly combined using weights estimated from the experimental frequency tuning curves on a neuron-by-neuron basis. Finally, the linear combination of narrowband inputs was passed through an exponential nonlinearity to simulate the responses in ICx. The response of ICx neurons was therefore given by the following: *y _{i}* =

*a*exp(

_{i}*y*is the response of modeled ICx neurons,

*x*are the modeled tonal ITD curves, and

_{n}*f*their weights given by the normalized frequency tuning of ICx neurons. The parameter

_{i,n}*a*

_{i}was adjusted for each neuron to minimize the squared error between the modeled and the experimental ITD-tuning curve. The model predicted the shape of the ITD-tuning curves well (Fig. 4

*C*,

*D*). In particular, the frequency tuning of the model ICx neurons predicted the width of the ITD-tuning (

*r*

^{2}= 0.95,

*p*< 0.01; Fig. 4

*C*,

*D*). When the frequency weights

*f*

_{i,n}were shuffled between neurons, the model failed to reproduce the widths even if the parameter

*a*

_{i}was readjusted (

*r*

^{2}<0.01,

*p*= 0.34; Fig. 4

*E*). Therefore, a weighting of the frequencies in the band driving ICx neurons is sufficient to explain the width of their ITD-tuning curves. These results show that, by its effect on the shape of tuning curves, frequency selectivity is sufficient to induce a progressive change in width that correlates with the spatial dependence of the IPD variability due to the presence of concurrent sources.

Next, we used the cross-correlation model to predict the changes in tuning width and gain when IPD reliability was manipulated by changing BC. Similarly to the experimental data, the model predicted that the width of ITD-tuning curves significantly increased only after BC dropped below 0.4 (Fig. 4*F*; Wilcoxon rank-sum test, *p* < 0.05), whereas the gain already decreased significantly for BCs <0.9 (Fig. 4*G*; Wilcoxon rank-sum test, *p* < 0.001). The model reproduced both the exponential broadening of ITD-tuning (*p* < 0.001, RMSE = 3.5%) and the linear increase in the inverse square root of the gain (*p* < 0.001, RMSE = 4.5%) induced by decreasing BC. These results suggest that cross-correlation and integration across frequency explain the variation in tuning width corresponding to both the dependency of cue reliability on horizontal directions and the dynamic fluctuations in cue reliability induced by changes in signal-to-noise ratio. In the model, the cross-correlation at low BC resulted in the activity spreading over more neurons in the population on a single-trial basis rather than a shift of the peak of the population activity across trials (Fig. 4*H*). This shows that the broadening of individual tuning curves is consistent with a decrease in selectivity at the population level. Finally, we calculated the linear Fisher information in the modeled population response and characterized the correlated variability of ICx neurons as reliability changes. Using the cross-correlation model allowed us to model realistic correlated variability between ICx neurons due to external shared noise and internal fluctuations. We found that the linear Fisher information declined at peripheral directions and at lower BC (Fig. 5*A*), showing, respectively, that the reliability of the ICx representation decreases with eccentricity and with BC. Therefore, the linear Fisher information was correlated with the IPD reliability across azimuth (*r*^{2} = 0.79, *p* = 0.04) and BC (*r*^{2} = 0.62, *p* = 0.007). We also computed the linear Fisher information in a population with a uniform distribution of preferred directions and found a similar decrease in linear Fisher information (Fig. 5*B*) and thus a similar correlation with IPD reliability across azimuth (*r*^{2} = 0.92, *p* = 0.01) and BC (*r*^{2} = 0.48, *p* = 0.02). Because the gain is constant across preferred directions, this shows that the direction-dependent change in tuning width is the primary factor causing the linear Fisher information to decrease in the periphery and to correlate with IPD reliability.

We also found that the differential correlations, which occur because of higher levels of shared noise, increased at eccentric directions and as BC decreases (Fig. 5*C*). However, the differential correlations for peripheral stimulus directions computed in a model population with a uniform distribution of preferred directions (Fig. 5*D*) started increasing at a lower BC than for the nonuniform population. This is not surprising because the effect of differential correlations depends on the number of neurons in the population (Moreno-Bote et al., 2014). Here, the nonuniform population contained fewer peripheral neurons than the uniform population. However, the change in the variance of the information-limiting noise (ε) computed from the uniform and nonuniform populations had strongly correlated patterns across azimuth (*r*^{2} = 0.98, *p* = 0.002) and BC (*r*^{2} = 0.68, *p* = 0.003). This shows that the variation in the differential correlations across space largely depends on the shape of the tuning curves. The increase in differential correlations when curves broaden predicts that behavioral performance must decrease at eccentric directions and at low BC (Pitkow et al., 2015), which we examined next.

### Decoding the owl's behavior

To test which parameters of the neural responses influence the behavioral performance, we tested how different decoders predict the owl's behavioral responses at different azimuths and while varying BC. Specifically, we tested the suboptimal PV (Fischer and Peña, 2011) and the optimal MAP estimator (Seung and Sompolinsky, 1993; Salinas and Abbott, 1995; Ma et al., 2006). To test these decoders, we generated a large population of model neurons (*n* = 500) with preferred directions covering the frontal hemisphere. The population of neurons had greater density of preferred directions near the center of gaze, as seen in the owl's midbrain (Knudsen, 1982; Fischer and Peña, 2011). The cross-correlation model produced spatial tuning curves with widths that increased with eccentricity and displayed constant maximum firing rates, as observed in the data (Fig. 2*B*). In addition, the model introduces realistic correlated variability between neurons (Fig. 5).

We compared the owl's behavior with the directions predicted by the PV and MAP decoders. We found that the PV predicted well the systematic underestimation of the owl's at eccentric directions (Knudsen et al., 1979) (Fig. 6*A*, correlation of PV and behavioral accuracy: *r*^{2} = 0.67, *p* = 0.01; correlation of PV and behavioral precision: *r*^{2} = 0.87, *p* < 10^{−3}). The MAP estimates were slightly more accurate and precise than the behavior and the PV (Fig. 6*A*, correlation of MAP and behavioral accuracy: *r*^{2} = 0.63, *p* = 0.02; correlation of MAP and behavioral precision: *r*^{2} = 0.68, *p* = 0.01).

Saberi et al. (1998) previously characterized the localization behavior of the owl for different BCs by measuring the head-turn accuracy in response to four different sound locations. They showed that the owl's direction estimates remained relatively constant for BCs from 1 down to 0.4 or 0.3. However, as BC declined further, the mean direction estimates became biased toward the front (Fig. 6*B*). Both the PV and MAP estimates matched the owl's behavior as BC changed (Fig. 6*B–D*, correlation of PV and behavioral accuracy: *r*^{2} = 0.85, *p* = 0.01; correlation of PV and behavioral precision: *r*^{2} = 0.74, *p* = 0.04; correlation of MAP and behavioral accuracy: *r*^{2} = 0.88, *p* = 0.008; correlation of MAP and behavioral precision: *r*^{2} = 0.75, *p* = 0.04). Therefore, both the behavior and the PV approached optimality. In addition, we have shown that correlations among neurons emerged when the reliability of IPD was degraded. Because the PV is based solely on the signal strength of individual neurons without removing the correlated fluctuations present in the population, this shows that the correlated variability that does not limit information does not play a major role in the representation of IPD reliability. Next, we investigated why the suboptimal PV approximated the optimal MAP.

The MAP is an estimate from the posterior distribution over directions, given the population response *p*_{Θ|R}(θ|**r**). If the conditions are right in the nonuniform population code (Fischer and Peña, 2011; Rich et al., 2015), the PV can approximate a Bayesian estimate based on the posterior *p*_{Θ|S}(θ|*S*), where *S* is the sensory cues used; for example, ITD or the cross-correlation response. The PV and MAP will be similar if these two posteriors are similar. The relationship between the posteriors is defined by *p*_{Θ|R}(θ|**r**) = ∫*p*_{Θ|S}(θ|*S*)*p*_{S|R}(*S*|**r**)*dS*. There will be equality, *p*_{Θ|R}(θ|**r**) = *p*_{Θ|S}(θ|*S*), under two conditions. The first is that the mapping between the sensory cue *S* and the response *r* is unique so that *p*_{S|R}(*S*|**r**) is a delta function; that is, the sensory cues *S* can be uniquely identified from the pattern of activity across the population. This happens when there are enough neurons firing enough spikes so that the Poisson noise does not cause the patterns of activity for different stimuli to be indistinguishable. We have found that, in a population of 500 ICx neurons each firing an average of at least 2 spikes, the likelihood *p*_{S|R}(*S*|**r**) is approximately a delta function at BCs from 1 down to 0.3 (Fig. 6*E*, top and middle). The second condition for the equality of the posteriors is when the sensory cue is not informative about the stimulus and both posteriors are uniform distributions. We found that this is approximately the case at BCs near 0, where the cross-correlation response is almost flat for any target direction (Fig. 6*E*, bottom).

To further parse out the key parameters that influence behavior, we artificially kept the gain or the shape of the tuning curves constant across BC. Because other stimulus features, such as sound level, can affect the neuronal gain, we tested the decoders both at high and low gain (similar to the gains at BC = 1 and BC = 0.1, respectively). Preventing a change in gain did not significantly affect the predictions of the two decoders (Fig. 6*F*). This means that the width change with BC was sufficient to code for IPD reliability. However, preventing the broadening of the tuning curves and also eliminating the differential correlations (Moreno-Bote et al., 2014) greatly reduced the predictive power of both decoders (Fig. 6*F*). In fact, the decoders' estimates did not change with BC at any direction (Fig. 6*G*), being markedly inconsistent with the owl's behavior. Therefore, the broadening of the tuning curves and the corresponding increase in differential correlations were also necessary to code for reliability.

## Discussion

We have identified a behaviorally relevant representation of the reliability of a sound localization cue in the owl's midbrain. We showed that the mechanism that accounts for spatial cue selectivity also explained how neural responses change with degraded signals, allowing for the cue reliability to be encoded in parallel with the value of the cue. We demonstrated that the owl's behavior as reliability changes reflects a near optimal decoding of the midbrain population. The effect of signal degradation on the shape of the tuning curves was necessary and sufficient to predict the owl's behavior from stochastically firing neurons. This finding demonstrates that neurons' selectivity can act as a coding dimension independent from the gain.

The auditory space map in the owl's midbrain is a classic example of a sensory representation of the environment relevant for behavior. This study demonstrates that this map encodes, not only space, but also the reliability of spatial cues, along with the selectivity of the stimulus itself. Broadening tuning curves across a population has been observed in other sensory (Johansson and Vallbo, 1980; Wilson and Sherman, 1976) and high-order cortical systems (Harvey et al., 2013), where it was linked to a decrease in sensory precision. In the owl, the dependence of the spatial tuning on frequency tuning is crucial for the width to correlate with signal degradation. In the frontal space, where cues are the most reliable, high frequencies dominate and provide the ICx neurons with narrow spatial tuning. In the periphery, where cues are inherently less reliable, low frequencies dominate and yield broader spatial tuning. Interestingly, the auditory system is not the only sensory system where neurons' selectivity is correlated with frequency tuning. In fact, in the primary visual cortex of the macaque, neurons with narrow orientation tuning prefer high spatial frequencies and neurons with broader orientation tuning prefer lower spatial frequencies (Xing et al., 2004; Zhu et al., 2010). Remarkably, correlation in natural images, a proxy for image variability (Field, 1987), is also tightly related to spatial frequency. Robust invariance exists across natural images at high spatial frequencies due to the stronger correlation of nearby points (Field, 1987; Simoncelli and Olshausen, 2001). Therefore, a representation of stimulus reliability in the neurons' selectivity, achieved by regulating frequency tuning, could also take place in the visual system.

Our results may establish links between previous theories for the estimation of cue reliability by populations of neurons. Probabilistic population codes assume that information about cue reliability is carried in the neural variability of stochastic populations (Ma et al., 2006; Beck et al., 2008), whereas nonuniform population codes hypothesize that sensory reliability is represented in the shape of tuning curves (Fischer and Peña, 2011). Here, we show that both schemes converge because changes in neurons' selectivity by degraded signals and the corresponding changes in differential correlations (Moreno-Bote et al., 2014) are the parameters linking neural noise to cue reliability and ultimately driving performance (Pitkow et al., 2015).

We demonstrated that the midbrain population is decoded near optimally to guide behavior. Moreover, we showed that the suboptimal PV approximated the performance of an optimal decoder. Unlike the optimal decoder (Ma et al., 2006; Beck et al., 2008), a PV uses a low-dimensional probability distribution over sound direction. This probability distribution depends only on the statistical relationship between the sound source direction and the localization cue IPD and is therefore more plausible for learning and implementation than the high-dimensional probability distribution used in the MAP framework. Furthermore, the PV discards the information contained in the response gain from the decoding. Because the gain is not only correlated with stimulus reliability, but also varies with sound level, removing the gain from the decoding of sound direction could eliminate this confound. This would allow “where” to be decoded independently from “what” the stimulus is. However, this simplified code may result in some accuracy errors that would not occur if stimulus information was taken into account in the decoding, as it is in the MAP framework.

Our findings lend support to Bayesian theories of brain processing (Knill and Pouget, 2004; Ma et al., 2006; Fischer and Peña, 2011). Bayesian inference allows sensory information and prior knowledge to be combined optimally. In the Bayesian framework, a likelihood function captures the reliability of the sensory input and a prior distribution represents previous knowledge. The combination of the prior and likelihood are used to produce an estimate of the variable of interest that reflects a probabilistic model of the environment. Typically, the noisier the sensory inputs, the more the final decision will rely on prior knowledge because, at high noise levels, the likelihood function widens (Knill and Pouget, 2004; Fischer and Peña, 2011). A computational model of the barn owl's auditory midbrain proposed that the sound localization behavior is a Bayesian inference process (Fischer and Peña, 2011). In ICx, information about the stimulus reliability combined with prior information could allow for Bayesian inference (Fischer and Peña, 2011). However, until now, this coding had not been verified or the mechanism of its emergence understood. By showing that the shape of tuning curves explains how reliability is incorporated into the code and that the regulation of frequency tuning is the underlying mechanism, our study demonstrates the plausibility of an additional framework for how the brain can support perceptual judgments operating in situations of uncertainty.

## Footnotes

This work was funded by the National Institutes of Health (Grant DC012949). We thank Kourosh Saberi for providing the data used in Figure 6

*B*, Clifford Keller for providing data and for analysis of the head-related impulse responses, Ruben Coen-Cagli for helpful discussion, and Adam Kohn and Odelia Schwartz for comments on early versions of the manuscript.The authors declare no competing financial interests.

- Correspondence should be addressed to Fanny Cazettes, Department of Neuroscience, Albert Einstein College of Medicine, 1410 Pelham Parkway South, Bronx, NY 10461. fanny.cazettes{at}phd.einstein.yu.edu