On the precision of neural computation with interaural time differences in the medial superior olive
    
Petr Marsalek, Marek Hajny, Zbynek Bures
Institute of Pathological Physiology, First Medical Faculty, Charles University in Prague, U Nemocnice 5/478, 128 53, Praha 2, Czech Republic
Department of Radioelectronics, Faculty of Electrical Engineering, Czech Technical University in Prague, Technická 2/1902, 16627, Praha 6, Czech Republic
Czech Institute of Informatics, Robotics and Cybernetics, Czech Technical University in Prague, Jugoslávských partyzánů 4/1903, 166 36, Praha 6, Czech Republic 
College of Polytechnics, Tolstého 16/1556, 586 01, Jihlava, Czech Republic
Corresponding authors: Petr.Marsalek@LF1.CUni.CZ, Marek.Hajny@LF1.CUni.CZ 
Keywords: binaural hearing, cochlear implant, coincidence detection, ergodic hypothesis, ideal observer, interaural time difference, just noticeable difference, lateral and medial superior olive, neuronal arithmetic, psychophysics, sound localization. 
Abstract
Abstract will be rewritten after the manuscript completion. Paper should explain analytical MSO model, all figures and tables.
We study a model of mammalian sound source localization in horizontal plane. (Experiments on small rodents indicate that some mammals use broadly tuned channels of azimuth for localization in horizontal plane.) In mammals this neural computation is implemented by medial superior olive for low frequency sounds up to 1500 Hz. It has been shown previously that spike timing jitter, coincidence detection window length, sound frequency, among other input parameters, influence the output precision, measured by the just noticeable difference in output of the circuit. We use stochastic model with spiking neurons. We explore variations of the parameters mentioned above. We calculate properties of the model with the use of analytical methods. Predictions of this model have straightforward applications in testing and designing stimulation protocols used in cochlear implants.  
Abbreviations and symbols
F, , sound frequency; , critical sound frequency value; j, sound phase; ILD, interaural level difference; ITD, interaural time difference; ISI, inter-spike interval; JND, just noticeable difference; LSO, lateral superior olive; , , A, B, C, proportionality constants; l, sound level; MSO, medial superior olive; r(V), , VS, vector strength; , firing rate; s, standard deviation; , timing jitter; t, time; T, , sound period;
1  Introduction
Mammalian sound localization circuits contain two nuclei in the auditory brainstem, the medial and the lateral superior olive (MSO and LSO). Neurons in these nuclei are the first binaural neurons in the auditory pathway, they are connected to both ears.
This article presents description of information encoding and neural computation in the MSO obtained mostly with analytical computations. Using the analytical tools we extend quantitative results obtained by numerical computations in [Sanda2012]. We compare this analytical MSO description to the LSO description in [Bures2013] to arrive to unified description of neural circuits in the superior olive. We use this description to find maximum spike timing precision attained by the following two respective modalities (low and high frequency sound localization) of natural and electrical hearing.
Due to the physical nature of the binaural sound, the MSO neurons process spike timing differences in the range of tens of microseconds, ms. The MSO processes low frequency sounds, in human this is from 20 Hz to not more than 2 kHz. The LSO processes high frequency sounds, in human this is from 1 kHz up to 20 kHz. The overlapping region is known to have a sensitivity drop at 1.5 kHz [Mills1958].
Firing rate, first spike latency and individual spike timings are used in neural system coding, especially in the auditory pathway. The highest known spike timing precision that has been proven to be utilized in the mammalian auditory system is the motor response in a bat. In behavioral experiment, the big brown bat, Eptesicus fuscus catches tidbit larvae of the meal-worm beetle Tenebrio molitor. The motor precision is in the range of hundreds of nanoseconds, ns, and auditory separation of ultrasound echoes is in the order of microsecond units, 1 ms [Simmons1998]. Comparable performance of the human MSO is in the range of microseconds, 10 ms and can be improved after several hours of training.
Different neural mechanisms are employed in the two nuclei. It has been reported that computation in the MSO is independent on sound intensity [Grothe2010]. With higher sound intensity, first spike latency is shortening. Relation of this dependence to ITD and ILD has been described, yet it is difficult to interpret. [Michelet2012] have shown the latencies in experiment on the domestic cat, Felis silvestris.
The LSO and MSO extract location information with the use of different physical cues. The sensitivity of the system in dependence on the main sound frequency to sounds of different frequencies should be different.
In extant human, the MSO is the larger of the two nuclei and contains approximately 10000 - 11000 neurons and the human LSO contains 5600 neurons [Moore2000]. As loudnes is implemented simpler that microsecond time delay, most of currend auditory technology works as if the more important of the two nuclei is the LSO in Homo sapiens.
We have previously studied parameter free models of human sound localization. Current development of binaural cochlear implant technology enables to stimulate auditory periphery with rich repertoire of patterns. Better understanding to the MSO mechanisms will help to improve electrical hearing.
( Note: First look at the data enables us to generalize when we ... Speeding up, or slowing down the system time can be extended to two octaves up and down the regular time. Table ~\ref{tab:speed-up} has been deleted. )
The range of sound frequencies processed in the MSO circuit is limited to the low frequency band. In human, where audible sound frequencies range from 20 Hz to 20 kHz, this low frequency band spans from 20 Hz to 1.5 kHz. Moreover, in the marginal values of this interval, performance of sound azimuth discrimination is worse. It is generally agreed that the main reason, why the workings of the MSO circuit deteriorate towards higher frequency is lowering of the synchronization of spike trains in the circuit with the sound source phase. The synchronization between two corresponding series of point events can be expressed as a (discrete version of) vector strength, defined below in equation (?).
2  Methods
2.1  Preliminaries
Our procedures are mostly in using computational models, experimental data and statistical methods. Here we list several assumptions used in the paper. In any model of ITD based azimuth discrimination, time differences have to be computed with high accuracy (in the range of tens of microseconds).
[(1)] 
1. For single MSO neuron constants, see [Toth2015]. Membrane time constants are lower than 1 ms. 
2. For the synaptic machinery, see [Toth2015]. 
3. For the ergodic assumption see [Sanda2012], [Toth2018].  The basic set of parameters is shown in Table ?.
2.2  Model of the MSO neural circuit
Our model MSO circuit is based on connected phenomenological neurons. Input sound to left and right sides is transformed by the auditory periphery module into spike trains. Spikes in these trains are point events, only spike times matter and the details of spike numerical implementation do not make any differences in model output. These spike trains converge and diverge into higher order neurons (from the auditory nerve and cochlear nucleus through the medial and lateral nuclei of trapezoid body) up to the neurons of medial superior olive, which are first binaural neurons. Output of these neurons is the azimuth signal encoded in a spike train. Basic set of parameters is shown in Table ?.
2.3  ITD readout curve
Let us have a monotonous function with firing rate as an input, which outputs azimuth. We will call it the ITD readout curve. In the paper by [Sanda2012] this curve is constructed by curve fitting to simulated points. Here we construct the curve based on assumption that the main frequency of sound input exist, is unique and is known. In addition to this known frequency, other parameters of the readout curve are set to make the fitting well posed and to obtain correct position of the curve maximum.
2.4  Vector strength
The vector strength is defined as follows, by [Goldberg1969]. Let us have sample (spike) phases , i=1,2,…,N relative to phases of a given input master periodic function, which does not enter the formula. (Such function can be for example sound stimulus eliciting the spike train.) (Discrete sum) vector strength of sample  is defined as

		(1)

2.5  JND and ideal observer
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 different spike count, then to distinguish between two different stimuli, the associated spike count change must be larger than the spike count variability. This way we determine the rate code JND (in other words, the precision of rate coding) using theoretical tools.
We will ask whether it is possible to distinguish between two random processes with rates  and , . If we count events in a given counting window, we get counts  and . The probability that the observer obtains a result that  equals to the probability that . Let’s assume that the random variables , i=1,2, have probability distribution  with means  and equal standard deviations s. A detection distance is then defined as 

		(2)
This definition expresses the fact that the larger is the variance of the spike count, the worse is 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, . Assuming that both  and  are Gaussian distributed, the 75% probability of  corresponds to d'=1. To obtain the JND of firing rate, we scale the detection distance with DL and put d'=d'/DL. Then, the JND of firing rate is

		(3)

3  Results
We investigate how circuit output and overall performance depend on sound frequency and sound intensity. All Figures show sample data of numerical simulation as black line, lower estimate of data as blue line and upper estimate of data as red line. Figure ? shows an example of typical vector strength (which cannot be larger than 1) lowering towards higher frequencies, as it can be observed in the numerical module of the auditory periphery (auditory nerve and cochlear nucleus). The prevailing majority of neurons in the auditory pathway has vector strength spike train statistics sigmoidally dropping towards higher sound frequencies as it is in this example. In this figure, simulated data (Note: recorded by [Joris1996] from the MSO in cat.) were fitted to the sigmoidal curve with the general formula of the Boltzmann function used in [Marsalek2005]. The curve fit of vector strength  in dependence on sound frequency  is:

		(4)
where  (sound coefficient) and  (critical coefficient),  (critical half frequency) are proportionality constants.
Figure ? shows the curves limiting the ITD obtained with the basic parameter set in dependence on the sound frequency.
The quadratic curve fit of the JND denoted  is:

		(5)

where upper and lower bound parameters A and B are fitted constants and the other parameters and variables are as above. See also [Zwislocki1956].
Figure ? shows how the (JND) of ITD depends on timing jitter . There is one critical  and one optimal value (normalized to ) of timing jitter related to the output JND and firing rate. Simulations show that with lowering timing jitter the circuit output is virtually more and more precise. Yet, when the jitter is lower than critical value , determined by intrinsic noise, duration of coincidence detection window, and by other time constants, the precision lowers again. The two curves fitted to the simulation are:
1. fit of exponential function to simulations (red curve): 

		(6)
where ,  and  are fitted parameters. This is shown conveniently by the logarithmic y-axis in this figure.
2. Another fit, which also takes into account shot noise in lower jitter values, is to a quadratic function (blue curve):

		(7)

where , , . There is only one parameter (A) fitted to data, as the point  has been chosen to be a unit, this fit is a normalized fit of the model. (Note: this must be explained - that this fit is after normalization! ).
Logarithm of the simulated JND lowers with the exponential curve (?), which is concave function of jitter  as the jitter gets lower. The trend towards higher accuracy diverges from the parabolic fit (?), when jitter reaches critical value between  and .

		(8)

Beyond that point towards the lower jitter values, the neural circuit cannot function properly, as too low jitter prevents the interaction of spikes from the left and right side within the coincidence detection mechanism.
(Note: The other (optimal) jitter value notable in this Figure ? is firing rate in dependency on the jitter values.)
(Note: Check this.) This is an analytical dependency obtained in [Salinas2000] for a perfect integrator model with several inputs. The mechanism studied thereof is close to the MSO neural mechanism studied here.
The firing rate changes in dependency on the input spike timing variability of partially correlated definition? / description?  input spike trains.
Figure ? shows the ITD readout curve. The rising slope of this curve is used as a readout function yielding the firing rate in dependency on the ITD, which in turn signals the sound azimuth to the next nuclei of the auditory pathway.

4  Discussion
4.1  Things to do
(First/ we may skip this, since this adds no explanatory value to the text: Reproduce the curve in panel with sound frequency on x-axis and JND on the y-axis in Fig. ?.) 
Second: Present the curve from the right panel of Fig.4 from Sanda P., Jitter Effect on the Performance..., as part of Fig. ?. 
Some other tentative questions: What is the highest slope of the ITD interpolation curve, such that it gives the resolution of 4 (angular degrees) in the ITD = 0 direction?  (This we will do later: We should also discuss the relation of the MSO neural circuit to its neighboring nucleus, the LSO, and its output.)
Other periodic functions can be used as the ITD readout curves. In [Toth2018] we have compared the Sine function with the circular Beta density, and as other alternatives, regular Normal density and circular Normal density have been used and tested in previous versions of our model. No differences between these densities with proper parameters have been shown by common testing differences between two probability densities, ibidem.
(Note: Below are suggestions what papers to cite in Discussion. )
Cite: [Mlynarski2014], this is a paper about statistics of natural environmental sounds in relation to low frequency sound localization...
Cite: [Bures2019], this is a paper submitted this spring to Biosystems, about encoding of low frequency sounds...

(Note: Should we have here a schema showing the anatomical wiring of the mammalian brainstem circuit together with the equivalent circuit used in the model. Without the loss of generality, for the sake of simplicity, only one half of the circuit is shown. Ipsi-lateral and contra-lateral sides are shown as a left and right sides, respectively...)

Cite [Franken2014], this is an experimental description of MSO workings in line with our concepts and coincidence detection theories...
Cite [Silver2010], because this is a prospect of neuronal arithmetic...
To do: check, whether both [Smith1991] and [Smith1993] (papers from the Joris group) are relevant.
The ergodic paper of [Toth2018] is already cited, but in other context it may also go to Discussion...
Another MSO model, already studied in 2005 is: [Zhou2005], this is a great example of simplistic model, I (Petr) am much inspired by this model...
Cochlear implants and MSO, see [Kerber2012]...
(No)Cite?  Check [Wu1992]: I thought we should cite this because of high precision of some synapses, but maybe this is not necessary...
(No)Cite?  Here are the three LSO experimental papers, which have detailed methodology applicable to overlap of sound frequency ranges with the MSO; maybe none of these is needed... [Joris1995]; [Joris1996]; [Joris1998\s\do7(e)nvelope].
Cite [Jeffress1962] should go here instead of Jeffress, 1948, which is a bit too notorious...
(No)Cite?  [Bures2012] is neuronal arithmetic in Biological Cybernetics.
(No)Cite?  [Bures2013] is “dual/ LSO” paper in Brain Research
(No)Cite?  Do we still need to cite this review [Grothe2010]? ...
(No)Cite?  There are the classics: [Goldberg1969]...
Cite: Cochlear delays, very important: [Joris2006].
Yet another neuronal arithmetic: [Koutsou2012].
(No)Cite?  This is about other dimension; if we cite this, then we must discuss all the three, [Kopco2012]. Shall we do that? 
(No)Cite?  Jitter and variability by [Kostal2010]...
(No)Cite?  [Krips2009]...
(No)Cite?  [Laback2008]...
(No)Cite?  This is a chronological record of the neuronal arithmetic ideas, starting from: [Marsalek2000], through [Marsalek2001] through [Marsalek2004] to the [Marsalek2005\s\do7(K)ofranek]; watch the years ...
No(Cite)?  Some more detailed account of the diversity of cell functions and shapes: [Reed2002]...
Next possible citations: [Cariani2011], [Bouse2019], [Bures2014], [Michelet2012], [Srinivasan2018].
(Note: Other notes: Figure 4 should correct the omission in the paper by Sanda and Marsalek. When we cover everything else, we should add this: we urgently need binaural data/ psychophysics/ or electrophysiology, all the same.)
Conclusions
This theoretical paper is a followup of sound localization precision descriptions in the MSO [Sanda2012] and in the LSO [Bures2013]. Novel results are two: 1) analytical estimates of results obtained previously only by numerical simulation and 2) comparison of natural and electrical hearing modalities and demonstrating what is the maximum timing precision which can be attained by the two modalities.
Other Notes
Things to do to co-authors.
Petr Marsalek. Clean up and finish first pass writing.
Marek Hajny is joining the author team as: middle position author with tasks: 1) to compile Discussion and write and 2) to rewrite parts of the manuscript.
Zbynek Bures. Should do final reading of the manuscript before submission to the journal.
Acknowledgments
PM and MH were supported by: Charles University graduate students research program (acronym SVV) No. 260371/ 2017-2019. Special thanks to Pavel Sanda for permission to use his simulations.
Compilation notes
LATEX compilation date is: .

List of Figures

[Sorry. Ignored \begin{centering} ... \end{centering}]

Figure 1:  Vector strength of auditory nerve spike trains in dependence on sound frequency. X-axis shows sound frequency in Hz in a logarithmic scale and y-axis shows the vector strength. Even though in some nuclei up the auditory pathway the synchronization can be maintained towards higher frequencies than shown here, the decrease of the vector strength towards higher frequencies is a general property of all neurons in the auditory pathway. Red curve shows upper theoretical limit and blue curve shows lower limit. Black points are data simulated with the use of point-process-spikes train generation with the use of the dead time Poisson process.

  
[Sorry. Ignored \begin{centering} ... \end{centering}]

Figure 2:  The shortest ITD detected in the dependence on sound frequency. X-axis shows sound frequency in Hz in a logarithmic scale and y-axis the shortest JND in ms. This is a theoretical prediction based on the analytical model and basic parameter set used in simulations. 

  
[Sorry. Ignored \begin{centering} ... \end{centering}]

Figure 3:  The JND values in basic parameter set in dependence on the jitter magnitude. This plot in semi-logarithmic y-scale shows just noticeable difference of interaural time difference depending on variation of the spike timing jitter. Jagged line: simulated data, solid line: an exponential fit to the simulations under the assumption of arbitrary time precision in the model circuit, line with circles: a quadratic function estimate of spike timing precision in a system with addition of noise. 

  
[Sorry. Ignored \begin{centering} ... \end{centering}]

Figure 4:  Fit of example functions to firing rate slope ITD readout curve. X-axis shows ITD in ms and y-axis shows corresponding firing rate in action potentials per second. Note that the curve peak is offset from the origin of coordinates at . Red curve is sine function fir and blue curve is fit by the Normal density (Gauss function) with the variance set to correspond to the known sound main period. Green curve shows the difference between the red and blue curves. 

List of Tables

[Sorry. Ignored \begin{centering} ... \end{centering}]

Table 1:  The basic set of parameters. (Note: the table should be here, in analogy to the tables in [Sanda2012].)