\documentclass[a4paper,11pt]{article}

% Use \includegraphics{file} if figures are needed
\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage[sumlimits,intlimits]{amsmath}
\usepackage{amsfonts,amssymb}

\oddsidemargin 5mm
\evensidemargin -5mm
\topmargin -15mm
\textwidth 160mm
\textheight 24.5cm

\newcommand{\df}{\ \textup{d}}

\begin{document}

\begin{center}
\textbf{{\Large 
On the precision of neural computation with interaural level differences in the lateral superior olive
}}

\bigskip

\textit{Zbynek Bures} \\
College of Polytechnics \\
Tolsteho 16, Jihlava, Czech Republic \\
buresz@vspj.cz 
\end{center}

keywords{binaural hearing, ideal observer, subtraction of firing rates, coincidence detection, neuronal arithmetic}

\section*{Abstract}

Interaural level difference (ILD) is one of basic binaural clues supporting the spatial localization of a sound source. Due to acoustic shadow cast by the head, a deviation of a sound source from the medial plane results in an increased sound level at the nearer ear and a decreased sound level at the averted ear. In mammalian auditory brainstem, the ILD is evaluated by a specialized neuronal circuit, comprising particularly the binaural neurons in the lateral superior olive (LSO) which receive excitatory projections from the ipsilateral side and inhibitory projections from the contralateral side. As the sound level is encoded predominantly by neuronal discharge rate, the principal function of LSO neurons is to estimate and encode the difference between discharge rates of the excitatory and inhibitory input. However, the exact mechanism of this operation is not known; furthermore, it is not clear how the evaluator attains the remarkable precision of ILD assessment observed experimentally. We employ a probabilistic model and complementary computer simulations to explore how various properties of the system affect the just noticeable differences (JND) of ILD. Introducing a concept of an ideal observer, we assess the lowest possible JNDs of ILD depending on the statistics of interacting spike trains, overall firing rate, detection time, number of converging fibers, and finally on the evaluation mechanism itself. The results show that the JNDs of ILD strongly rely on the precision of spike timing, however, with an appropriate parameter setting, the lowest theoretical values are similar or better than the experimental values. Furthermore, the mechanism based on excitatory-inhibitory detection of coincidence may give better results than ideal subtraction of firing rates.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Introduction}

Mammals, including humans, use the difference in sound levels at the ears (interaural level difference, ILD) as a clue for the spatial localization of the sound source \cite{Mills58, Blauert}. The detection and coding of the ILD are supported by the cells in the lateral superior olive (LSO) that scale their response magnitude according to the ILD \cite{Tollin2003}. The LSO neurons are innervated from both ipsi- and contralateral loci, receiving information on sound intensities in both ears encoded predominantly by means of the instantaneous firing rate \cite{goldberg_brown, Tollin2003}; their task is to evaluate the input spike rates and translate their difference (induced by the ILD) to the output discharge rate. 

The LSO cells receive ascending excitatory inputs from the ispilateral cochlear nucleus (CN) and inhibitory inputs coming from the contralateral CN via the medial nucleus of the trapezoid body \cite{Pickles}. The mean number of excitatory and inhibitory afferents terminating on an LSO neuron is 9.6 and 8.2, respectively \cite{Sanes90}. Projections to LSO neurons may originate either in bushy cells in the ipsi- and contralateral CN \cite{cant86, smith91, smith93} or in planar multipolar cells in the ipsilateral CN \cite{cant86, doucet2003}. The former pathway preserves precisely the timing of spikes in the auditory nerve (AN) -- the inputs of this origin (both excitatory and inhibitory) are thus nearly identical to the AN activity. The major domain of LSO operation lies in the high frequency band, where the AN activity does not phase-lock to the stimulus waveform and the spike trains resemble rather a Poisson random process \cite{young_barta}. On the contrary, the planar cells (which serve as another possible source of excitatory spikes) exhibit a chopping response only weakly related to the instantaneous AN activity \cite{Palmer2003}. For these reasons, two cases of LSO inputs are analyzed: random primary-like activity based on the dead-time Poisson process and the nearly-periodic (chopping) behavior of planar cells. 

Due to the variability of neuronal firing, the correct functioning and precision of the ILD evaluation is affected by the stochastic properties of the input spike trains \cite{Bures2012}. Even if an ideal subtraction of firing rates is performed, the stochastic nature of inputs brings considerable uncertainty into the ILD evaluation. On the other hand, it is known from both psychophysical and neurophysiological studies that the ILD may be assessed with a remarkable precision, the just noticeable differences (JND) ranging from 0.5 dB to 4 dB \cite{Tollin2008}. The purpose of the current work is to bridge the gap between the assumed neural mechanism and psychophysics by assessing the theoretical limits of the precision of ILD evaluation depending on various parameters of the system and by relating the outcomes to the psychophysically observed values. Two models of evaluation mechanism are considered and compared: ideal subtraction of firing rates (SFR) in a chosen counting window, and coincidence detection (CD) mechanism with an ideal counter on its output.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Methods}

\subsection{\label{Meth_JND}Just noticeable differences of ILD}

A higher variability of firing leads to a lower precision of the rate code. Intuitively, if a repeated presentation of the same stimulus evokes each time a different spike count, then to distinguish between two different stimuli, the associated spike count change must be larger than the spike count variability. Mathematically, the current work determines the JND of ILD (in other words, the precision of ILD coding) using theoretical tools derived in \cite{Green1966, Teich1979, Colburn2003}. We will ask whether it is possible to distinguish between two random processes with rates $L_1$ and $L_2$, $L_2 > L_1$. If we count events in a given counting window, we get counts $n_1$ and $n_2$. The probability that the observer obtains a result that $L_2 > L_1$ equals to the probability that $n_2 > n_1$. Let's assume that the random variables $n_i$, $i=1, 2$, have probability distribution $p(n_i)$ with means $\mu_i$ and roughly equal standard deviation $\sigma$. A detection distance is then defined as 
%
\begin{equation}
d' = \frac{\mu_2 - \mu_1}{\sigma}
\label{R_SDT1}
\end{equation}
%
This definition expresses the fact that the larger the variance of the spike count, the worse the detection capability. 

In psychophysics, a threshold value is commonly defined as that value for which the percentage of correct answers equals 75. In our case, the examined value is the just-noticeable change of firing rate, $\Delta L = L_2 - L_1$. Assuming that both $p(n_1)$ and $p(n_2)$ are Gaussian, the 75\% probability of $n_2 - n_1 > 0$ corresponds to $d' = 1$. To obtain the JND of firing rate, we scale the detection distance with $\Delta L$ and put $\delta' = d'/\Delta L$. Then, the JND of firing rate is 
%
\begin{equation}
\Delta L_\text{JND} = 1/\delta' = \frac{L_2 - L_1}{\mu_2 - \mu_1}\sigma
\label{R_SDT2}
\end{equation}
%

The JND of ILD may be evaluated, e.g., by fixing the contralateral sound level and varying the ipsilateral sound level. However, the relationship between the sound level at the ear and the discharge rate at the LSO input depends on many factors (e.g., auditory nerve rate-level function, RLF). For this reason, we evaluate primarily the just noticeable changes of firing rate of the excitatory LSO input. The corresponding threshold ILDs are estimated by considering that the maximum slope of auditory nerve RLFs is on average 5 spikes/sec/dB \cite{young_barta}, that means that an incerase in sound level equal to 1~dB increases mean firing rate by 5 spikes/sec. 


\subsection{Model of an LSO cell}

We introduce two very simple abstractions of an LSO unit that represent the two ILD-evaluating mechanisms considered. Similarly as in \cite{Bures2012}, the cell receives one excitatory and one inhibitory input and its output results from the interaction of both inputs. The inputs and the output are modeled as one-dimensional random point processes. For the mathematical analysis, it is not necessary that the events of these processes be action potentials in the physiological sense; rather, the events may be regarded as a temporary increase (excitatory) or decrease (inhibitory) of the probability of an output discharge. The detailed biophysical mechanisms involved in real neural systems are neglected in order to concentrate on the combinatorial and probabilistic phenomena, such an approach allows to lower the number of free variables in the model so that the influence of the individual parameters is more understandable. Nevertheless, the principles explored here are inevitably embedded in real complex systems; our modeling results set the theoretical bounds for the precision of ILD coding for the extreme case of ideal observer. Therefore we are also not interested in the exact timing of individual output events, we only derive the mean and variance of the output event count in a given counting interval. 

As the first model, a cell performing an ideal subtraction of firing rates is introduced. If $n_\text{E}$ and $n_\text{I}$ are the event counts of the excitatory and inhibitory inputs, respectively, in a given counting window, then the output event-count of this cell equals $n_\text{SFR} = \textrm{max}(0,n_\text{E}-n_\text{I})$. Depending on the parameters of the model, the just-noticeable changes of ILD are assessed according to Section \ref{Meth_JND}.

As the second model, a cell performing the coincidence detection is introduced. The details of this model may be found in \cite{Bures2012}; in principle, an excitatory input event may be cancelled by an inhibitory input event on condition that the inhibitory event arrives simultaneously with, or within a short time interval $\tau$ before, the excitatory event. The inhibitory effects of more than one inhibitory event do not summate. The output of this cell is given by the uncancelled excitatory events, these events are counted in a chosen counting window and the mean and variance of the count is used for the purposes of the assessment of JND of the ILD. 

The parameters of the LSO model that are varied during the evaluation are: mechanism of operation (SFR or CD), length of the counting window (detection time) $\Delta$, and the coincidence interval $\tau$ (only for the CD mechanism).


\subsection{Stochastic inputs to an LSO cell}

Each afferent LSO input (one excitatory and one inhibitory) is modeled as a random point process. For reasons explained in the Introduction, we employ two basic classes of inputs: a primary-like input modeled as a dead-time Poisson process with shifted exponential distribution of inter-spike intervals (ISI), and a chopper input modeled by a process with a uniform distribution of ISI. Only stationary asynchronous firing with no phase-lock to the stimulus waveform is considered.

A key parameter which profoundly influences the JND of ILD is the variability of inputs. The temporal variability of neuronal discharge patterns may be expressed using the coefficient of variation (CV: standard deviation of ISI divided by mean ISI) or Fano factor (FF: spike count variance divided by mean spike count over some counting interval). The free parameters of the input processes that are varied during the evaluation are: ISI distribution (shifted exponential or uniform), CV, and firing rate (for the excitatory input only).

Besides the inputs consisting of a single random point process, we also considered a possibility of convergence of multiple excitatory and inhibitory fibers at one cell. The corresponding model input is given by a sum of multiple independent random renewal processes with identical statistics. The properties of the individual processes are set analogously to the case of a single point process: the ISI distribution, CV, and firing rate are fixed to appropriate values. The free parameter which is varied in this case is only the number of converging fibers, $c$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Results}

Unless stated otherwise, the parameters used in the simulations were set to the following values: 

\begin{itemize}
\itemsep0em
\item excitatory firing rate: 200 events per second
\item inhibitory firing rate: 100 events per second
\item CV of excitatory fiber: 0.25
\item CV of inhibitory fiber: 0.5
\item counting window $\Delta$: 200 ms
\item coincidence interval $\tau$: 6 ms
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Each input is a single renewal process}

To explore the relationship between the JND of excitatory rate and various model parameters, we start with the most simple case of both inputs being described by a single random renewal process with a given distribution $f(t)$ of ISI ($f_\text{E}$ and $f_\text{I}$ for the excitatory and inhibitory process, respectively). The results were obtained using numerical simulations in MATLAB environment.

In the case of ideal subtraction of firing rates (SFR), the results show that the JND increases nearly linearly with increasing CV of the excitatory or inhibitory input (see Fig. \ref{F_1}A), the distribution of ISI of the inputs does not play a role. 

Fig. \ref{F_1}B shows the dependence of the JND on $\Delta$. The length of the counting window $\Delta$ strongly influences the JND -- the longer the $\Delta$, the smaller the JND. Explicitly, $\text{JND} = K/\sqrt{\Delta}$, where $K$ is some positive constant. 

In the next experiment, the excitatory fiiring rate was varied between 100 and 2000 events per second; the inhibitory firing rate was always set to one half of the excitatory rate. In this case, we have found a dependence inverse to the influence of $\Delta$: given a constant CV of the inputs, the JND grows with the square root of the input rates (see Fig. \ref{F_1}C). This leads to a paradox that to keep the JND constant, higher input rates must be accompanied by \emph{longer} counting windows.
%
\begin{figure}
\centering
\resizebox*{1\textwidth}{!}{\includegraphics{Fig1}}
\caption{JND depending of the excitatory CV (A), counting window $\Delta$ (B), and on excitatory firing rate (C).}
\label{F_1}
\end{figure} 
%

The second model of ILD evaluation, excitatory-inhibitory coincidence detection (CD), gives results very similar to the SFR model. For the CD mechanism, however, it is necessary to ensure that the algorithm really performs subtraction -- see \cite{Bures2012}. It is noteworthy that the obtained JNDs may be smaller than those resulting from ideal SFR of the same inputs. An important parameter is the coincidence window $\tau$: the longer the $\tau$, the smaller the JNDs, see Fig. \ref{F_2}. For comparison, the JNDs obtained by SFR are shown as well -- the line is nearly constant as the coincidence window $\tau$ does not affect the SFR mechanism. The figure does not show JNDs for $\tau<5$, as for the firing rates used (excitatory rate 200 events per second, inhibitory rate 100 events per second), the CD mechanism performs subtraction only when $\tau \ge 5$ \cite{Bures2012}.
%
\begin{figure}
\centering
\resizebox*{.35\textwidth}{!}{\includegraphics{Fig2}}
\caption{Comparison of the SFR and CD algorithm: JND depending of the coincidence interval $\tau$.}
\label{F_2}
\end{figure} 
%

For the first model mechanism (ideal SFR), an exact mathematical description is possible. A crucial task is to find the mean $\mu$ and variance $\sigma^2$ of number of events $n$ occuring in a window $\Delta$, or, in a more general case, the probability distribution $p(n)$ of the event-count $n$, given the distribution of inter-event intervals of the process. The $p(n)$ may be obtained by expressing the individual probabilities $P(n=M), M=1, 2, \dots$. Let $T_1$ be the time from the beginning of the counting window to the first event occuring in that window, and let $T_i, i=2, 3, \dots$, be the time intervals between the $i$-th and $(i-1)$-th event occuring in the counting window. Then the probability that $n=M$ is 
%
\begin{equation}
p(M) = P(n=M) = P\left(\sum_{i=1}^M T_i <\Delta\right) - P\left(\sum_{i=1}^{M+1} T_i <\Delta\right)
\label{R_SFR1}
\end{equation}
%

To express the probabilities at the right side of Eq. \ref{R_SFR1}, we need to get the probability distributions of sums of inter-event intervals $T_i$. The distribution of a sum of random variables is equal to the convolution of the distributions of the individual variables. We know the distribution of $T_i$ for $i>1$ -- this is the distribution of the ISI of the renewal process, $f(t)$. However, the distribution of $T_1$ may be different -- we do not know in which instant of the process the counting window starts and unless the beginning of the window coincides with some event of the process, the distribution $f_s(t)$ of the random variable $T_1$ is not equal to $f(t)$. 

The probability distribution $f_s(t)$ depends on $f(t)$. If we start counting during the course of the process, the beginning of the counting window will lie either at the instant of some event, or somewhere between two successive events. We seek the probability distribution of the interval $T_1$ between the beginning of the counting window and the first succeeding event. This interval is certainly non-negative, hence $f_s(t)=0$ for $t<0$. At the same time, $T_1$ can not be larger than the longest possible ISI given by $f(t)$. This means that the distance $T_1$ can be equal to some $t_0$ only when the inter-event interval $T$ during which the counting started is at least the same or larger than $t_0$, $T \ge t_0$. Therefore, $f_s(t)$ is proportional to the probability that $T \ge t$:
%
\begin{equation}
f_s(t) = \frac{1}{K_1}P(T \ge t) = \frac{1}{K_1} \int_{t}^\infty f(x) \df x
\label{R_SFR2}
\end{equation}
%
where $K_1$ is set so that
%
\begin{equation}
\int_{-\infty}^\infty f_s(t) \df t = 1
\label{R_SFR3}
\end{equation}
%

We can now convolve the distribution $f_s(t)$ with an appropriate number of distributions $f(t)$, and we can rewrite Eq. \ref{R_SFR1} as
%
\begin{equation}
p(M) = \frac{1}{K_2} \int_{0}^\Delta f_s(x) * f^{(M-1)}(x) \df x - \frac{1}{K_3} \int_{0}^\Delta f_s(x) * f^{(M)}(x) \df x
\label{R_SFR4}
\end{equation}
%
where $f(x)*g(x)$ stands for convolution of functions $f$ and $g$ and the notation $f^{(M)}$ stands for $M$-fold convolution of funcion $f$ with itself. The normalization constants $K_2$ and $K_3$ are set so that
%
\begin{equation}
\frac{1}{K_2} \int_{-\infty}^\infty f_s(x) * f^{(M-1)}(x) \df x = 1
\label{R_SFR5}
\end{equation}
%
and 
%
\begin{equation}
\frac{1}{K_3} \int_{-\infty}^\infty f_s(x) * f^{(M)}(x) \df x = 1
\label{R_SFR6}
\end{equation}
%

The next step is to count the mean and variance of the probability distribution $p(n)$ and use the tools from Section \ref{Meth_JND} to evaluate the JND of firing rate and the JND of ILD. The results of the analytical model are shown along with the results obtained by numerical simulations in Fig. \ref{F_1}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Each input is a superposition of independent renewal processes}

The above results hold for excitatory and inhibitory inputs each comprising a single random renewal process with a given ISI distribution. As multiple excitatory and inhibitory fibers may converge at one LSO cell, such situation was considered by modeling each input as a sum of multiple independent random renewal processes with identical statistics. Each process then represents one of the converging fibers. 

A superposition of multiple independent random processes results in a process which has nearly an exponential distribution of inter-event intervals. Hence, a common expectation would be that the resulting process is close to Poissonian and that the JND would not improve with an increasing number of converging fibers. However, despite that the superposed processes are independent and the resulting ISI distribution resembles more and more an exponential distribution as the number of fibers increases, the JND \emph{decreases} for an increasing number of fibers, see Fig. \ref{F_3}A. This behaviour is caused by the fact that despite a nearly exponential distribution of ISI, the sum of multiple non-Poisson processes is not a Poisson process; it is not even a renewal process. The stochastic properties of this process differ for different time scales; in particular, on longer time scales, the superposed process preserves the properties of the individual component processes \cite{Lindner2006}. The simulations show that in the case of a renewal process, Fano factor is roughly proportional to the square of CV (see Fig. \ref{F_3}B). However, if we sum up multiple processes each having CV equal to, e.g., 0.25, the resulting CV quickly approaches 1, yet the resulting FF corresponds to the CV of an individual component process (i.e., $\textrm{FF}\sim0.25^2$), leading to markedly lower JNDs, see Fig. \ref{F_3}C. Furthermore, as the summary process is not renewal, the spike counts (and thus also the JNDs) obtained at the beginning of the process are different from those obtained during the course of the process (see Fig. \ref{F_3}A). This may have implications for ILD perception: JNDs of ILD may be lower when evaluated using short tone pips than when using continuous signals.
%
\begin{figure}
\centering
\resizebox*{1\textwidth}{!}{\includegraphics{Fig3}}
\caption{Convergence of multiple fibers. A: JND depending on the number of converging fibers, counting window positioned at the beginning of the process or at an arbitrary instant of the process. B: Relationship of coefficient of variation and Fano factor. C: Coefficient of variation and Fano factor for increasing number of converging fibers.}
\label{F_3}
\end{figure} 
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Discussion}


%%%% dodelano sem
DODELANO SEM, DISKUSI NUTNO JESTE NAPSAT






Concluded, the lowest theoretical JNDs of ILD are similar or better than the experimental values. However, the LSO cells hardly behave as an ideal detector; a certain worsening of JND is to be expected. Given that the lowest excitatory and inhibitory CVs are approx. 0.2 and 0.7, respectively \cite{Palmer2003, Young1986}, then to attain the desired precision, the counting window should be at least 200~ms long, ca. 10 excitatory and 10 inhibitory fibers should converge at one cell, and the input firing rates should be kept low. Furthermore, spike timing must be conveyed accurately from the auditory afferents to the LSO.

Supported by projects M00176 "Elektronicko-biomedicinska kooperace" and "Podpora a individualni rozvoj perspektivnich akademickych pracovniku na VSPJ".

\medskip\noindent\textbf{\textit{Keywords:}}
Binaural interaction, Lateral superior olive, Subtraction of firing rates.



\begin{thebibliography}{99}

\bibitem[1]{Palmer2003} Palmer, A.R., Wallace, M.N., Arnott, R.H., Shackleton, T.M. (2003)
Morphology of physiologically characterised ventral cochlear nucleus stellate cells. \emph{Exp Brain Res}, \textbf{153}: 418--426.

\bibitem[2]{Tollin2008} Tollin, D.J., Koka, K., Tsai, J.J. (2008)
Interaural level difference discrimination thresholds for single neurons in the lateral superior olive. \emph{J Neurosci}, \textbf{28}(19): 4848--4860.

%\bibitem[3]{Heinz2005} Heinz, M.G., Issa, J.B., Young, E.D. (2005)
%Auditory-nerve rate responses are inconsistent with common hypotheses for the neural correlates of loudness recruitment. \emph{JARO}, \textbf{6}: 91--105.

\bibitem[3]{Young1986} Young, E.D., Barta, P.E. (1986)
Rate responses of auditory nerve fibers to tones in noise near masked threshold. \emph{J Acoust Soc Am}, \textbf{79}(2): 426--442.

\bibitem[4]{Lindner2006} Lindner, B. (2006)
Superposition of many independent spike trains is generally not a poisson process. \emph{ Phys Rev E}, \textbf{73}: 022901.

%\bibitem[5]{Avissar2007} Avissar,M., Furman, A.C., Saunders, J.C., Parsons, T.D. (2007)
%Adaptation reduces spike-count reliability, but not spike-timing precision, of auditory nerve responses. \emph{ J Neurosci}, \textbf{27}(24): 6461--6472.


\end{thebibliography}

\end{document}
