benda/susceptibility1
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

fixes

Browse Source
This commit is contained in:
Jan Benda 2024-06-07 17:38:16 +02:00
parent 42174e8d29
commit 475053ddb2
1 changed files with 69 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

140
susceptibility1.tex
View File

@ -405,19 +405,19 @@ Neuronal processing is inherently nonlinear --- spiking thresholds or rectificat
\begin{figure*}[t]
\includegraphics[width=\columnwidth]{plot_chi2}
\caption{\label{fig:lifresponse} First- (linear) and second-order response functions of the leaky integrate-and-fire model. \figitem{A} Magnitude of the first-order response function $|\chi_1(f_1)|$, also known as ``gain'' quantifies the response amplitude relative to the stimulus amplitude, both measured at the stimulus frequency. \figitem{B} Magnitude of the second-order response function $|\chi_2(f_1, f_2)|$ quantifies the response at the sum of two stimulus frequencies. For linear systems, the second-order response functions is zero, because linear systems do not create new frequencies and thus there is no response at the summed frequencies. The plots show the analytical solutions from \cite{Lindner2001} and \cite{Voronenko2017} with $\mu = 1.1$ and $D = 0.001$.}
\caption{\label{fig:lifresponse} First- (linear) and second-order response functions of the leaky integrate-and-fire model. \figitem{A} Magnitude of the first-order response function $|\chi_1(f_1)|$, also known as ``gain'' quantifies the response amplitude relative to the stimulus amplitude, both measured at the stimulus frequency. \figitem{B} Magnitude of the second-order response function $|\chi_2(f_1, f_2)|$ quantifies the response at the sum of two stimulus frequencies. For linear systems, the second-order response functions is zero, because linear systems do not create new frequencies and thus there is no response at the summed frequencies. The plots show the analytical solutions from \citep{Lindner2001} and \citep{Voronenko2017} with $\mu = 1.1$ and $D = 0.001$.}
\end{figure*}
Nonlinear processes are key to neuronal information processing. Decision making is a fundamentally nonlinear process on the systemic level. On the cellular level, spiking neurons are inherently nonlinear: whether an action potential is elicited depends on the membrane potential to exceed a threshold\cite{Hodgkin1952,Koch1995}. Because of such nonlinearities, understanding and predicting neuronal responses to sensory stimuli is in general a difficult task.
Nonlinear processes are key to neuronal information processing. Decision making is a fundamentally nonlinear process on the systemic level. On the cellular level, spiking neurons are inherently nonlinear: whether an action potential is elicited depends on the membrane potential to exceed a threshold \citep{Hodgkin1952,Koch1995}. Because of such nonlinearities, understanding and predicting neuronal responses to sensory stimuli is in general a difficult task.
At the heart of nonlinear system identification is the Volterra series\cite{Rieke1999}. Second-order kernels have been used to predict firing rate responses of catfish retinal ganglion cells \cite{Marmarelis1972}.
In the frequency domain second-order kernels are known as second-order response functions or susceptibilities. They quantify the amplitude of the response at the sum and difference of two stimulus frequencies. Adding also third-order kernels, spike trains of a spider mechanorecptors have beend predicted from sensory stimuli \cite{French2001}. The nonlinear nature of Y cells in contrast to the more linear responses of X cells in cat retinal ganglion cells has been demonstrated by means of second-order kernels\cite{Victor1977}. Interactions between different frequencies in the response of neurons in visual cortices of cats and monkeys have been studied using bispectra, the crucial constituent of the second-order susceptibliity \cite{Schanze1997}. Locking of chinchilla auditory nerve fibres to pure tone stimuli is captured by second-order kernels\cite{Temchin2005}. In paddlefish ampullary afferents, bursting in response to strong, natural sensory stimuli boost nonlinear responses in the bicoherence, the bispectrum normalized by stimulus and response spectra \cite{Neiman2011}.
At the heart of nonlinear system identification is the Volterra series \citep{Rieke1999}. Second-order kernels have been used to predict firing rate responses of catfish retinal ganglion cells \citep{Marmarelis1972}.
In the frequency domain second-order kernels are known as second-order response functions or susceptibilities. They quantify the amplitude of the response at the sum and difference of two stimulus frequencies. Adding also third-order kernels, spike trains of a spider mechanorecptors have beend predicted from sensory stimuli \citep{French2001}. The nonlinear nature of Y cells in contrast to the more linear responses of X cells in cat retinal ganglion cells has been demonstrated by means of second-order kernels \citep{Victor1977}. Interactions between different frequencies in the response of neurons in visual cortices of cats and monkeys have been studied using bispectra, the crucial constituent of the second-order susceptibliity \citep{Schanze1997}. Locking of chinchilla auditory nerve fibres to pure tone stimuli is captured by second-order kernels \citep{Temchin2005}. In paddlefish ampullary afferents, bursting in response to strong, natural sensory stimuli boost nonlinear responses in the bicoherence, the bispectrum normalized by stimulus and response spectra \citep{Neiman2011}.
Noise linearizes nonlinear systems\cite{Yu1989,Chialvo1997} and therefore noisy neural systems can be well described by linear response theory in the limit of small stimulus amplitudes \cite{Roddey2000,Doiron2004,Rocha2007,Sharafi2013, }. When increasing stimulus amplitude, at first the contribution of the second-order kernel of the Volterra series becomes relevant in addition to the linear one. For these weakly nonlinear responses of leaky-integrate-and fire (LIF) neurons an analytical expression for the second-order susceptibility has been derived \cite{Voronenko2017} in addition to its linear response function\cite{Lindner2001}. In the superthreshold regime, where the LIF generates a baseline firing rate in the absense of an external stimulus, the linear response function has a peak at the baseline firing rate and its harmonics (\subfigrefb{fig:lifresponse}{A}) and the second-order susceptibility shows very distinct ridges of elevated nonlinear responses, exactly where two stimulus frequencies equal or add up to the neuron's baseline firing rate (\subfigrefb{fig:lifresponse}{B}). In experimental data such structures in the second-order susceptibility have not been reported yet.
Noise linearizes nonlinear systems \citep{Yu1989,Chialvo1997} and therefore noisy neural systems can be well described by linear response theory in the limit of small stimulus amplitudes \citep{Roddey2000,Doiron2004,Rocha2007,Sharafi2013, }. When increasing stimulus amplitude, at first the contribution of the second-order kernel of the Volterra series becomes relevant in addition to the linear one. For these weakly nonlinear responses of leaky-integrate-and fire (LIF) neurons an analytical expression for the second-order susceptibility has been derived \citep{Voronenko2017} in addition to its linear response function \citep{Lindner2001}. In the superthreshold regime, where the LIF generates a baseline firing rate in the absense of an external stimulus, the linear response function has a peak at the baseline firing rate and its harmonics (\subfigrefb{fig:lifresponse}{A}) and the second-order susceptibility shows very distinct ridges of elevated nonlinear responses, exactly where two stimulus frequencies equal or add up to the neuron's baseline firing rate (\subfigrefb{fig:lifresponse}{B}). In experimental data such structures in the second-order susceptibility have not been reported yet.
Here we study weakly nonlinear repsonses in two electrosensory systems in the wave-type electric fish \textit{Apteronotus leptorhynchus}. These fish generate a quasi-sinusoidal dipolar electric field (electric organ discharge, EOD). In communication contexts\cite{Walz2014, Henninger2018} the EODs of close-by fish superimpose and lead to amplitude modulations (AMs), called beats (two-fish interaction), and modulations of beats, called envelopes (multiple-fish interaction) \cite{Yu2005,Fotowat2013}. Therefore, stimuli with multiple distinct frequencies are part of the everyday life of wave-type electric fish\cite{Benda2020} and interactions of these frequencies in the electrosensory periphery are to be expected. P-type electroreceptor afferents of the tuberous electrosensory system, the P-units, use nonlinearities to extract and encode these AMs in their time-dependent firing rates \cite{Bastian1981a,Walz2014,Middleton2006,Stamper2012Envelope,Savard2011,Barayeu2023}. P-units are heterogeneous in their baseline firing rates as well as in their intrinsic noise levels, as quantified by the coefficient of variation (CV) of the interspike intervals (ISI) \cite{Grewe2017, Hladnik2023}. Low-CV P-units have a less noisy firing pattern that is closer to pacemaker firing, whereas high-CV P-units show a more irregular firing pattern that is closer to a Poisson process. On the other hand, ampullary cells of the passive electrosensory system are homogeneous in ther response properties and have very low CVs \cite{Grewe2017}.
Here we study weakly nonlinear repsonses in two electrosensory systems in the wave-type electric fish \textit{Apteronotus leptorhynchus}. These fish generate a quasi-sinusoidal dipolar electric field (electric organ discharge, EOD). In communication contexts \citep{Walz2014, Henninger2018} the EODs of close-by fish superimpose and lead to amplitude modulations (AMs), called beats (two-fish interaction), and modulations of beats, called envelopes (multiple-fish interaction) \citep{Yu2005,Fotowat2013}. Therefore, stimuli with multiple distinct frequencies are part of the everyday life of wave-type electric fish \citep{Benda2020} and interactions of these frequencies in the electrosensory periphery are to be expected. P-type electroreceptor afferents of the tuberous electrosensory system, the P-units, use nonlinearities to extract and encode these AMs in their time-dependent firing rates \citep{Bastian1981a,Walz2014,Middleton2006,Stamper2012Envelope,Savard2011,Barayeu2023}. P-units are heterogeneous in their baseline firing rates as well as in their intrinsic noise levels, as quantified by the coefficient of variation (CV) of the interspike intervals (ISI) \citep{Grewe2017, Hladnik2023}. Low-CV P-units have a less noisy firing pattern that is closer to pacemaker firing, whereas high-CV P-units show a more irregular firing pattern that is closer to a Poisson process. On the other hand, ampullary cells of the passive electrosensory system are homogeneous in ther response properties and have very low CVs \citep{Grewe2017}.
Field observations have shown that courting males were able to react to distant intruder males despite the strong EOD of the nearby female\cite{Henninger2018}. Weakly nonlinear responses are of direct relevance in this setting since they could boost responses to the faint signal of the distant intruder and thus improve detection \cite{Mascha}.
Field observations have shown that courting males were able to react to distant intruder males despite the strong EOD of the nearby female \citep{Henninger2018}. Weakly nonlinear responses are of direct relevance in this setting since they could boost responses to the faint signal of the distant intruder and thus improve detection \citep{Mascha}.
@ -439,14 +439,14 @@ Theoretical work on leaky integrate-and-fire and conductance-based models sugges
\end{figure*}
Without any external stimulation, a P-unit fires action potentials at a spontaneous baseline rate \fbase{} to the fish's own EOD of frequency \feod{}. Accordingly, a peak at \fbase{} is present in the power spectrum of this baseline activity (\subfigrefb{fig:motivation}{A}). Superposition of the receiver's EOD with an EOD of another fish with frequency $f_1$ results in a beat, a periodic amplitude modulation of the receiver's EOD. The frequency of the beat is given by the difference frequency $\Delta f_1 = f_1 - \feod$ between the two fish. P-units encode this beat in their firing rate \cite{Bastian1981a,Barayeu2023} and consequently a peak at this beat frequency appears the the power spectrum of the response (\subfigrefb{fig:motivation}{B}). A second peak at the first harmonic indicates a nonlinear response that here is easily identified by the clipping of the P-unit's firing rate at zero. Pairing the fish with another fish with a higher beat frequency $\Delta f_2 = f_2 - \feod$ results in a weaker response with a single peak in the response power spectrum (\subfigrefb{fig:motivation}{C}). Note that $\Delta f_2$ has been choosen to match the P-unit's baseline firing rate.
Without any external stimulation, a P-unit fires action potentials at a spontaneous baseline rate \fbase{} to the fish's own EOD of frequency \feod{}. Accordingly, a peak at \fbase{} is present in the power spectrum of this baseline activity (\subfigrefb{fig:motivation}{A}). Superposition of the receiver's EOD with an EOD of another fish with frequency $f_1$ results in a beat, a periodic amplitude modulation of the receiver's EOD. The frequency of the beat is given by the difference frequency $\Delta f_1 = f_1 - \feod$ between the two fish. P-units encode this beat in their firing rate \citep{Bastian1981a,Barayeu2023} and consequently a peak at this beat frequency appears the the power spectrum of the response (\subfigrefb{fig:motivation}{B}). A second peak at the first harmonic indicates a nonlinear response that here is easily identified by the clipping of the P-unit's firing rate at zero. Pairing the fish with another fish with a higher beat frequency $\Delta f_2 = f_2 - \feod$ results in a weaker response with a single peak in the response power spectrum (\subfigrefb{fig:motivation}{C}). Note that $\Delta f_2$ has been choosen to match the P-unit's baseline firing rate.
When stimulating the fish with both frequencies, additional peaks appear in the response power spectrum at the sum \bsum{} and the difference frequency \bdiff{} (\subfigrefb{fig:motivation}{D}). Thus, the response is not equal to sum of the responses to the two beats presented separately. These peaks at the sum and the difference of the two stimulus frequencies are a hallmark of nonlinear interactions that by definition are absent in linear systems.
For this example we have chosen two specific stimulus (beat) frequencies. However, for a full characterization of nonlinear repsonses we would need to measure the response of the P-units to many different combinations of stimulus frequencies. In addition, the beat stimuli in the example were strong and partially revealed saturating nonlinearities of the P-units. For weakly nonlinear responses we need to use stimuli of much lower amplitudes.
\subsection{Nonlinear signal transmission in low-CV P-units}
Weakly nonlinear responses are expected in cells with sufficiently low intrinsic noise levels, i.e. low baseline CVs \cite{Voronenko2017}. P-units fire action potentials probabilistically phase-locked to the self-generated EOD \cite{Bastian1981a}. Skipping of EOD cycles leads to the characteristic multimodal ISI distribution with maxima at integer multiples of the EOD period (\subfigrefb{fig:cells_suscept}{A}). In this example the baseline ISI distribution has a CV$_{\text{base}}$ of 0.2, which is at the lower end of the P-unit population \cite{Hladnik2023}. Spectral analysis of the baseline activity shows two major peaks: the first is located at the baseline firing rate \fbase, the second is located at the discharge frequency \feod{} of the electric organ and is flanked by two smaller peaks at $\feod \pm \fbase{}$ (\subfigref{fig:cells_suscept}{B}).
Weakly nonlinear responses are expected in cells with sufficiently low intrinsic noise levels, i.e. low baseline CVs \citep{Voronenko2017}. P-units fire action potentials probabilistically phase-locked to the self-generated EOD \citep{Bastian1981a}. Skipping of EOD cycles leads to the characteristic multimodal ISI distribution with maxima at integer multiples of the EOD period (\subfigrefb{fig:cells_suscept}{A}). In this example the baseline ISI distribution has a CV$_{\text{base}}$ of 0.2, which is at the lower end of the P-unit population \citep{Hladnik2023}. Spectral analysis of the baseline activity shows two major peaks: the first is located at the baseline firing rate \fbase, the second is located at the discharge frequency \feod{} of the electric organ and is flanked by two smaller peaks at $\feod \pm \fbase{}$ (\subfigref{fig:cells_suscept}{B}).
\begin{figure*}[t]
@ -456,9 +456,9 @@ Weakly nonlinear responses are expected in cells with sufficiently low intrinsic
Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:cells_suscept}{C}, top trace, red line), are commonly used to characterize stimulus driven responses of sensory neurons by means of transfer functions (first-order susceptibility), spike-triggered averages, or stimulus-response coherences. Here, we additionally estimate the second-order susceptibility to quantify nonlinear encoding. P-unit spikes align more or less clearly to fluctuations in the RAM stimulus. A higher stimulus intensity, here a higher contrast of the RAM relative to the EOD amplitude (see methods), entrains the P-unit response more clearly (light and dark purple for low and high contrast stimuli, \subfigrefb{fig:cells_suscept}{C}). Linear encoding, quantified by the transfer function \Eqnref{linearencoding_methods}, is similar for the two RAM contrasts in this low-CV P-unit (\subfigrefb{fig:cells_suscept}{D}), as expected for a linear system. The first-order susceptibility is low for low frequencies, peaks in the range below 100\,Hz and then falls off again \notejb{Cite Moe paper?}.
The second-order susceptibility, \Eqnref{eq:susceptibility}, quantifies amplitude and phase of the stimulus-evoked response at the sum \fsum{} for each combination of two stimulus frequencies \fone{} and \ftwo{}. Large values of the second-order susceptibility indicate stimulus-evoked peaks in the response spectrum at the summed frequency that cannot be explained by linear response theory. Similar to the first-order susceptibility, the second-order susceptibility can be estimated directly from the response evoked by a RAM stimulus that simulates the neuron with a whole range of frequencies simultaneously (\subfigsref{fig:cells_suscept}{E, F}). For a LIF driven in the super-threshold regime, theory predicts nonlinear interactions between the two stimulus frequencies, when the two frequencies \fone{} and \ftwo{} or their sum \fsum{} excatly match the neuron's baseline firing rate \fbase{} \cite{Voronenko2017}. Only then, additional stimulus-evoked peaks appear in the spectrum of the spiking response that would show up in the second-order susceptibility as a horizontal, a vertical, and an anti-diagonal line (\subfigrefb{fig:lifresponse}{B}, pink triangle in \subfigsref{fig:cells_suscept}{E, F}).
The second-order susceptibility, \Eqnref{eq:susceptibility}, quantifies amplitude and phase of the stimulus-evoked response at the sum \fsum{} for each combination of two stimulus frequencies \fone{} and \ftwo{}. Large values of the second-order susceptibility indicate stimulus-evoked peaks in the response spectrum at the summed frequency that cannot be explained by linear response theory. Similar to the first-order susceptibility, the second-order susceptibility can be estimated directly from the response evoked by a RAM stimulus that simulates the neuron with a whole range of frequencies simultaneously (\subfigsref{fig:cells_suscept}{E, F}). For a LIF driven in the super-threshold regime, theory predicts nonlinear interactions between the two stimulus frequencies, when the two frequencies \fone{} and \ftwo{} or their sum \fsum{} excatly match the neuron's baseline firing rate \fbase{} \citep{Voronenko2017}. Only then, additional stimulus-evoked peaks appear in the spectrum of the spiking response that would show up in the second-order susceptibility as a horizontal, a vertical, and an anti-diagonal line (\subfigrefb{fig:lifresponse}{B}, pink triangle in \subfigsref{fig:cells_suscept}{E, F}).
For the low-CV P-unit we observe a band of slightly elevated second-order susceptibility for the low RAM contrast at \fsumb{} (yellowish anti-diagonal between pink edges, \subfigrefb{fig:cells_suscept}{E}). This structure vanishes for the stronger stimulus (\subfigref{fig:cells_suscept}{F}). Further, the overall level of the second-order susceptibility is reduced with increasing stimulus strength. To quantify the structural changes in the susceptibility matrices we projected susceptibility values onto the diagonal by taking the mean over all anti-diagonals (\subfigrefb{fig:cells_suscept}{G}). For the low RAM contrast this projected second-order susceptibility indeed has a small peak at \fbase{} (\subfigrefb{fig:cells_suscept}{G}, dot on top line). For the higher RAM contrast, however, this peak vanishes and the overall level of the second-order susceptibility is reduced (\subfigrefb{fig:cells_suscept}{G}). The reason behind this reduction is that a RAM with a higher contrast is not only a stimulus with an increased amplitude, but also increases the total noise in the system. Increased noise is known to linearize signal transmission \cite{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017} and thus the second-order susceptibility is expected to decrease.
For the low-CV P-unit we observe a band of slightly elevated second-order susceptibility for the low RAM contrast at \fsumb{} (yellowish anti-diagonal between pink edges, \subfigrefb{fig:cells_suscept}{E}). This structure vanishes for the stronger stimulus (\subfigref{fig:cells_suscept}{F}). Further, the overall level of the second-order susceptibility is reduced with increasing stimulus strength. To quantify the structural changes in the susceptibility matrices we projected susceptibility values onto the diagonal by taking the mean over all anti-diagonals (\subfigrefb{fig:cells_suscept}{G}). For the low RAM contrast this projected second-order susceptibility indeed has a small peak at \fbase{} (\subfigrefb{fig:cells_suscept}{G}, dot on top line). For the higher RAM contrast, however, this peak vanishes and the overall level of the second-order susceptibility is reduced (\subfigrefb{fig:cells_suscept}{G}). The reason behind this reduction is that a RAM with a higher contrast is not only a stimulus with an increased amplitude, but also increases the total noise in the system. Increased noise is known to linearize signal transmission \citep{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017} and thus the second-order susceptibility is expected to decrease.
In contrast, a high-CV P-unit (CV$_{\text{base}}=0.4$) does not exhibit pronounced nonlinearities even at low stimulus contrasts (\figrefb{fig:cells_suscept_high_CV}).
@ -471,23 +471,23 @@ In contrast, a high-CV P-unit (CV$_{\text{base}}=0.4$) does not exhibit pronounc
}
\end{figure*}
Irrespective of the CV, neither of the two example P-units shows the complete expected structure of nonlinear interactions. Electric fish posses an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogeneous electric stimuli. The population of ampullary afferents is much less heterogeneous, and known for the much lower CVs of their baseline ISIs (CV$_{\text{base}}=0.06$--$0.22$)\cite{Grewe2017}. Ampullary cells do not phase-lock to the EOD and the ISIs are unimodally distributed (\subfigrefb{fig:ampullary}{A}). As a consequence of the high regularity of their baseline spiking activity, the corresponding power spectrum shows distinct peaks at the baseline firing rate \fbase{} and harmonics of it. Since the cells do not fire phase locked to the EOD, there is no peak at \feod{} (\subfigrefb{fig:ampullary}{B}). When driven by a low-contrast noise stimulus (note: this is not an AM but a stimulus that is added to the self-generated EOD, \subfigref{fig:ampullary}{C}), ampullary cells exhibit very pronounced bands in the second-order susceptibility, where \fsum{} is equal to \fbase{} or its harmonic (yellow diagonals in \subfigrefb{fig:ampullary}{E}), implying strong nonlinear response components at these frequency combinations (\subfigrefb{fig:ampullary}{G}, top). With higher stimulus contrasts these bands disappear (\subfigrefb{fig:ampullary}{F}), the projection onto the diagonal looses its distinct peak at \fsum{} and its overal level is reduced (\subfigrefb{fig:ampullary}{G}, bottom).
Irrespective of the CV, neither of the two example P-units shows the complete expected structure of nonlinear interactions. Electric fish posses an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogeneous electric stimuli. The population of ampullary afferents is much less heterogeneous, and known for the much lower CVs of their baseline ISIs (CV$_{\text{base}}=0.06$--$0.22$) \citep{Grewe2017}. Ampullary cells do not phase-lock to the EOD and the ISIs are unimodally distributed (\subfigrefb{fig:ampullary}{A}). As a consequence of the high regularity of their baseline spiking activity, the corresponding power spectrum shows distinct peaks at the baseline firing rate \fbase{} and harmonics of it. Since the cells do not fire phase locked to the EOD, there is no peak at \feod{} (\subfigrefb{fig:ampullary}{B}). When driven by a low-contrast noise stimulus (note: this is not an AM but a stimulus that is added to the self-generated EOD, \subfigref{fig:ampullary}{C}), ampullary cells exhibit very pronounced bands in the second-order susceptibility, where \fsum{} is equal to \fbase{} or its harmonic (yellow diagonals in \subfigrefb{fig:ampullary}{E}), implying strong nonlinear response components at these frequency combinations (\subfigrefb{fig:ampullary}{G}, top). With higher stimulus contrasts these bands disappear (\subfigrefb{fig:ampullary}{F}), the projection onto the diagonal looses its distinct peak at \fsum{} and its overal level is reduced (\subfigrefb{fig:ampullary}{G}, bottom).
\subsection{Model-based estimation of the nonlinear structure}
In the example recordings shown above (\figsrefb{fig:cells_suscept} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagnal of the second-order susceptibility, where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, which is in line with theoretical expectations\cite{Voronenko2017}. However, a pronounced nonlinear response at frequencies \foneb{} or \ftwob{}, although predicted by theory, cannot be observed. Here we investigate how these discrepancies can be understood.
In the example recordings shown above (\figsrefb{fig:cells_suscept} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagnal of the second-order susceptibility, where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, which is in line with theoretical expectations \citep{Voronenko2017}. However, a pronounced nonlinear response at frequencies \foneb{} or \ftwob{}, although predicted by theory, cannot be observed. Here we investigate how these discrepancies can be understood.
\begin{figure*}[t]
\includegraphics[width=\columnwidth]{model_and_data}
\caption{\label{model_and_data} Estimation of second-order susceptibilities in the limit of weak stimuli. \figitem{A} \suscept{} estimated from $N=11$ trials \notejb{of duration XXX?} of an electrophysiological recording of another low-CV P-unit (cell 2012-07-03-ak, $\fbase=120$\,Hz, CV=0.20) driven with a weak RAM stimulus with contrast 2.5\,\%. Pink edges mark baseline firing rate where enhanced nonlinear responses are expected. \figitem[i]{B} \textit{Standard condition} of model simulations with intrinsic noise (bottom) and a RAM stimulus (top). \notejb{Since the model overestimated the sensitivity of the real P-unit, we adjusted the RAM contrast to 0.009\,\%, such that the resulting spike trains had the same CV as the electrophysiolgical recorded P-unit during the 2.5\,\% contrast stimulation (see table~\ref{modelparams} for model parameters).}\notejb{Specify fbase and baseline CV} \figitem[ii]{B} \suscept{} estimated from simulations of the cell's LIF model counterpart (cell 2012-07-03-ak, table~\ref{modelparams}) based on the same number of trials as in the electrophysiological recording. \figitem[iii]{B} Same as \panel[ii]{B} but using $10^6$ stimulus repetitions. \figitem[i-iii]{C} Same as in \panel[i-iii]{B} but in the \textit{noise split} condition: there is no external RAM signal driving the model. Instead, a large part (90\,\%) of the total intrinsic noise is treated as signal and is presented as an equivalent amplitude modulation (\signalnoise, center), while the intrinsic noise is reduced to 10\,\% of its original strength (see methods for details). In addition to one million trials, this reveals the full expected structure of the second-order susceptibility.}
\caption{\label{model_and_data} Estimation of second-order susceptibilities in the limit of weak stimuli. \figitem{A} \suscept{} estimated from $N=11$ trials \notejb{of duration XXX?} of an electrophysiological recording of another low-CV P-unit (cell 2012-07-03-ak, $\fbase=120$\,Hz, CV=0.20) driven with a weak RAM stimulus with contrast 2.5\,\%. Pink edges mark baseline firing rate where enhanced nonlinear responses are expected. \figitem[i]{B} \textit{Standard condition} of model simulations with intrinsic noise (bottom) and a RAM stimulus (top). \notejb{Since the model overestimated the sensitivity of the real P-unit, we adjusted the RAM contrast to 0.009\,\%, such that the resulting spike trains had the same CV as the electrophysiolgical recorded P-unit during the 2.5\,\% contrast stimulation (see table~\ref{modelparams} for model parameters).} \figitem[ii]{B} \suscept{} estimated from simulations of the cell's LIF model counterpart (cell 2012-07-03-ak, table~\ref{modelparams}) based on the same number of trials as in the electrophysiological recording. \figitem[iii]{B} Same as \panel[ii]{B} but using $10^6$ stimulus repetitions. \figitem[i-iii]{C} Same as in \panel[i-iii]{B} but in the \textit{noise split} condition: there is no external RAM signal driving the model. Instead, a large part (90\,\%) of the total intrinsic noise is treated as signal and is presented as an equivalent amplitude modulation (\signalnoise, center), while the intrinsic noise is reduced to 10\,\% of its original strength (see methods for details). In addition to one million trials, this reveals the full expected structure of the second-order susceptibility.}
\end{figure*}
One reason could be simply too little data for a good estimate of the second-order susceptibility. Electrophysiological recordings are limited in time, and therefore responses to only a limited number of trials, i.e. repetitions of the same RAM stimulus, are available. As a consequence, the cross-spectra, \Eqnref{eq:crosshigh}, are insufficiently averaged and the full structure of the second-order susceptibility might be hidden in finite-data noise. This experimental limitation can be overcome by using a computational model for the P-unit, a stochastic leaky integrate-and-fire model with adaptation current and dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \cite{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of a low-CV cell estimated from the same low number of trials as in the experiment ($\n{}=11$, compare \panel{A} and \panel[ii]{B} in \figrefb{model_and_data}).
One reason could be simply too little data for a good estimate of the second-order susceptibility. Electrophysiological recordings are limited in time, and therefore responses to only a limited number of trials, i.e. repetitions of the same RAM stimulus, are available. As a consequence, the cross-spectra, \Eqnref{eq:crosshigh}, are insufficiently averaged and the full structure of the second-order susceptibility might be hidden in finite-data noise. This experimental limitation can be overcome by using a computational model for the P-unit, a stochastic leaky integrate-and-fire model with adaptation current and dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of a low-CV cell estimated from the same low number of trials as in the experiment ($\n{}=11$, compare \panel{A} and \panel[ii]{B} in \figrefb{model_and_data}).
In simulations of the model we can increase the number of trials beyond what would be experimentally possible, here to one million (\subfigrefb{model_and_data}\,\panel[iii]{B}). The estimate of the second-order susceptibility indeed improves. It gets less noisy and the diagonal at \fsum{} is emphasized. However, the expected vertical and horizontal lines at \foneb{} and \ftwob{} are still missing.
Using a broadband stimulus increases the effective input-noise level and this may linearize signal transmission and suppress potential nonlinear responses \cite{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017}. Assuming that the intrinsic noise level in this P-unit is small enough, the full expected structure of the second-order susceptibility should appear in the limit of weak AMs. Again, this cannot be done experimentally, because the problem of insufficient averaging becomes even more severe for weak AMs (low contrast). In the model, however, we know the time course of the intrinsic noise and can use this knowledge to determine the susceptibilities by input-output correlations via the Furutsu-Novikov theorem \cite{Furutsu1963, Novikov1965}. This theorem, in its simplest form, states that the cross-spectrum $S_{x\eta}(\omega)$ of a Gaussian noise $\eta(t)$ driving a nonlinear system and the system's output $x(t)$ is proportional to the linear susceptibility according to $S_{x\eta}(\omega)=\chi(\omega)S_{\eta\eta}(\omega)$. Here $\chi(\omega)$ characterizes the linear response to an infinitely weak signal $s(t)$ in the presence of the background noise $\eta(t)$. Likewise, the nonlinear susceptibility can be determined in an analogous fashion from higher-order input-output cross-spectra (see methods, equations \eqref{eq:crosshigh} and \eqref{eq:susceptibility}) \cite{Egerland2020}. In line with an alternative derivation of the Furutsu-Novikov theorem \cite{Lindner2022}, we can split the total noise and consider a fraction of it as stimulus. This allows to calculate the susceptibility from the cross-spectrum between the output and this stimulus fraction of the noise. Adapting this approach to our P-unit model (see methods), we replace the intrinsic noise by an approximately equivalent RAM stimulus $s_\xi(t)$ and a weak remaining intrinsic noise $\sqrt{2D \, c_{\rm{noise}}} \cdot \xi(t)$ with $c_\text{noise} = 0.1$ (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \subfigrefb{model_and_data}\,\panel[i]{C}). We tune the amplitude of the RAM stimulus $s_\xi(t)$ such that the output firing rate and variability (CV) are the same as in the baseline activity (i.e. full intrinsic noise $\sqrt{2D}\xi(t)$ in the voltage equation but no RAM) and compute the cross-spectra between the RAM part of the noise $s_\xi(t)$ and the output spike train. This procedure has two consequences: (i) by means of the cross-spectrum between the output and \signalnoise, which is a large fraction of the noise, the signal-to-noise ratio of the measured susceptibilities is drastically improved; (ii) the total noise in the system has been reduced (by what was before the external RAM stimulus $s(t)$), which makes the system more nonlinear. For both reasons we now see the expected nonlinear features in the second-order susceptibility for a sufficient number of trials (\subfigrefb{model_and_data}\,\panel[iii]{C}), but not for a number of trials comparable to the experiment (\subfigrefb{model_and_data}\,\panel[ii]{C}). In addition to the strong response for \fsumb{}, we now also observe pronounced nonlinear responses at \foneb{} and \ftwob{} (vertical and horizontal lines, \subfigrefb{model_and_data}\,\panel[iii]{C}).
Using a broadband stimulus increases the effective input-noise level and this may linearize signal transmission and suppress potential nonlinear responses \citep{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017}. Assuming that the intrinsic noise level in this P-unit is small enough, the full expected structure of the second-order susceptibility should appear in the limit of weak AMs. Again, this cannot be done experimentally, because the problem of insufficient averaging becomes even more severe for weak AMs (low contrast). In the model, however, we know the time course of the intrinsic noise and can use this knowledge to determine the susceptibilities by input-output correlations via the Furutsu-Novikov theorem \citep{Furutsu1963, Novikov1965}. This theorem, in its simplest form, states that the cross-spectrum $S_{x\eta}(\omega)$ of a Gaussian noise $\eta(t)$ driving a nonlinear system and the system's output $x(t)$ is proportional to the linear susceptibility according to $S_{x\eta}(\omega)=\chi(\omega)S_{\eta\eta}(\omega)$. Here $\chi(\omega)$ characterizes the linear response to an infinitely weak signal $s(t)$ in the presence of the background noise $\eta(t)$. Likewise, the nonlinear susceptibility can be determined in an analogous fashion from higher-order input-output cross-spectra (see methods, equations \eqref{eq:crosshigh} and \eqref{eq:susceptibility}) \citep{Egerland2020}. In line with an alternative derivation of the Furutsu-Novikov theorem \citep{Lindner2022}, we can split the total noise and consider a fraction of it as stimulus. This allows to calculate the susceptibility from the cross-spectrum between the output and this stimulus fraction of the noise. Adapting this approach to our P-unit model (see methods), we replace the intrinsic noise by an approximately equivalent RAM stimulus $s_\xi(t)$ and a weak remaining intrinsic noise $\sqrt{2D \, c_{\rm{noise}}} \cdot \xi(t)$ with $c_\text{noise} = 0.1$ (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \subfigrefb{model_and_data}\,\panel[i]{C}). We tune the amplitude of the RAM stimulus $s_\xi(t)$ such that the output firing rate and variability (CV) are the same as in the baseline activity (i.e. full intrinsic noise $\sqrt{2D}\xi(t)$ in the voltage equation but no RAM) and compute the cross-spectra between the RAM part of the noise $s_\xi(t)$ and the output spike train. This procedure has two consequences: (i) by means of the cross-spectrum between the output and \signalnoise, which is a large fraction of the noise, the signal-to-noise ratio of the measured susceptibilities is drastically improved; (ii) the total noise in the system has been reduced (by what was before the external RAM stimulus $s(t)$), which makes the system more nonlinear. For both reasons we now see the expected nonlinear features in the second-order susceptibility for a sufficient number of trials (\subfigrefb{model_and_data}\,\panel[iii]{C}), but not for a number of trials comparable to the experiment (\subfigrefb{model_and_data}\,\panel[ii]{C}). In addition to the strong response for \fsumb{}, we now also observe pronounced nonlinear responses at \foneb{} and \ftwob{} (vertical and horizontal lines, \subfigrefb{model_and_data}\,\panel[iii]{C}).
Note, that the increased number of trials goes along with a substantial reduction of second-order susceptibility values (\subfigrefb{model_and_data}\,\panel[iii]{C}), that saturate in its peak values for $N>10^5$ (\figrefb{fig:trialnr}). This demonstrates the limited reliability of an estimate of the second-order susceptibility that is based on 11 trials only. However, we would like to point out that already the limited number of trials as used in the experiments reveals key features of the nonlinear response.
@ -502,10 +502,11 @@ We estimated the second-order susceptibility of P-unit responses using RAM stimu
\notejb{zero missing on y-axis in B}
\notejb{add tick at 100Hz in E \& F}
\notejb{C-F tags a bit lower, D\&F tag closer to yaxis}
\notejb{What do we want in B? square or linear or dB, noise split as in A, or stimulus as in C-F }
\caption{\label{fig:model_full} Using second-order susceptibility to predict responses to sine-wave stimuli. \figitem[]{A} Absolute value of the second-order susceptibility, \Eqnref{eq:susceptibility}, for both positive and negative frequencies. \susceptf{} was estimated from $N=10^6$ trials of model simulations in the noise-split condition (cell 2012-07-03-ak, see table~\ref{modelparams} for model parameters). White lines indicate zero frequencies. Nonlinear responses at \fsum{} are quantified in the upper right and lower left quadrants. Nonlinear responses at \fdiff{} are quantified in the upper left and lower right quadrants. Baseline firing rate of this cell was at $\fbase=120$\,Hz. The position of the orange/red letters correspond to the beat frequencies used for the stimulation with pure sine-waves in the subsequent panels and indicate the sum/difference of those beat frequencies. \figitem[]{B} Absolute value of the first-order susceptibility. \figitem{C--F} Black line -- power spectral density of model simulations in response to stimulation with two pure sine waves, \fone{} and \ftwo, in addition to the receiving fish's own EOD (three fish scenario). The contrasts of beat beats is 0.0065. Colored circles highlight the height of selected peaks in the power spectrum. Black circles highlight the peak height that can be predicted from \panel{A, B}. Grey line -- power spectral density of model in the baseline condition. \figitem{B} The sum of the two beat frequencies match \fbase{}. \figitem{C} The difference of \fone{} and \ftwo{} match \fbase{}. \figitem{D} Only the first beat frequency matches \fbase{}. \figitem{C} None of the two beat frequencies matches \fbase{}. \noteab{Für die Transfer Funktion habe ich jetzt einen Faktor 1, für die Nichtlinearität einen Faktor 30, aber vielleicht wenn ich über mehrere Punkte mitteln muss und das alles so noisy ist das eben noch keine Gute Abschätzung in der Stauration?}}
\end{figure*}
However, the second-order susceptibility \Eqnref{eq:susceptibility} is a spectral measure that is based on Fourier transforms and thus is also defined for negative stimulus frequencies. The full \susceptf{} matrix is symmetric with respect to the origin. In the upper-right and lower-left quadrants of \susceptf{}, stimulus-evoked responses at \fsum{} are shown, whereas in the lower-right and upper-left quadrants nonlinear responses at the difference \fdiff{} are shown (\figref{fig:model_full}). The vertical and horizontal lines at \foneb{} and \ftwob{} are very pronounced in the upper-right quadrant of \subfigrefb{fig:model_full}{A} for nonlinear responses at \fsum{} and extend into the upper-left quadrant (representing \fdiff) where they fade out towards more negative $f_1$ frequencies. The peak in the response power-spectrum at \fdiff{} evoked by pure sine-wave stimulation (\subfigrefb{fig:motivation}{D}) is predicted by the horizontal line in the upper-left quadrant (\subfigrefb{fig:model_full}{A}, \cite{Schlungbaum2023}).
However, the second-order susceptibility \Eqnref{eq:susceptibility} is a spectral measure that is based on Fourier transforms and thus is also defined for negative stimulus frequencies. The full \susceptf{} matrix is symmetric with respect to the origin. In the upper-right and lower-left quadrants of \susceptf{}, stimulus-evoked responses at \fsum{} are shown, whereas in the lower-right and upper-left quadrants nonlinear responses at the difference \fdiff{} are shown (\figref{fig:model_full}). The vertical and horizontal lines at \foneb{} and \ftwob{} are very pronounced in the upper-right quadrant of \subfigrefb{fig:model_full}{A} for nonlinear responses at \fsum{} and extend into the upper-left quadrant (representing \fdiff) where they fade out towards more negative $f_1$ frequencies. The peak in the response power-spectrum at \fdiff{} evoked by pure sine-wave stimulation (\subfigrefb{fig:motivation}{D}) is predicted by the horizontal line in the upper-left quadrant (\subfigrefb{fig:model_full}{A}, \citep{Schlungbaum2023}).
Is it possible based on the second-order susceptibility estimated by means of RAM stimuli (\subfigrefb{fig:model_full}{A}) to predict nonlinear responses in a three-fish setting? We can test this by stimulating the same model with two beats with weak amplitudes (\subfigrefb{fig:model_full}{C--F}). If we choose a frequency combination where the sum of the two beat frequencies is equal to the model's baseline firing rate \fbase{}, a peak at the sum of the two beat frequencies appears in the power spectrum of the response (\subfigrefb{fig:model_full}{C}), as expected from \suscept. If instead we choose two beat frequencies that differ by \fbase{}, a peak is present at the difference frequency (\subfigrefb{fig:model_full}{D}). If only one beat frequency is equal to \fbase{}, both a peak at the sum and at the difference frequency is present in the P-unit response (\subfigrefb{fig:model_full}{E}). And if none of these conditions are met, neither a peak at the sum nor at the difference of the two beat frequencies appears (\subfigrefb{fig:model_full}{F}).
\begin{figure*}[tp]
@ -521,7 +522,7 @@ Is it possible based on the second-order susceptibility estimated by means of RA
\subsection{Low CVs and weak stimuli are associated with strong nonlinearity}
All the statements about nonlinear encoding in p-type and ampullary electroreceptor afferents based on single cell examples shown above are supported by the analysis of our pool of 221 P-units and 47 ampullary afferents recorded in 71 specimen. For a comparison across cells we summarize the structure of the second-order susceptibilities in a single number, the peakedness of the nonlinearity \nli{} \Eqnref{eq:nli_equation}, that characterizes the size of the expected peak of the projections of a \suscept{} matrix onto its diagonal at the baseline firing rate (e.g. \subfigref{fig:cells_suscept}{G}). \nli{} assumes high values when the peak at \fbase{} is pronounced relative to the median of projections onto the diagonal and is small when there is no distinct peak. The \nli{} values of the P-unit population depend weakly on the CV of the baseline ISI distribution. Cells with lower baseline CVs tend to have more pronounced peaks in their projections than those that have high baseline CVs (\subfigrefb{fig:data_overview}{A}). This negative correlation is more pronounced against the CV measured during stimulation (\subfigrefb{fig:data_overview}{C}).
The effective stimulus strength also plays an important role. We quantify the effect a stimulus has on a cell's response by the response modulation --- the standard deviation of a cell's firing rate in response to a RAM stimulus. P-units are heterogeneous in their sensitivity, their response modulations to the same stimulus contrast vary a lot \cite{Grewe2017}. Cells with weak responses to a stimulus, be it an insensitive cell or a weak stimulus, have higher \nli{} values and thus a more pronounced ridge in the second-order susceptibility at \fsumb{} in comparison to strongly responding cells that basically have flat second-order susceptibilities (\subfigrefb{fig:data_overview}{E}). How pronounced nonlinear response components are in P-units thus depends on the baseline CV (a proxy for the internal noise level), and both the CV and response strength during stimulation (effective output noise).
The effective stimulus strength also plays an important role. We quantify the effect a stimulus has on a cell's response by the response modulation --- the standard deviation of a cell's firing rate in response to a RAM stimulus. P-units are heterogeneous in their sensitivity, their response modulations to the same stimulus contrast vary a lot \citep{Grewe2017}. Cells with weak responses to a stimulus, be it an insensitive cell or a weak stimulus, have higher \nli{} values and thus a more pronounced ridge in the second-order susceptibility at \fsumb{} in comparison to strongly responding cells that basically have flat second-order susceptibilities (\subfigrefb{fig:data_overview}{E}). How pronounced nonlinear response components are in P-units thus depends on the baseline CV (a proxy for the internal noise level), and both the CV and response strength during stimulation (effective output noise).
%(Pearson's $r=-0.35$, $p<0.001$)221 P-units and 47 (Pearson's $r=-0.16$, $p<0.01$)
%In a P-unit population where each cell is represented not by several contrasts but by the lowest recorded contrast, \nli{} significantly correlates with the CV during baseline ($r=-0.17$, $p=0.01$), the response modulation ($r=-0.35$, $p<0.001$) and \fbase{} ($r=-0.32$, $p<0.001$).%, $\n{}=221$*, $\n{}=221$******, $\n{}=221$
@ -531,42 +532,39 @@ The population of ampullary cells is generally more homogeneous, with lower base
\section{Discussion}
\notejb{Leftovers from discussion}
Nonlinearities are ubiquitous in nervous systems, they are essential to extract certain features such as amplitude modulations of a carrier. Here, we analyzed the nonlinearity in primary electroreceptor afferents of the weakly electric fish \lepto{}. Under natural stimulus conditions, when the electric fields of two or more animals produce interference patterns that contain information of the context of the encounter, such nonlinear encoding is required to enable responding to AMs or also so-called envelopes. The work presented here is inspired by observations of electric fish interactions in the natural habitat. In these observed electrosensory cocktail-parties, a second male is intruding the conversations of a courting dyad. The courting male detects the intruder at a large distance where the foreign electric signal is strongly attenuated due to spatial distance and tiny compared to the signal emitted by the close-by female. Such a three-animal situation is exemplified in \figref{fig:motivation} for a rather specific combination of interfering EODs. There, we observe that the electroreceptor response contains components that result from nonlinear interference (\subfigref{fig:motivation}{D}) that might help to solve this detection task. Based on theoretical work we work out the circumstances under which electroreceptor afferents show such nonlinearities.
\notejb{Even though the extraction of the AM itself requires a \notejb{nonlinearity} \citep{Middleton2006,Stamper2012Envelope,Savard2011,Barayeu2023} encoding the time-course of the AM is linear over a wide range of AM amplitudes and frequencies \citep{Xu1996,Benda2005,Gussin2007,Grewe2017,Savard2011}. In the context of social signalling among three fish we observe an AM of the AM, also referred to as second-order envelope or just social envelope \citep{Middleton2006, Savard2011, Stamper2012Envelope}. Encoding this again requires nonlinearities \citep{Middleton2006} and it was shown that a subpopulation of P-units are sensitive to envelopes \citep{Savard2011} and exhibit nonlinearities e.g. when driven by strong stimuli \citep{Nelson1997,Chacron2004}.}
\notejb{Leftovers from introduction}
Phase-locking to the own field also leads to a representation of \feod{} in the P-unit firing rate (see \figref{fig:cells_suscept}\panel{B}) \citep{Sinz2020}.
\notejb{Even though the extraction of the AM itself requires a \notejb{nonlinearity}\cite{Middleton2006,Stamper2012Envelope,Savard2011,Barayeu2023} encoding the time-course of the AM is linear over a wide range of AM amplitudes and frequencies\cite{Xu1996,Benda2005,Gussin2007,Grewe2017,Savard2011}. In the context of social signalling among three fish we observe an AM of the AM, also referred to as second-order envelope or just social envelope\cite{Middleton2006, Savard2011, Stamper2012Envelope}. Encoding this again requires nonlinearities\cite{Middleton2006} and it was shown that a subpopulation of P-units are sensitive to envelopes\cite{Savard2011} and exhibit nonlinearities e.g. when driven by strong stimuli\cite{Nelson1997,Chacron2004}.}
Phase-locking to the own field also leads to a representation of \feod{} in the P-unit firing rate (see \figref{fig:cells_suscept}\panel{B})\cite{Sinz2020}.
The appearing difference peak is known as the social envelope\cite{Stamper2012Envelope, Savard2011}. The neuron shown here, clearly encodes the envelope. Whether P-units in general encode envelopes has been subject of controversy, some works do not consider P-units as envelope encoders\cite{Middleton2006}, while others identify some P-units as envelope encoders\cite{Savard2011}.
The appearing difference peak is known as the social envelope \citep{Stamper2012Envelope, Savard2011}. The neuron shown here, clearly encodes the envelope. Whether P-units in general encode envelopes has been subject of controversy, some works do not consider P-units as envelope encoders \citep{Middleton2006}, while others identify some P-units as envelope encoders \citep{Savard2011}.
\notejb{Weakly nonlinear responses versus saturation regime}
\notejb{Estimating the infinite Volterra series from limited experimental data is usually limited to the first two or three kernels, that then might not be sufficient for a proper prediction of the neuronal response \cite{French2001}. Making assumptions about the nonlinearities in a system reduces the amount of data needed for parameter estimation. In particular, models combining linear filtering with static nonlinearities\cite{Chichilnisky2001}, have been successful in capturing functionally relevant neuronal computations in visual \cite{Gollisch2009} as well as auditory systems\cite{Clemens2013}. On the other hand, linear methods based on backward models for estimating the stimulus from neuronal responses, have been extensively used to quantify information transmission in neural systems \cite{Theunissen1996,Borst1999,Wessel1996,Chacron?}, because backward models do not need to generate action potentials \cite{Rieke1999}.}
\notejb{Estimating the infinite Volterra series from limited experimental data is usually limited to the first two or three kernels, that then might not be sufficient for a proper prediction of the neuronal response \citep{French2001}. Making assumptions about the nonlinearities in a system reduces the amount of data needed for parameter estimation. In particular, models combining linear filtering with static nonlinearities \citep{Chichilnisky2001}, have been successful in capturing functionally relevant neuronal computations in visual \citep{Gollisch2009} as well as auditory systems \citep{Clemens2013}. On the other hand, linear methods based on backward models for estimating the stimulus from neuronal responses, have been extensively used to quantify information transmission in neural systems \citep{Theunissen1996,Borst1999,Wessel1996,Chacron?}, because backward models do not need to generate action potentials \citep{Rieke1999}.}
\notejb{Note in the discussion that \suscept{} predicts stimulus evoked responses whereas in the power spectrum we see peaks even if they are unrelated to the stimulus. Like the baseline firing rate.}
% TO DISCUSSION:
%Even though the second-order susceptibilities here were estimated from data and models with an modulated (EOD) carrier (\figrefb{fig:model_full}) they are in good accordance with the second-order susceptibilities found in LIF models without a carrier\cite{Voronenko2017, Schlungbaum2023}.
%Even though the second-order susceptibilities here were estimated from data and models with an modulated (EOD) carrier (\figrefb{fig:model_full}) they are in good accordance with the second-order susceptibilities found in LIF models without a carrier \citep{Voronenko2017, Schlungbaum2023}.
%\,\panel[iii]{C}
\subsection{Theory applies to systems with and without carrier}
Theoretical work\cite{Voronenko2017} explained analytically the occurrence of nonlinear products when a LIF model neuron is stimulated with pure sine-waves. To investigate whether the same mechanisms occur in electroreceptor afferents which are driven by AMs of a carrier and not by pure sine-waves, we followed the previous approach and quantified the second-order susceptibility from responses to white-noise stimuli \cite{Voronenko2017, Egerland2020, Neiman2011fish,Nikias1993}. We expected to see elevated second-order susceptibility where either of the foreign signals matches the baseline firing rate ($f_1=f_{base}$ or $f_2=f_{base}$) or when the sum equals the baseline firing rate of the neuron (\fsumb{}) creating a triangular pattern of elevated \suscept{} e.g.\,\subfigref{model_and_data}\,\panel[iii]{C}. Indeed, we find traces of the same nonlinearities in the neuronal responses of p-type electroreceptor afferents. The nonlinear pattern observed in the experimental data, however, matches to the expectations only partially and only in a subset of neurons (\figsref{fig:cells_suscept} and\,\ref{fig:ampullary}). Nevertheless, the theory holds also for systems that are driven by AMs of a carrier and is thus more widely applicable.
Theoretical work \citep{Voronenko2017} explained analytically the occurrence of nonlinear products when a LIF model neuron is stimulated with pure sine-waves. To investigate whether the same mechanisms occur in electroreceptor afferents which are driven by AMs of a carrier and not by pure sine-waves, we followed the previous approach and quantified the second-order susceptibility from responses to white-noise stimuli \citep{Voronenko2017, Egerland2020, Neiman2011fish,Nikias1993}. We expected to see elevated second-order susceptibility where either of the foreign signals matches the baseline firing rate ($f_1=f_{base}$ or $f_2=f_{base}$) or when the sum equals the baseline firing rate of the neuron (\fsumb{}) creating a triangular pattern of elevated \suscept{} e.g.\,\subfigref{model_and_data}\,\panel[iii]{C}. Indeed, we find traces of the same nonlinearities in the neuronal responses of p-type electroreceptor afferents. The nonlinear pattern observed in the experimental data, however, matches to the expectations only partially and only in a subset of neurons (\figsref{fig:cells_suscept} and\,\ref{fig:ampullary}). Nevertheless, the theory holds also for systems that are driven by AMs of a carrier and is thus more widely applicable.
\subsection{Intrinsic noise limits nonlinear responses}
Only those P-units that exhibit low coefficients of variation (CV) of the interspike-interval distribution (\subfigref{fig:cells_suscept}{A}) in their unperturbed baseline response show the expected nonlinerities (\subfigref{fig:data_overview}{A}). Such low-CV cells are rare among the 221 P-units that we used in this study. The afferents of the passive electrosensory system, the ampullary cells, however have generally lower CVs and show a much clearer nonlinearity pattern than the low-CV P-unit exemplified here (compare \figsref{fig:cells_suscept} and \ref{fig:ampullary}). The single ampullary cell featured in \figref{fig:ampullary} is a representative of the majority of ampullary cells analyzed here. All ampullary cells have CVs below 0.4 with a median around 0.12 and the observed \nli{}s are 10-fold higher than in P-units.
The CV serves as a proxy for the intrinsic noise in the cells. In both cell types, we observe a negative correlation between \nli{} and the CV, indicating that it is the level of intrinsic noise that plays a role here. These findings are in line with previous studies that propose that noise linearizes the system\cite{Roddey2000, Chialvo1997, Voronenko2017}. More intrinsic noise has been demonstrated to increase the CV and reduce nonlinear phase-locking in vestibular afferents\cite{Schneider2011}. Reduced noise, on the other hand, has been associated with stronger nonlinearity in pyramidal cells of the ELL\cite{Chacron2006}. Further support for the notion of noise limiting the nonlinearity comes from our P-unit LIF model that faithfully reproduces P-unit activity\cite{Barayeu2023}. We can use this model and apply a noise-split \cite{Lindner2022} based on the to the Furutsu-Novikov theorem\cite{Novikov1965,Furutsu1963}, to increase the signal-to-noise ratio in the cell while keeping the overall response variability constant (see methods). Treating 90\,\% of the total noise as signal and simulating large numbers of trial uncovers the full nonlinear structure (\figref{model_and_data}) seen in LIF neurons and the analytical derivations when driven with sine-wave stimuli\cite{Voronenko2017}.
The CV serves as a proxy for the intrinsic noise in the cells. In both cell types, we observe a negative correlation between \nli{} and the CV, indicating that it is the level of intrinsic noise that plays a role here. These findings are in line with previous studies that propose that noise linearizes the system \citep{Roddey2000, Chialvo1997, Voronenko2017}. More intrinsic noise has been demonstrated to increase the CV and reduce nonlinear phase-locking in vestibular afferents \citep{Schneider2011}. Reduced noise, on the other hand, has been associated with stronger nonlinearity in pyramidal cells of the ELL \citep{Chacron2006}. Further support for the notion of noise limiting the nonlinearity comes from our P-unit LIF model that faithfully reproduces P-unit activity \citep{Barayeu2023}. We can use this model and apply a noise-split \citep{Lindner2022} based on the to the Furutsu-Novikov theorem \citep{Novikov1965,Furutsu1963}, to increase the signal-to-noise ratio in the cell while keeping the overall response variability constant (see methods). Treating 90\,\% of the total noise as signal and simulating large numbers of trial uncovers the full nonlinear structure (\figref{model_and_data}) seen in LIF neurons and the analytical derivations when driven with sine-wave stimuli \citep{Voronenko2017}.
%
\subsection{Noise stimulation approximates the real three-fish interaction}
Our analysis is based on the neuronal responses to white noise stimulus sequences. For the P-units, the stimulus was a random amplitude modulation (RAM) while it was a direct noise stimulus for the ampullary cells. These broad-band stimuli have the advantage that all behaviorally relevant frequencies can be measured with a single stimulus presentation and it is a widely used approach to characterize sensory coding\cite{French1973,Marmarelis1999, Borst1999, Chacron2005, Grewe2017}. However, these stimuli also increase the total level of noise in the system and may have a linearizing effect on signal transmission. In a our P-unit models, we were able to make use of the Furutsu-Novikov theorem to estimate nonlinear signal transmission for zero-amplitude stimulation. Only with this procedure and many trials, we find for low-CV P-units all the pronounced ridges in the second-order susceptibility that we would expect according to theory \cite{Voronenko2017}.
Our analysis is based on the neuronal responses to white noise stimulus sequences. For the P-units, the stimulus was a random amplitude modulation (RAM) while it was a direct noise stimulus for the ampullary cells. These broad-band stimuli have the advantage that all behaviorally relevant frequencies can be measured with a single stimulus presentation and it is a widely used approach to characterize sensory coding \citep{French1973,Marmarelis1999, Borst1999, Chacron2005, Grewe2017}. However, these stimuli also increase the total level of noise in the system and may have a linearizing effect on signal transmission. In a our P-unit models, we were able to make use of the Furutsu-Novikov theorem to estimate nonlinear signal transmission for zero-amplitude stimulation. Only with this procedure and many trials, we find for low-CV P-units all the pronounced ridges in the second-order susceptibility that we would expect according to theory \citep{Voronenko2017}.
In the natural situation, the stimuli are periodic signals defined by the difference frequencies. How well can we extrapolate from the white noise analysis to the pure sinewave situation? \notejg{Predictions from the X2 matrix and the equations in Voronekov}
@ -574,28 +572,28 @@ In contrast to the situation with individual frequencies (direct sine-waves or s
% The nonlinearity of ampullary cells in paddlefish \cite{Neiman2011fish} has been previously accessed with bandpass limited white noise.
% The nonlinearity of ampullary cells in paddlefish \citep{Neiman2011fish} has been previously accessed with bandpass limited white noise.
% Here it was demonstrated that the second-order susceptibility for the two RAM noise input frequencies \fone{} and \ftwo{} can approximate a three-fish setting, where the driving force for the P-unit are two beats with frequencies \bone{} and \btwo{}. This was confirmed by a low-CV P-unit, where nonlinearities in the P-unit response occurred at the sum and difference of the beat frequencies for pure sine-wave stimulation (\figrefb{fig:motivation}). In this P-unit the nonlinearity appeared in a three-wave setting with \bone{} being close to \fbase{}, corresponding to a frequency combination on the vertical line at \foneb{} in the second-order susceptibility (\subfigrefb{fig:model_full}{B}). This implies that even if only the diagonal structure can be accessed with noise stimulation in the second-order susceptibility in an experiment (\subfigrefb{fig:model_full}{A}), it can be taken as an indicator that the whole nonlinear structure should be present during pure sine-wave stimulation in the same cell. With this RAM stimulation was demonstrated to be an effective method to scan the three-fish or two-beat plane and estimate the whole theoretically predicted nonlinear structure in experimentally recorded cells.
\subsection{Selective readout versus integration of heterogeneous populations}% Nonlinearity might be influenced once integrating from a P-unit population with heterogeneous baseline properties}%Heterogeneity of P-units might influence
The observed nonlinear effects might facilitate the detectability of faint signals during a three fish setting, the electrosensory cocktail party. These nonlinear effects are, however, very specific with respect to the relation of stimulus frequencies and the P-unit baseline frequency. The EOD frequencies of the interacting fish would be drawn from the distributions of EOD frequencies found in male and female fish\cite{Hopkins1974Eigen, Meyer1987, Henninger2018, Henninger2020}. To be behaviorally relevant the faint signal detection would require reliable neuronal signaling irrespective of the individual EOD frequencies.
P-units, however are very heterogeneous in their baseline firing properties\cite{Grewe2017, Hladnik2023}. The baseline firing rates vary in wide ranges (50--450\,Hz). This range covers substantial parts of the beat frequencies that may occur during animal interactions which is limited to frequencies below the Nyquist frequency (0 -- \feod/2)\cite{Barayeu2023}. It is thus likely that there are P-units which match approximately to the specificities of the different encounters.
The observed nonlinear effects might facilitate the detectability of faint signals during a three fish setting, the electrosensory cocktail party. These nonlinear effects are, however, very specific with respect to the relation of stimulus frequencies and the P-unit baseline frequency. The EOD frequencies of the interacting fish would be drawn from the distributions of EOD frequencies found in male and female fish \citep{Hopkins1974Eigen, Meyer1987, Henninger2018, Henninger2020}. To be behaviorally relevant the faint signal detection would require reliable neuronal signaling irrespective of the individual EOD frequencies.
P-units, however are very heterogeneous in their baseline firing properties \citep{Grewe2017, Hladnik2023}. The baseline firing rates vary in wide ranges (50--450\,Hz). This range covers substantial parts of the beat frequencies that may occur during animal interactions which is limited to frequencies below the Nyquist frequency (0 -- \feod/2) \citep{Barayeu2023}. It is thus likely that there are P-units which match approximately to the specificities of the different encounters.
On the other hand, the nonlinearity was found only in low-CV P-units (with white noise stimulation). The CVs are also very heterogeneous (0.1--1.4, \figref{fig:data_overview}\panel{A}) in our sample. Only a small fraction of the P-units has a sufficiently low level of intrinsic noise and will exhibit nonlinear responses. The P-units project to the ELL\cite{Krahe2014} and the integrating pyramidal cells in the different segments receive inputs in the range of 10s to 1000s of neurons\cite{Maler2009a}. Since the receptive fields of the pyramidal neurons are patches of adjacent receptors on the fish's body surface \cite{Bastian2002, Maler2009a, Haggard2023} and the input heterogeneity does not depend on the location of the receptor on the fish body \cite{Hladnik2023} the pyramidal cell input population will be heterogeneous. Heterogeneity was shown to be generally advantageous for the encoding in this\cite{Hladnik2023} and other systems\cite{Padmanabhan2010, Beiran2018}. At the same time it contradicts the apparent need for a selective readout of low-CV cells to maintain information arising through nonlinear interactions.
On the other hand, the nonlinearity was found only in low-CV P-units (with white noise stimulation). The CVs are also very heterogeneous (0.1--1.4, \figref{fig:data_overview}\panel{A}) in our sample. Only a small fraction of the P-units has a sufficiently low level of intrinsic noise and will exhibit nonlinear responses. The P-units project to the ELL \citep{Krahe2014} and the integrating pyramidal cells in the different segments receive inputs in the range of 10s to 1000s of neurons \citep{Maler2009a}. Since the receptive fields of the pyramidal neurons are patches of adjacent receptors on the fish's body surface \citep{Bastian2002, Maler2009a, Haggard2023} and the input heterogeneity does not depend on the location of the receptor on the fish body \citep{Hladnik2023} the pyramidal cell input population will be heterogeneous. Heterogeneity was shown to be generally advantageous for the encoding in this \citep{Hladnik2023} and other systems \citep{Padmanabhan2010, Beiran2018}. At the same time it contradicts the apparent need for a selective readout of low-CV cells to maintain information arising through nonlinear interactions.
A possible readout mechanism should be the topic of future studies that need to take into account that the nonlinearities are stronger in pure sine-wave stimulation and the fraction of cells that show it under naturalistic stimulation might be larger than expected from the distribution of CVs.
\subsection{Behavioral relevance of nonlinear interactions}
The behavioral relevance of the weak signal detection in P-units is evident from the courtship context observed in freely interacting animals\cite{Henninger2018}. Outside courtship behavior, the encoding of secondary or social envelopes is a common need\cite{Stamper2012Envelope}. In a previous study it was demonstrated that information about low-frequency secondary envelopes would not be present in P-units' responses but would arise thorough nonlinear processing downstream in the ELL \cite{Middleton2006,Middleton2007}. Based on our work we would predict that only a small subset of cells, with low CVs, should encode the social envelopes under white noise stimulation. An absence of low-CVs cells in the population analyzed in the previous studies could explain their conclusions. On the other hand, another study showed that P-units with strong nonlinearities, low firing rates and high CVs could indeed encode social envelopes\cite{Savard2011}. These findings are in contrast to the previously mentioned work\cite{Middleton2007} and, at first glance, also to our results. The missing link, that has not been considered in this work, might be the bursting of P-units, the repeated firing of spikes after one EOD period interleaved with longer intervals of quietness\cite{Chacron2004}. Bursting was not explicitly addressed in the previous work, still the reported high CVs of the envelope encoding P-units indicate a higher rate of bursting\cite{Savard2011}. How bursts influence the second-order susceptibility of P-units will be addressed in following works (in preparation). Note that in this work we operated in a regime of weak stimuli and that the envelope encoding addressed in \cite{Savard2011,Middleton2007} operates in a regime of strong stimuli, where the firing rate is saturated. The exact transition from the nonlinearities in a regime of weak stimuli to a regime of strong stimuli could be addressed in further P-unit studies.
The behavioral relevance of the weak signal detection in P-units is evident from the courtship context observed in freely interacting animals \citep{Henninger2018}. Outside courtship behavior, the encoding of secondary or social envelopes is a common need \citep{Stamper2012Envelope}. In a previous study it was demonstrated that information about low-frequency secondary envelopes would not be present in P-units' responses but would arise thorough nonlinear processing downstream in the ELL \citep{Middleton2006,Middleton2007}. Based on our work we would predict that only a small subset of cells, with low CVs, should encode the social envelopes under white noise stimulation. An absence of low-CVs cells in the population analyzed in the previous studies could explain their conclusions. On the other hand, another study showed that P-units with strong nonlinearities, low firing rates and high CVs could indeed encode social envelopes \citep{Savard2011}. These findings are in contrast to the previously mentioned work \citep{Middleton2007} and, at first glance, also to our results. The missing link, that has not been considered in this work, might be the bursting of P-units, the repeated firing of spikes after one EOD period interleaved with longer intervals of quietness \citep{Chacron2004}. Bursting was not explicitly addressed in the previous work, still the reported high CVs of the envelope encoding P-units indicate a higher rate of bursting \citep{Savard2011}. How bursts influence the second-order susceptibility of P-units will be addressed in following works (in preparation). Note that in this work we operated in a regime of weak stimuli and that the envelope encoding addressed in \citep{Savard2011,Middleton2007} operates in a regime of strong stimuli, where the firing rate is saturated. The exact transition from the nonlinearities in a regime of weak stimuli to a regime of strong stimuli could be addressed in further P-unit studies.
Sinusoidal AMs are relevant in interactions with a few fish. We can understand the noise as the presence of many animals with individual EOD frequencies at the same time. Under noise stimulation, nonlinearities were demonstrated to be strong for weak stimuli but were shown to decrease for stronger noise stimuli (\figrefb{fig:cells_suscept}). As long as the noise signal is weak, those fish are distant and the nonlinearity is maintained. An increasing stimulus amplitude would indicate that many fish are close to the receiver and a decrease of nonlinear effects can be observed. These findings imply that the nonlinear effects arising in the presence of three fish decline the more fish join. \lepto{} usually prefers small groups of fish\cite{Stamper2010}. Thus, the described second-order susceptibility might still be behaviorally relevant under natural conditions. The decline of nonlinear effects when several fish are present might be an adaptive process reducing the number of frequencies represented in its primary sensory afferents to a minimum. Such representation would still leave room to create nonlinear effects at later processing steps in higher-order neurons.
Sinusoidal AMs are relevant in interactions with a few fish. We can understand the noise as the presence of many animals with individual EOD frequencies at the same time. Under noise stimulation, nonlinearities were demonstrated to be strong for weak stimuli but were shown to decrease for stronger noise stimuli (\figrefb{fig:cells_suscept}). As long as the noise signal is weak, those fish are distant and the nonlinearity is maintained. An increasing stimulus amplitude would indicate that many fish are close to the receiver and a decrease of nonlinear effects can be observed. These findings imply that the nonlinear effects arising in the presence of three fish decline the more fish join. \lepto{} usually prefers small groups of fish \citep{Stamper2010}. Thus, the described second-order susceptibility might still be behaviorally relevant under natural conditions. The decline of nonlinear effects when several fish are present might be an adaptive process reducing the number of frequencies represented in its primary sensory afferents to a minimum. Such representation would still leave room to create nonlinear effects at later processing steps in higher-order neurons.
The afferents of the passive electrosensory system, the ampullary cells, were found to exhibit much stronger nonlinearities than P-units (\figref{fig:data_overview}). The adequate stimulus for this system is a direct stimulation not an amplitude modulation. In this sense the ampullary cells are closer to the LIF models used by Voroneko and colleagues\cite{Voronenko2017} and we can thus expect that the same nonlinear mechanisms are at work here. For the ampullary system, the sinewave stimulation is not as relevant as for the P-unit system. Ampullary cells encode low-frequency exogenous electric signals such as muscle potentials induced by prey movement\cite{Kalmijn1974, Engelmann2010, Neiman2011fish}. The simultaneous muscle activity of a swarm of prey (such as \textit{Daphnia}) resembles Gaussian white noise\cite{Neiman2011fish}, similar to the stimuli used here. Our results show some similarities with the analyses by Neiman and Russel\cite{Neiman2011fish} who study the encoding in ampullary afferents in the paddlefish. There, the power spectrum of the spontaneous activity also shows a peak at the baseline frequency (internal oscillator) but also at the oscillation frequency measured at the epithelium and interactions of both. Most of the latter disappear in the case of external stimulation, though. Here we find only peaks at the baseline frequency of the neuron and its harmonics. There are interesting similarities and dissimilarities; stimulus encoding in the paddlefish as well as in the brown ghost is very linear for low frequencies and there are nonlinearities in both systems. Linear encoding in the paddlefish shows a gap in the spectrum at the frequency of the epitheliar oscillation, instead the nonlinear response is very pronounced there. In \lepto{}, the dominating frequency under baseline conditions is the baseline firing rate, and we do see increased nonlinearity in this frequency range. The baseline frequency, however, is outside the linear coding range\cite{Grewe2017} while it is within the linear coding range in paddlefish\cite{Neiman2011fish}. Interestingly, the nonlinear response in the paddlefish ampullaries increases with stimulus intensity while it disappears in our case (\subfigrefb{fig:data_overview}{F}). The population of ampullary cells is much more homogeneous with respect to the baseline rate (131$\pm$29\,Hz) and stimulus encoding properties than the P-units\cite{Grewe2017}. This implies that, if the stimulus contains the appropriate frequency components that sum up to the baseline rate, there should be a nonlinear response at a frequency that is similar for the full population of ampullary cells (the baseline frequency) that is outside the linear coding range. Postsynaptic cells integrating ampullary input might be able to extract this nonlinear response from the input population. How such nonlinear effects in ampullary cells might influence prey detection should be addressed in further studies.
The afferents of the passive electrosensory system, the ampullary cells, were found to exhibit much stronger nonlinearities than P-units (\figref{fig:data_overview}). The adequate stimulus for this system is a direct stimulation not an amplitude modulation. In this sense the ampullary cells are closer to the LIF models used by Voroneko and colleagues \citep{Voronenko2017} and we can thus expect that the same nonlinear mechanisms are at work here. For the ampullary system, the sinewave stimulation is not as relevant as for the P-unit system. Ampullary cells encode low-frequency exogenous electric signals such as muscle potentials induced by prey movement \citep{Kalmijn1974, Engelmann2010, Neiman2011fish}. The simultaneous muscle activity of a swarm of prey (such as \textit{Daphnia}) resembles Gaussian white noise \citep{Neiman2011fish}, similar to the stimuli used here. Our results show some similarities with the analyses by Neiman and Russel \citep{Neiman2011fish} who study the encoding in ampullary afferents in the paddlefish. There, the power spectrum of the spontaneous activity also shows a peak at the baseline frequency (internal oscillator) but also at the oscillation frequency measured at the epithelium and interactions of both. Most of the latter disappear in the case of external stimulation, though. Here we find only peaks at the baseline frequency of the neuron and its harmonics. There are interesting similarities and dissimilarities; stimulus encoding in the paddlefish as well as in the brown ghost is very linear for low frequencies and there are nonlinearities in both systems. Linear encoding in the paddlefish shows a gap in the spectrum at the frequency of the epitheliar oscillation, instead the nonlinear response is very pronounced there. In \lepto{}, the dominating frequency under baseline conditions is the baseline firing rate, and we do see increased nonlinearity in this frequency range. The baseline frequency, however, is outside the linear coding range \citep{Grewe2017} while it is within the linear coding range in paddlefish \citep{Neiman2011fish}. Interestingly, the nonlinear response in the paddlefish ampullaries increases with stimulus intensity while it disappears in our case (\subfigrefb{fig:data_overview}{F}). The population of ampullary cells is much more homogeneous with respect to the baseline rate (131$\pm$29\,Hz) and stimulus encoding properties than the P-units \citep{Grewe2017}. This implies that, if the stimulus contains the appropriate frequency components that sum up to the baseline rate, there should be a nonlinear response at a frequency that is similar for the full population of ampullary cells (the baseline frequency) that is outside the linear coding range. Postsynaptic cells integrating ampullary input might be able to extract this nonlinear response from the input population. How such nonlinear effects in ampullary cells might influence prey detection should be addressed in further studies.
\subsection{Conclusion}
We have demonstrated that there are pronounced nonlinear responses in the primary electrosensory afferences of the weakly electric fish \lepto{}, systems that are very often characterized using linear methods. The observed nonlinearities match the expectations from previous theoretical studies\cite{Voronenko2017}. We can confirm that the theory applies also to systems that are encoding amplitude modulations of a carrier signal. Comparisons of P-units and ampullary cells showed that it is the level of intrinsic noise that determines how strongly nonlinear the system acts. Using the second-oder susceptibility estimated from the responses to white noise stimuli provides an easy way to determine the nonlinearity of the system under study. P-units share several features with mammalian
auditory nerve fibers and such nonlinear effects might also be expected in the auditory system during the encoding of amplitude modulations \cite{Joris2004}.
We have demonstrated that there are pronounced nonlinear responses in the primary electrosensory afferences of the weakly electric fish \lepto{}, systems that are very often characterized using linear methods. The observed nonlinearities match the expectations from previous theoretical studies \citep{Voronenko2017}. We can confirm that the theory applies also to systems that are encoding amplitude modulations of a carrier signal. Comparisons of P-units and ampullary cells showed that it is the level of intrinsic noise that determines how strongly nonlinear the system acts. Using the second-oder susceptibility estimated from the responses to white noise stimuli provides an easy way to determine the nonlinearity of the system under study. P-units share several features with mammalian
auditory nerve fibers and such nonlinear effects might also be expected in the auditory system during the encoding of amplitude modulations \citep{Joris2004}.
\section{Methods}
@ -612,7 +610,7 @@ Respiration was then switched to normal tank water and the fish was transferred
For the recordings fish were positioned centrally in the experimental tank, with the major parts of their body submerged into the water. Those body parts that were above the water surface were covered with paper tissue to avoid drying of the skin. Local analgesia was refreshed in intervals of two hours by cutaneous reapplication of Lidocaine (2\,\%; bela-pharm, Vechta, Germany) around the surgical wounds. Electrodes (borosilicate; 1.5\,mm outer diameter; GB150F-8P; Science Products, Hofheim, Germany) were pulled to a resistance of 50--100\,\mega\ohm{} (model P-97; Sutter Instrument, Novato, CA) and filled with 1\,M KCl solution. Electrodes were fixed in a microdrive (Luigs-Neumann, Ratingen, Germany) and lowered into the nerve (\figrefb{fig:setup}, blue triangle). Recordings of electroreceptor afferents were amplified and lowpass filtered at 10\,kHz (SEC-05, npi-electronics, Tamm, Germany, operated in bridge mode). All signals, neuronal recordings, recorded EOD and the generated stimulus, were digitized with sampling rates of 20 or 40\,kHz (PCI-6229, National Instruments, Austin, TX). RELACS (\url{www.relacs.net}) running on a Linux computer was used for online spike and EOD detection, stimulus generation, and calibration. Recorded data was then stored on the hard drive for offline analysis.
\subsection{Identification of P-units and ampullary cells}
The neurons were classified into cell types during the recording by the experimenter. P-units were classified based on baseline firing rates of 50--450\,Hz and a clear phase-locking to the EOD and their responses to amplitude modulations of their own EOD\cite{Grewe2017, Hladnik2023}. Ampullary cells were classified based on firing rates of 80--200\,Hz absent phase-locking to the EOD and responses to low-frequency sinusoidal stimuli\cite{Grewe2017}. We here selected only those cells of which the neuron's baseline activity as well as the responses to frozen noise stimuli were recorded.
The neurons were classified into cell types during the recording by the experimenter. P-units were classified based on baseline firing rates of 50--450\,Hz and a clear phase-locking to the EOD and their responses to amplitude modulations of their own EOD \citep{Grewe2017, Hladnik2023}. Ampullary cells were classified based on firing rates of 80--200\,Hz absent phase-locking to the EOD and responses to low-frequency sinusoidal stimuli \citep{Grewe2017}. We here selected only those cells of which the neuron's baseline activity as well as the responses to frozen noise stimuli were recorded.
\subsection{Electric field recordings}
The electric field of the fish was recorded in two ways: 1. we measured the so-called \textit{global EOD} with two vertical carbon rods ($11\,\centi\meter$ long, 8\,mm diameter) in a head-tail configuration (\figrefb{fig:setup}, green bars). The electrodes were placed isopotential to the stimulus. This signal was differentially amplified with a factor between 100 and 500 (depending on the recorded animal) and band-pass filtered (3 to 1500\,Hz pass-band, DPA2-FX; npi electronics, Tamm, Germany). 2. The so-called \textit{local EOD} was measured with 1\,cm-spaced silver wires located next to the left gill of the fish and orthogonal to the fish's longitudinal body axis (amplification 100 to 500 times, band-pass filtered with 3 to 1\,500\,Hz pass-band, DPA2-FX; npi-electronics, Tamm, Germany, \figrefb{fig:setup}, red markers). This local measurement recorded the combination of the fish's own field and the applied stimulus and thus serves as a proxy of the transdermal potential that drives the electroreceptors.
@ -630,9 +628,9 @@ The fish were stimulated with band-limited white noise stimuli with a cut-off fr
% and between 2.5 and 40\,\% for \eigen
\subsection{Data analysis} Data analysis was performed with Python 3 using the packages matplotlib\cite{Hunter2007}, numpy\cite{Walt2011}, scipy\cite{scipy2020}, pandas\cite{Mckinney2010}, nixio\cite{Stoewer2014}, and thunderfish (\url{https://github.com/bendalab/thunderfish}).
\subsection{Data analysis} Data analysis was performed with Python 3 using the packages matplotlib \citep{Hunter2007}, numpy \citep{Walt2011}, scipy \citep{scipy2020}, pandas \citep{Mckinney2010}, nixio \citep{Stoewer2014}, and thunderfish (\url{https://github.com/bendalab/thunderfish}).
%sklearn\cite{scikitlearn2011},
%sklearn \citep{scikitlearn2011},
\paragraph{Baseline analysis}\label{baselinemethods}
The baseline firing rate \fbase{} was calculated as the number of spikes divided by the duration of the baseline recording (on average 18\,s). The coefficient of variation (CV) was calculated as the standard deviation of the interspike intervals (ISI) divided by the average ISI: $\rm{CV} = \sqrt{\langle (ISI- \langle ISI \rangle) ^2 \rangle} / \langle ISI \rangle$. If the baseline was recorded several times in a recording, the average \fbase{} and CV were calculated.
@ -717,7 +715,7 @@ If the same frozen noise was recorded several times in a cell, each noise repeti
\subsection{Leaky integrate-and-fire models}\label{lifmethods}
Leaky integrate-and-fire (LIF) models with a carrier were constructed to reproduce the specific firing properties of P-units \cite{Chacron2001,Sinz2020}. The sole driving input into the P-unit model during baseline, i.e. when no external stimulus was given, is the fish's own EOD modeled as a cosine wave
Leaky integrate-and-fire (LIF) models with a carrier were constructed to reproduce the specific firing properties of P-units \citep{Chacron2001,Sinz2020}. The sole driving input into the P-unit model during baseline, i.e. when no external stimulus was given, is the fish's own EOD modeled as a cosine wave
\begin{equation}
\label{eq:eod}
\carrierinput = y_{EOD}(t) = \cos(2\pi f_{EOD} t)
@ -763,7 +761,7 @@ Whenever the membrane voltage $V_m(t)$ crossed the spiking threshold $\theta=1$
% \label{eifnl}
% f(V_m)= \Delta_V \text{e}^{\frac{V_m-1}{\Delta_V}}
% \end{equation}
% \cite{Fourcaud-Trocme2003}, where $\Delta_V$ was varied from 0.001 to 0.1.
% \citep{Fourcaud-Trocme2003}, where $\Delta_V$ was varied from 0.001 to 0.1.
%, \figrefb{eif}
\begin{figure*}[t]
@ -797,7 +795,7 @@ From each simulation run, the first second was discarded and the analysis was ba
% The second-order susceptibility in the model was calculated with \Eqnref{eq:susceptibility}, resulting in matrices as in \figrefb{model_and_data} and \figrefb{fig:model_full}. For this, the model neuron was presented the input $x(t)$ for 2\,s, with the first second being dismissed as the transient. The second-order susceptibility calculation was performed on the last second, resulting in a frequency resolution of 1\,Hz.
\subsection{Model noise split into a noise and a stimulus component}\label{intrinsicsplit_methods}% the Furutsu-Novikov Theorem with the same correlation function
According to previous works \cite{Lindner2022} the total noise of a LIF model ($\xi$) can be split up into several independent noise processes. Here we split the internal noise into two parts: (i) One part is treated as a driving input signal $s_\xi(t)$, a RAM stimulus where frequencies above 300\,Hz are discarded (\Eqnref{eq:ram_split}), and used to calculate the cross-spectra in \Eqnref{eq:crosshigh} and (ii) the remaining noise $\sqrt{2D \, c_{\rm{noise}}} \cdot \xi(t)$ that is treated as pure noise (\Eqnref{eq:Noise_split_intrinsic}). In this way the effective signal-to-noise ratio can be increased while maintaining the total noise in the system.
According to previous works \citep{Lindner2022} the total noise of a LIF model ($\xi$) can be split up into several independent noise processes. Here we split the internal noise into two parts: (i) One part is treated as a driving input signal $s_\xi(t)$, a RAM stimulus where frequencies above 300\,Hz are discarded (\Eqnref{eq:ram_split}), and used to calculate the cross-spectra in \Eqnref{eq:crosshigh} and (ii) the remaining noise $\sqrt{2D \, c_{\rm{noise}}} \cdot \xi(t)$ that is treated as pure noise (\Eqnref{eq:Noise_split_intrinsic}). In this way the effective signal-to-noise ratio can be increased while maintaining the total noise in the system.
%\sqrt{\rho \, 2D \,c_{\rm{signal}}} \cdot \xi(t)
%(1-c_{\rm{signal}})\cdot\xi$c_{\rm{noise}} = 1-c_{\rm{signal}}$
@ -826,17 +824,17 @@ According to previous works \cite{Lindner2022} the total noise of a LIF model ($
In the here used model a small portion of the original noise was assigned to the noise component ($c_{\rm{noise}} = 0.1$, \subfigrefb{flowchart}\,\panel[iii]{C}) and a big portion used as the signal component ($c_{\rm{signal}} = 0.9$, \subfigrefb{flowchart}\,\panel[iii]{A}). For the noise split to be valid \cite{Lindner2022} it is critical that both components add up to the initial 100\,\% of the total noise and the baseline properties as the firing rate and the CV of the model are maintained. This is easily achieved in a model without a carrier if the condition $c_{\rm{signal}}+c_{\rm{noise}}=1$ is satisfied. The situation here is more complicated. After the original noise was split into a signal component with $c_{\rm{signal}}$, the frequencies above 300\,Hz were discarded and the signal strength was reduced after the dendritic low pass filtering. To compensate for these transformations the initial signal component was multiplying with the factor $\rho$, keeping the baseline CV (only carrier) and the CV during the noise split comparable, and resulting in $s_\xi(t)$. $\rho$ was found by bisecting the plane of possible factors and minimizing the error between the CV during baseline and stimulation.
In the here used model a small portion of the original noise was assigned to the noise component ($c_{\rm{noise}} = 0.1$, \subfigrefb{flowchart}\,\panel[iii]{C}) and a big portion used as the signal component ($c_{\rm{signal}} = 0.9$, \subfigrefb{flowchart}\,\panel[iii]{A}). For the noise split to be valid \citep{Lindner2022} it is critical that both components add up to the initial 100\,\% of the total noise and the baseline properties as the firing rate and the CV of the model are maintained. This is easily achieved in a model without a carrier if the condition $c_{\rm{signal}}+c_{\rm{noise}}=1$ is satisfied. The situation here is more complicated. After the original noise was split into a signal component with $c_{\rm{signal}}$, the frequencies above 300\,Hz were discarded and the signal strength was reduced after the dendritic low pass filtering. To compensate for these transformations the initial signal component was multiplying with the factor $\rho$, keeping the baseline CV (only carrier) and the CV during the noise split comparable, and resulting in $s_\xi(t)$. $\rho$ was found by bisecting the plane of possible factors and minimizing the error between the CV during baseline and stimulation.
%that was found by minimizing the error between the
%Furutsu-Novikov Theorem \cite{Novikov1965, Furutsu1963}\Eqnref{eq:ram_split}, (red in \subfigrefb{flowchart}\,\panel[iii]{A}) bisecting the space of possible $\rho$ scaling factors
%Furutsu-Novikov Theorem \citep{Novikov1965, Furutsu1963}\Eqnref{eq:ram_split}, (red in \subfigrefb{flowchart}\,\panel[iii]{A}) bisecting the space of possible $\rho$ scaling factors
%$\rho$ a scaling factor that compensates (see below) for the signal transformations the amplitude modulation stimulus undergoes in the model, i.e. the threshold and the dendritic lowpass.
%In our case the model has a carrier (the fish's self-generated EOD) and we thus want to drive the model with an amplitude modulation stimulus
% See section \ref{lifmethods} for model and parameter description.
\begin{table*}[hp!]
\caption{\label{modelparams} Model parameters of LIF models, fitted to 2 electrophysiologically recorded P-units \cite{Ott2020}.}
\caption{\label{modelparams} Model parameters of LIF models, fitted to 2 electrophysiologically recorded P-units \citep{Ott2020}.}
\begin{center}
\begin{tabular}{lrrrrrrrr}
\hline
@ -874,7 +872,7 @@ In the here used model a small portion of the original noise was assigned to the
\section{Supporting information}
\subsection{S1 Second-order susceptibility of high-CV P-unit}
CVs in P-units can range up to 1.5 \cite{Grewe2017, Hladnik2023}. We show the same analysis as in \figrefb{fig:cells_suscept} for an example higher-CV P-unit. Similar to the low-CV cell, high-CV P-units fire at multiples of the EOD period (\subfigrefb{fig:cells_suscept_high_CV}{A}). In contrast to low-CV P-units, however, the higher CV characterizes the noisier baseline firing pattern and the peak at \fbase{} is less pronounced in the power spectrum of the baseline activity (\subfigrefb{fig:cells_suscept_high_CV}{B}). High-CV P-units do not exhibit a clear nonlinear structure related to \fbase{} neither in the second-order susceptibility matrices (\subfigrefb{fig:cells_suscept_high_CV}{E--F}), nor in the projected diagonals (\subfigrefb{fig:cells_suscept_high_CV}{G}). The overall level of nonlinearity, however shows the same dependence on the stimulus contrast. It is much reduced for high-contrast stimuli that drive the neuron much stronger (\subfigrefb{fig:cells_suscept_high_CV}{F}).
CVs in P-units can range up to 1.5 \citep{Grewe2017, Hladnik2023}. We show the same analysis as in \figrefb{fig:cells_suscept} for an example higher-CV P-unit. Similar to the low-CV cell, high-CV P-units fire at multiples of the EOD period (\subfigrefb{fig:cells_suscept_high_CV}{A}). In contrast to low-CV P-units, however, the higher CV characterizes the noisier baseline firing pattern and the peak at \fbase{} is less pronounced in the power spectrum of the baseline activity (\subfigrefb{fig:cells_suscept_high_CV}{B}). High-CV P-units do not exhibit a clear nonlinear structure related to \fbase{} neither in the second-order susceptibility matrices (\subfigrefb{fig:cells_suscept_high_CV}{E--F}), nor in the projected diagonals (\subfigrefb{fig:cells_suscept_high_CV}{G}). The overall level of nonlinearity, however shows the same dependence on the stimulus contrast. It is much reduced for high-contrast stimuli that drive the neuron much stronger (\subfigrefb{fig:cells_suscept_high_CV}{F}).
\label{S1:highcvpunit}
\begin{figure*}[!ht]
@ -883,12 +881,12 @@ CVs in P-units can range up to 1.5 \cite{Grewe2017, Hladnik2023}. We show the sa
\end{figure*}
\begin{figure*}[hp]%hp!
\begin{figure*}[t]
\includegraphics[width=\columnwidth]{trialnr}
\notejb{ylabel: absolute value of chi_2!}
\notejb{ylabel: absolute value of $chi_2$!}
\notejb{sort the legend labels!}
\notejb{one of the 99.99\% percentiles is enough}
\notejb{add 10^7 tick on x-axis}
\notejb{add $10^7$ tick on x-axis}
\caption{\label{fig:trialnr} Dependence of the estimate of the second-order susceptibility on the number of trials $\n$. While the estimate of the noise floor (10th and 90th percentile) of the $|\chi_2(f_1, f_2)|$ matrix does not saturate yet, the estimates of the high values in the matrix that make up the characteristic ridges saturate for $N>10^6$.
}
\end{figure*}
@ -899,24 +897,24 @@ CVs in P-units can range up to 1.5 \cite{Grewe2017, Hladnik2023}. We show the sa
\begin{itemize}
\item \notejb{\cite{French1973} Derivation of the Fourier transformed kernels measured with white noise.}
\item \notejb{\cite{French1976} Technical issues and tests of Fourier transformed kernels measured with white noise.}
\item \notejb{\cite{Victor1977} Cat retinal ganglion cells, gratings with sum of 6 or 8 sinusoids. X - versus Y cells. Peak at f1 == f2 in Y cells. X-cells rather linear. Discussion of mechanism, where a nonlinearity comes in along the pathway}
\item \notejb{\cite{Marmarelis1972} Temporal 2nd order kernels, how well do kernels predict responses, catfish retinal ganglion cells}
\item \notejb{\cite{Marmarelis1973} Temporal 2nd order kernels, how well do kernels predict responses}
\item \notejb{\cite{Victor1988} Cat retinal ganglion cells, sum of sinusoids, very technical, one measurement similar to \cite{Victor1977}.}
\item \notejb{ \citep{French1973} Derivation of the Fourier transformed kernels measured with white noise.}
\item \notejb{ \citep{French1976} Technical issues and tests of Fourier transformed kernels measured with white noise.}
\item \notejb{ \citep{Victor1977} Cat retinal ganglion cells, gratings with sum of 6 or 8 sinusoids. X - versus Y cells. Peak at f1 == f2 in Y cells. X-cells rather linear. Discussion of mechanism, where a nonlinearity comes in along the pathway}
\item \notejb{ \citep{Marmarelis1972} Temporal 2nd order kernels, how well do kernels predict responses, catfish retinal ganglion cells}
\item \notejb{ \citep{Marmarelis1973} Temporal 2nd order kernels, how well do kernels predict responses}
\item \notejb{ \citep{Victor1988} Cat retinal ganglion cells, sum of sinusoids, very technical, one measurement similar to \citep{Victor1977}.}
\item \notejb{\citep{Nikias1993} Third order spectra or bispectra. Very technical overview to higher order spectra}
\item \notejb{\cite{Mitsis2007} Spider mechanorecptor. Linear filterss, multivariate nonlinearity, and threshold. Second order kernel needed for this. Gausian noise stimuli.}
\item \notejb{\cite{French2001} Time kernels up to 3rd order for predicting spider mechanorecptor responses (spikes!)}
\item \notejb{\cite{French1999} Review on time domain nonlinear systems identification}
\item \notejb{\cite{Temchin2005,RecioSpinosa2005} 2nd order Wiener kernel for predicting chinchilla auditory nerve fibre firing rate responses. Strong 2nd order blob at characteristic frequency}
\item \notejb{\cite{Schanze1997} lots of bispectra, visual cortex MUA recordings}
\item \notejb{\cite{Theunissen1996} Linear backward stimulus reconstruction in the context of information theory/signal-to-noise ratios}
\item \notejb{\cite{Wessel1996} Same as Theunissen1996 but for P-units}
\item \notejb{\cite{Neiman2011} cross bispectrum, bicoherence, mutual information, saturating nonlinearities, `` ampullary electroreceptors of paddlefish are perfectly suited to linearly encode weak low-frequency stimuli.''}
\item \notejb{\cite{Chichilnisky2001} Linear Nonlinear Poisson model}
\item \notejb{\cite{Gollisch2009} Linear Nonlinear models in retina}
\item \notejb{\cite{Clemens2013} Grasshoppper model for female preferences}
\item \notejb{ \citep{Mitsis2007} Spider mechanorecptor. Linear filterss, multivariate nonlinearity, and threshold. Second order kernel needed for this. Gausian noise stimuli.}
\item \notejb{ \citep{French2001} Time kernels up to 3rd order for predicting spider mechanorecptor responses (spikes!)}
\item \notejb{ \citep{French1999} Review on time domain nonlinear systems identification}
\item \notejb{ \citep{Temchin2005,RecioSpinosa2005} 2nd order Wiener kernel for predicting chinchilla auditory nerve fibre firing rate responses. Strong 2nd order blob at characteristic frequency}
\item \notejb{ \citep{Schanze1997} lots of bispectra, visual cortex MUA recordings}
\item \notejb{ \citep{Theunissen1996} Linear backward stimulus reconstruction in the context of information theory/signal-to-noise ratios}
\item \notejb{ \citep{Wessel1996} Same as Theunissen1996 but for P-units}
\item \notejb{ \citep{Neiman2011} cross bispectrum, bicoherence, mutual information, saturating nonlinearities, `` ampullary electroreceptors of paddlefish are perfectly suited to linearly encode weak low-frequency stimuli.''}
\item \notejb{ \citep{Chichilnisky2001} Linear Nonlinear Poisson model}
\item \notejb{ \citep{Gollisch2009} Linear Nonlinear models in retina}
\item \notejb{ \citep{Clemens2013} Grasshoppper model for female preferences}
\end{itemize}