jg edits and some comments

nonlinearbaseline.tex

@@ -230,7 +230,6 @@
 \newcommand{\notebl}[1]{\note[BL]{#1}}
 \newcommand{\notems}[1]{\note[MS]{#1}}
 
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{document}
@@ -252,8 +251,9 @@
 
 % Please keep the abstract below 300 words
 \section{Abstract}
-Neuronal processing is inherently nonlinear --- spiking thresholds or rectification at synapses is essential to neuronal computations. Nevertheless, linear response theory has been instrumental in understanding, for example, the impact of noise or neuronal synchrony on signal transmission, or the emergence of oscillatory activity, but is valid only at low stimulus amplitudes or large levels of intrinsic noise. At higher signal-to-noise ratios, however, nonlinear response components become relevant. Theoretical results for leaky integrate-and-fire neurons in the weakly nonlinear regime suggest strong responses at the sum of two input frequencies if these frequencies or their sum match the neuron's baseline firing rate.
+%Neuronal processing is inherently nonlinear --- spiking thresholds or rectification at synapses are essential to neuronal computations.}
-We here analyze nonlinear responses in two types of primary electroreceptor afferents, the P-units of the active and the ampullary cells of the passive electrosensory system of the wave-type electric fish \textit{Apteronotus leptorhynchus}. In our combined experimental and modeling approach we identify the predicted nonlinear responses in low-noise P-units and much stronger in more than half of the ampullary cells. Our results provide experimental evidence for nonlinear responses of spike generators in the weakly nonlinear regime. We conclude that such nonlinear responses occur in any sensory neuron that operates in similar regimes particularly at near-threshold stimulus conditions.
+Spiking thresholds in neurons or rectification at synapses are essential for neuronal computations rendering neuronal processing inherently nonlinear. Nevertheless, linear response theory has been instrumental for understanding, for example, the impact of noise or neuronal synchrony on signal transmission, or the emergence of oscillatory activity, but is valid only at low stimulus amplitudes or large levels of intrinsic noise. At higher signal-to-noise ratios, however, nonlinear response components become relevant. Theoretical results for leaky integrate-and-fire neurons in the weakly nonlinear regime suggest strong responses at the sum of two input frequencies if these frequencies or their sum match the neuron's baseline firing rate.
+We here analyze nonlinear responses in two types of primary electroreceptor afferents, the P-units of the active and the ampullary cells of the passive electrosensory system of the wave-type electric fish \textit{Apteronotus leptorhynchus}. In our combined experimental and modeling approach we identify these predicted nonlinear responses in low-noise P-units and, much stronger, in ampullary cells. Our results provide experimental evidence for nonlinear responses of spike generators in the weakly nonlinear regime. We conclude that such nonlinear responses occur in any sensory neuron that operates in similar regimes particularly at near-threshold stimulus conditions.
 
 % Please keep the Author Summary between 150 and 200 words
 % Use the first person. PLOS ONE authors please skip this step. 
@@ -268,7 +268,8 @@ We here analyze nonlinear responses in two types of primary electroreceptor affe
   \includegraphics[width=\columnwidth]{lifsuscept}
   \caption{\label{fig:lifresponse} First- and second-order response functions of the leaky integrate-and-fire model. \figitem{A} Magnitude of the first-order (linear) response function $|\chi_1(f)|$, also known as the ``gain'' function, quantifies the response amplitude relative to the stimulus amplitude, both measured at the same stimulus frequency. \figitem{B} Magnitude of the second-order (non-linear) 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 function is zero, because linear systems do not create new frequencies and thus there is no response at the sum of the two frequencies. The plots show the analytical solutions from \citet{Lindner2001} and \citet{Voronenko2017} with $\mu = 1.1$ and $D = 0.001$. Note that the leaky integrate-and-fire model is formulated without dimensions, frequencies are given in multiples of the inverse membrane time constant.}
 \end{figure*}
-We like to think about signal encoding in terms of linear relations with unique mapping of a given input value to a certain output of the system under consideration. Indeed, such linear methods, for example the transfer function or first-order susceptibility shown in figure~\ref{fig:lifresponse}, have been widely and successfully applied to describe and predict neuronal responses and are an invaluable tool to characterize neural systems \citep{Eggermont1983,Borst1999}. Nonlinear mechanisms, on the other hand, are key on different levels of neural processing. Deciding for one action over another is a 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.
+
+We like to think about signal encoding in terms of linear relations with unique mapping of a given input value to a certain output of the system under consideration. Indeed, such linear methods, for example the transfer function or first-order susceptibility shown in figure~\ref{fig:lifresponse}, have been widely and successfully applied to describe and predict neuronal responses and are an invaluable tools to characterize neural systems \citep{Eggermont1983,Borst1999}. Nonlinear mechanisms, on the other hand, are key on different levels of neural processing. Deciding for one action over another is a 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.
 
 The transfer function that describes the linear properties of a system, is the first-order term of a Volterra series. Higher-order terms successively approximate nonlinear features of a system \citep{Rieke1999}. Second-order kernels have been used in the time domain to predict visual responses in catfish \citep{Marmarelis1972}. In the frequency domain, second-order kernels are known as second-order response functions or susceptibilities. Nonlinear interactions of two stimulus frequencies generate peaks in the response spectrum at the sum and the difference of the two. Including higher-order terms of the Volterra series, the nonlinear nature of mammalian visual systems \citep{Victor1977, Schanze1997}, auditory responses in the Torus semicircularis of frogs \citep{Aertsen1981}, locking in chinchilla auditory nerve fibers \citep{Temchin2005}, spider mechanoreceptors \citep{French2001}, and bursting responses in paddlefish \citep{Neiman2011} have been demonstrated.
 
@@ -284,10 +285,10 @@ Here we search for such weakly nonlinear responses in electroreceptors of the tw
   }
 \end{figure*}
 
-We explored a large set of recordings of P-units and ampullary cells of the active and passive electrosensory systems of the brown ghost knifefish \textit{Apteronotus leptorhynchus} to search for weakly nonlinear responses as they have been predicted by theoretical work \citep{Voronenko2017}. Additional simulations of LIF-based models of P-unit spiking put the experimental findings into context. We start with demonstrating the basic concepts using example P-units and models.
+We explored a large set of electrophysiological data from primary afferents of the active and passive electrosensory system, P-units and ampullary cells, that were recorded in the brown ghost knifefish \textit{Apteronotus leptorhynchus} to search for the weakly nonlinear responses that have been predicted in previous theoretical work\citep{Voronenko2017}. Additional simulations of LIF-based models of P-unit spiking put the experimental findings into context\notejg{very unspecific}. We start with demonstrating the basic concepts using example P-units and respective models and then compare the population of recordings in both cell types.
 
 \subsection{Nonlinear responses in P-units stimulated with two frequencies}
-Without external stimulation, a P-unit is driven by the fish's own EOD alone (with a specific EOD frequency $f_{\rm EOD}$) and spontaneously fires action potentials at the baseline rate $r$. Accordingly, the power spectrum of the baseline activity has a peak at $r$ (\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 - f_{\rm EOD}$ between the two fish. P-units encode this beat in their firing rate \citep{Bastian1981a} and consequently the power spectrum of the response has a peak at the beat frequency (\subfigrefb{fig:motivation}{B}). A second peak at the first harmonic of the beat frequency is indicative of a nonlinear process that here is associated with the clipping of the P-unit's firing rate at zero \citep{Barayeu2023}. Pairing the fish with another fish at a higher beat frequency $\Delta f_2 = f_2 - f_{\rm EOD} > \Delta f_1$ results in a weaker response with a single peak in the response power spectrum, suggesting a linear response (\subfigrefb{fig:motivation}{C}). The weaker response to this beat can be explained by the beat tuning of the cell \citep{Walz2014}. Note, $\Delta f_2$ has been deliberately chosen to match the recorded P-unit's baseline firing rate.
+Without external stimulation, a P-unit is driven by the fish's own EOD alone (with a specific EOD frequency $f_{\rm EOD}$) and spontaneously fires action potentials at the baseline rate $r$. Accordingly, the power spectrum of the baseline activity has a peak at $r$ (\subfigrefb{fig:motivation}{A}). In the communication context, this animal (the receiver) is exposed to the EODs of one or many foreign fish. 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 - f_{\rm EOD}$ between the two fish. P-units encode this beat in their firing rate \citep{Bastian1981a} and consequently the power spectrum of the response has a peak at the beat frequency (\subfigrefb{fig:motivation}{B}). A second peak at the first harmonic of the beat frequency is indicative of a nonlinear process that here is associated with the clipping of the P-unit's firing rate at zero \citep{Barayeu2023}. Pairing the fish with another fish at a higher beat frequency $\Delta f_2 = f_2 - f_{\rm EOD} > \Delta f_1$ results in a weaker response with a single peak in the response power spectrum, suggesting a linear response (\subfigrefb{fig:motivation}{C}). The weaker response to this beat can be explained by the beat tuning of the cell \citep{Walz2014}. Note, $\Delta f_2$ has been deliberately chosen to match the recorded P-unit's baseline firing rate.
 
 When stimulating with both foreign signals simultaneously, additional peaks appear in the response power spectrum at the sum $\Delta f_1 + \Delta f_2$ and the difference frequency $\Delta f_2 - \Delta f_1$ (\subfigrefb{fig:motivation}{D}). Thus, the cellular response is not equal to the sum of the responses to the two beats presented separately. These additional 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.
 
@@ -303,7 +304,7 @@ The stimuli used in \figref{fig:motivation} had the same not-small amplitude. Wh
 
 At very low stimulus contrasts (in the example cell less than approximately 1.2\,\% relative to the receiver's EOD amplitude) the spectrum has small peaks only at the beat frequencies (\subfigref{fig:regimes}{A,B}, green and purple). The amplitudes of these peaks initially increase linearly with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines), an indication of the linear response at lowest stimulus amplitudes.
 
-This linear regime is followed by the weakly nonlinear regime (in the example cell between approximately 1.2\,\% and 3.5\,\% stimulus contrasts). In addition to the peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:regimes}{C}, orange and red). The amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines). Note, that we have chosen $\Delta f_2$ to match the baseline firing rate $f_{base}$ of the neuron.
+This linear regime is followed by the weakly nonlinear regime (in the example cell between approximately 1.2\,\% and 3.5\,\% stimulus contrast). In addition to the peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:regimes}{C}, orange and red). The amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines). Note, that we have chosen $\Delta f_2$ to match the baseline firing rate $f_{base}$ of the neuron.
 
 At higher stimulus amplitudes, the linear response and the weakly-nonlinear response begin to deviate from their linear and quadratic dependency on amplitude (\subfigrefb{fig:regimes}{E}) and additional peaks appear in the response spectrum (\subfigrefb{fig:regimes}{D}). At high stimulus contrasts, additional nonlinearities in the system, in particular clipping of the firing rate, shape the responses.
 
@@ -340,20 +341,21 @@ Electric fish possess an additional electrosensory system, the passive or ampull
 
 
 \subsection{Model-based estimation of the second-order susceptibility}
-In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagonal 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,Franzen2023}. However, a pronounced nonlinear response where one of the stimulus frequencies matches the baseline firing rate (horizontal and vertical ridges), although predicted by theory, cannot be observed. In the following, we investigate how these discrepancies can be understood.
+In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagonal 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,Franzen2023}. However, a pronounced nonlinear response where one of the stimulus frequencies matches the baseline firing rate (horizontal and vertical ridges), although predicted by theory, cannot be observed. In the following, we investigate how these discrepancies can be understood. \notejg{shorter, first sentence is mainly repetition} 
 
 \begin{figure*}[p]
   \includegraphics[width=\columnwidth]{noisesplit}
   \caption{\label{fig:noisesplit} Estimation of second-order susceptibilities. \figitem{A} $|\chi_2(f_1, f_2)|$ (right) estimated from $N=100$ 256\,ms long FFT segments of an electrophysiological recording of another P-unit (cell ``2017-07-18-ai'', $r=78$\,Hz, CV$_{\text{base}}=0.22$) driven with a RAM stimulus with contrast 5\,\% (left). \figitem[i]{B} \textit{Standard condition} of model simulations with intrinsic noise (bottom) and a RAM stimulus (top). \figitem[ii]{B} $|\chi_2(f_1, f_2)|$ estimated from simulations of the cell's LIF model counterpart (cell ``2017-07-18-ai'') based on the same RAM contrast and number of $N=100$ FFT segments. As in the electrophysiological recording only a weak anti-diagonal is visible. \figitem[iii]{B} Same as \panel[ii]{B} but using $10^6$ FFT segments. Now, the expected triangular structure is revealed. \figitem[iv]{B} Convergence of the $|\chi_2(f_1, f_2)|$ estimate as a function of FFT segments. \figitem{C} At a lower stimulus contrast of 1\,\% the estimate did not converge yet even for $10^6$ FFT segments. The triangular structure is not revealed yet. \figitem[i]{D} Same as in \panel[i]{B} but in the \textit{noise split} condition: there is no external RAM signal (red) driving the model. Instead, a large part (90\,\%) of the total intrinsic noise is treated as a signal and is presented as an equivalent amplitude modulation ($s_{\xi}(t)$, orange, 10.6\,\% contrast), while the intrinsic noise is reduced to 10\,\% of its original strength (bottom, see methods for details). \figitem[i]{D} 100 FFT segments are still not sufficient for estimating $|\chi_2(f_1, f_2)|$. \figitem[iii]{D} Simulating one million segments reveals the full expected triangular structure of the second-order susceptibility. \figitem[iv]{D} In the noise-split condition, the $|\chi_2(f_1, f_2)|$ estimate converges already at about $10^{4}$ FFT segments.}
 \end{figure*}
 
-One simple reason could be the lack of data, i.e. the estimation of the second-order susceptibility is not good enough. Electrophysiological recordings are limited in time, and therefore only a limited number of trials, here repeated presentations of the same frozen RAM stimulus, are available. In our data set we have 1 to 199 trials (median: 10) of RAM stimuli with a duration ranging from 2 to 5\,s (median: 8\,s), resulting in a total duration of 30 to 400\,s. Using a resolution of 0.5\,ms and FFT segments of 512 samples this yields 105 to 1520 available FFT segments for a specific RAM stimulus. As a consequence, the cross-spectra, \eqnref{crossxss}, 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, dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of the P-unit estimated from the same low number of FFT (fast fourier transform) segments as in the experiment ($N=100$, compare faint anti-diagonal in the bottom left corner of the second-order susceptibility in \panel[ii]{A} and \panel[ii]{B} in \figrefb{fig:noisesplit}).
+One simple reason could be the lack of data, i.e. the estimation of the second-order susceptibility is not good enough. Electrophysiological recordings are limited in time, and therefore only a limited number of trials, here repeated presentations of the same frozen RAM stimulus, are available. In our data set we have 1 to 199 trials (median: 10) of RAM stimuli with a duration ranging from 2 to 50\,s (median: 10\,s), total stimulation durations per cell range between 30 and 400\,s. Using a temporal resolution of 0.5\,ms and FFT segments of 512 samples this yields 105 to 1520 available FFT segments for a specific RAM stimulus. As a consequence, the cross-spectra, \eqnref{crossxss}, 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, dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of the P-unit estimated from the same low number of FFT (fast fourier transform) segments as in the experiment ($N=100$, compare faint anti-diagonal in the bottom left corner of the second-order susceptibility in \panel[ii]{A} and \panel[ii]{B} in \figrefb{fig:noisesplit}).
 
-In model simulations we can increase the number of FFT segments beyond what would be experimentally possible, here to one million (\figrefb{fig:noisesplit}\,\panel[iii]{B}). Then the estimate of the second-order susceptibility indeed improves. It gets less noisy, the diagonal at $f_ + f_2 = r$ is emphasized, and the vertical and horizontal ridges at $f_1 = r$ and $f_2 = r$ are revealed. Increasing the number of FFT segments also reduces the order of magnitude of the susceptibility estimate until close to one million the estimate levels out at a low values (\subfigrefb{fig:noisesplit}\,\panel[iv]{B}).
+In model simulations we can increase the number of FFT segments beyond what would be experimentally possible, here to one million (\figrefb{fig:noisesplit}\,\panel[iii]{B}). Then, the estimate of the second-order susceptibility indeed improves. It gets less noisy, the diagonal at $f_ + f_2 = r$ is emphasized, and the vertical and horizontal ridges at $f_1 = r$ and $f_2 = r$ are revealed. Increasing the number of FFT segments also reduces the order of magnitude of the susceptibility estimate until close to one million the estimate levels out at a low values (\subfigrefb{fig:noisesplit}\,\panel[iv]{B}).
 
 At a lower stimulus contrast of 1\,\% (\subfigrefb{fig:noisesplit}{C}), however, one million FFT segments are still not sufficient for the estimate to converge (\figrefb{fig:noisesplit}\,\panel[iv]{C}). Still only a faint anti-diagonal is visible (\figrefb{fig:noisesplit}\,\panel[iii]{C}).
 
-Using a broadband stimulus increases the effective input-noise level. 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.  As we just have seen, 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{crossxss} and \eqref{chi2}) \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 a stimulus. This allows us 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 \, \alpha_{\rm noise}}\;\xi(t)$ with $\alpha_{\rm noise} = 0.1$ (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \figrefb{fig:noisesplit}\,\panel[i]{D}). We tuned the amplitude of the RAM stimulus $s_{\xi}(t)$ such that the output firing rate and variability (CV of interspike intervals) 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 $s_\xi(t)$, which is a large fraction of the noise, the signal-to-noise ratio of the measured susceptibilities is drastically improved and thus the estimate converges already at about ten thousand FFT segments (\figrefb{fig:noisesplit}\,\panel[iv]{D}); (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 FFT segments (\figrefb{fig:noisesplit}\,\panel[iii]{D}), but not for a number of segments comparable to the experiment (\figrefb{fig:noisesplit}\,\panel[ii]{D}). In addition to the strong response at $f_1 + f_2 = r$, we now also observe pronounced nonlinear responses at $f_1=r$ and $f_2=r$ (vertical and horizontal lines, \figrefb{fig:noisesplit}\,\panel[iii]{D}).
+Using a broadband stimulus increases the effective input-noise level. 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.  As we just have seen, 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{crossxss} and \eqref{chi2}) \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 a stimulus. This allows us 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 \, \alpha_{\rm noise}}\;\xi(t)$ with $\alpha_{\rm noise} = 0.1$ (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \figrefb{fig:noisesplit}\,\panel[i]{D}). We tuned the amplitude of the RAM stimulus $s_{\xi}(t)$ such that the output firing rate and variability (CV of interspike intervals) 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 $s_\xi(t)$, which is a large fraction of the noise, the signal-to-noise ratio of the measured susceptibilities is drastically improved and thus the estimate converges already at about ten thousand FFT segments (\figrefb{fig:noisesplit}\,\panel[iv]{D}); (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 FFT segments (\figrefb{fig:noisesplit}\,\panel[iii]{D}), but not for a number of segments comparable to the experiment (\figrefb{fig:noisesplit}\,\panel[ii]{D}). In addition to the strong response at $f_1 + f_2 = r$, we now also observe pronounced nonlinear responses at $f_1=r$ and $f_2=r$ (vertical and horizontal lines, \figrefb{fig:noisesplit}\,\panel[iii]{D}).
 
 
 \begin{figure}[p]
@@ -368,11 +370,11 @@ By just looking at the second-order susceptibilities estimated using the noise-s
 
 This categorization is supported by the susceptibility index, SI($r$), \eqnref{siindex}, which quantifies the height of the ridge where the stimulus frequencies add up to the neuron's baseline firing rate relative to the background. Values above one indicate an elevated ridge. The absence of such a ridge results in values close to one. Indeed, the cells showing only a weak triangle (orange) arise out of values around one and the cells showing strong triangles (red) have consistently SI($r$) values exceeding 1.8 (\figrefb{fig:modelsusceptcontrasts}\,\panel[i]{E}).
 
-The SI($r$) correlates with the cell's CV of its baseline interspike intervals ($r=-0.60$, $p<0.001$). The lower the cell's CV$_{\text{base}}$, the higher the SI($r$) value and thus the stronger the triangular structure of its second-order susceptibility. The model cells with the most distinct triangular pattern in their second-order susceptibility are the ones with the lowest CVs, hinting at low intrinsic noise levels.
+The SI($r$) correlates with the CVs of the cell's baseline interspike intervals ($r=-0.60$, $p<0.001$). The lower the cell's CV$_{\text{base}}$, the higher the SI($r$) value and thus the stronger the triangular structure of its second-order susceptibility. The model cells with the most distinct triangular pattern in their second-order susceptibility are the ones with the lowest CVs, hinting at low intrinsic noise levels.
 
 
 \subsection{Weakly nonlinear interactions vanish for higher stimulus contrasts}
-As pointed out above, the weakly nonlinear regime can only be observed for sufficiently weak stimulus amplitudes. In the model cells we estimated second-order susceptibilities for RAM stimuli with a contrast of 1, 3, and 10\,\%. The estimates for 1\,\% contrast (\figrefb{fig:modelsusceptcontrasts}\,\panel[ii]{E}) were quite similar to the estimates from the noise-split method, corresponding to a stimulus contrast of 0\,\% ($r=0.97$, $p\ll 0.001$). Thus, RAM stimuli with 1\,\% contrast are sufficiently small to not destroy weakly nonlinear interactions by their linearizing effect. At this low contrast, 51\,\% of the model cells have an SI($r$) value greater than 1.2.
+As pointed out above, the weakly nonlinear regime can only be observed for sufficiently weak stimuli. In the model cells we estimated second-order susceptibilities for RAM stimuli with a contrast of 1, 3, and 10\,\%. The estimates for 1\,\% contrast (\figrefb{fig:modelsusceptcontrasts}\,\panel[ii]{E}) were quite similar to the estimates from the noise-split method, corresponding to a stimulus contrast of 0\,\% ($r=0.97$, $p\ll 0.001$). Thus, RAM stimuli with 1\,\% contrast are sufficiently small to not destroy weakly nonlinear interactions by their linearizing effect. At this low contrast, 51\,\% of the model cells have an SI($r$) value greater than 1.2.
 
 At a RAM contrast of 3\,\% the SI($r$) values become smaller (\figrefb{fig:modelsusceptcontrasts}\,\panel[iii]{E}). Only 7 cells (18\,\%) have SI($r$) values exceeding 1.2. Finally, at 10\,\% the SI($r$) values of all cells drop below 1.2, except for three cells (8\,\%, \figrefb{fig:modelsusceptcontrasts}\,\panel[iv]{E}). The cell shown in \subfigrefb{fig:modelsusceptcontrasts}{A} is one of them. At 10\,\% contrast the SI($r$) values are no longer correlated with the ones in the noise-split configuration ($r=0.32$, $p=0.05$). To summarize, the regime of distinct nonlinear interactions at frequencies matching the baseline firing rate extends in this set of P-unit model cells to stimulus contrasts ranging from a few percents to about 10\,\%.
 
@@ -382,9 +384,12 @@ At a RAM contrast of 3\,\% the SI($r$) values become smaller (\figrefb{fig:model
 \end{figure}
 
 \subsection{Weakly nonlinear interactions can be deduced from limited data}
-Estimating second-order susceptibilities reliably requires large numbers (millions) of FFT segments (\figrefb{fig:noisesplit}). Electrophysiological measurements, however, suffer from limited recording durations and hence limited numbers of available FFT segments and estimating weakly nonlinear interactions from just a few hundred segments appears futile. The question arises, to what extend such limited-data estimates are still informative?
+Estimating second-order susceptibilities reliably requires large numbers (millions) of FFT segments (\figrefb{fig:noisesplit}). Electrophysiological measurements, however, suffer from limited recording durations and %hence limited numbers of available FFT segments and
+estimating weakly nonlinear interactions from just a few hundred segments appears futile. To what extend are such limited-data estimates still informative?
 
-The second-order susceptibility matrices that are based on only 100 segments look flat and noisy, lacking the triangular structure (\subfigref{fig:modelsusceptlown}{B}). The anti-diagonal ridge, however, where the sum of the stimulus frequencies matches the neuron's baseline firing rate, seems to be present whenever the converged estimate shows a clear triangular structure (compare \subfigref{fig:modelsusceptlown}{B} and \subfigref{fig:modelsusceptlown}{A}). The SI($r$) characterizes the height of the ridge in the second-oder susceptibility plane at the neuron's baseline firing rate $r$. Comparing SI($r$) values based on 100 FFT segments to the ones based on one or ten million segments for all 39 model cells (\subfigrefb{fig:modelsusceptlown}{C}) supports this impression. They correlate quite well at contrasts of 1\,\% and 3\,\% ($r=0.9$, $p\ll 0.001$). At a contrast of 10\,\% this correlation is weaker ($r=0.38$, $p<0.05$), because there are only three cells left with SI($r$) values greater than 1.2. Despite the good correlations, care has to be taken to set a threshold on the SI($r$) values for deciding whether a triangular structure would emerge for a much higher number of segments. Because at low number of segments the estimates are noisier, there could be false positives for a too low threshold. Setting the threshold to 1.8 avoids false positives for the price of a few false negatives.
+The second-order susceptibility matrices that are based on only 100 segments look flat and noisy, lacking the triangular structure (\subfigref{fig:modelsusceptlown}{B}). The anti-diagonal ridge, however, %where the sum of the stimulus frequencies matches the neuron's baseline firing rate,
+seems to be present whenever the converged estimate shows a clear triangular structure (compare \subfigref{fig:modelsusceptlown}{B} and \subfigref{fig:modelsusceptlown}{A}). %The SI($r$) characterizes the height of the ridge in the second-oder susceptibility plane at the neuron's baseline firing rate $r$.
+Comparing SI($r$) values based on 100 FFT segments to the ones based on one or ten million segments for all 39 model cells (\subfigrefb{fig:modelsusceptlown}{C}) supports this impression. They correlate quite well at contrasts of 1\,\% and 3\,\% ($r=0.9$, $p\ll 0.001$). At a contrast of 10\,\% this correlation is weaker ($r=0.38$, $p<0.05$), because there are only three cells left with SI($r$) values greater than 1.2. Despite the good correlations, care has to be taken to set a threshold on the SI($r$) values for deciding whether a triangular structure would emerge for a much higher number of segments. Because at low number of segments the estimates are noisier, there could be false positives for a too low threshold. Setting the threshold to 1.8 avoids false positives for the price of a few false negatives.
 
 Overall, observing SI($r$) values greater than about 1.8, even for a number of FFT segments as low as one hundred, seems to be a reliable indication for a triangular structure in the second-order susceptibility at the corresponding stimulus contrast. Small stimulus contrasts of 1\,\% are less informative, because of their bad signal-to-noise ratio. Intermediate stimulus contrasts around 3\,\% seem to be optimal, because there, most cells still have a triangular structure in their susceptibility and the signal-to-noise ratio is better. At RAM stimulus contrasts of 10\,\% or higher the signal-to-noise ratio is even better, but only few cells remain with weak triangularly shaped susceptibilities that might be missed as a false positives.
 
@@ -433,7 +438,7 @@ Overall, observing SI($r$) values greater than about 1.8, even for a number of F
 \end{figure*}
 
 \subsection{Low CVs and weak responses predict weakly nonlinear responses}
-Now we are prepared to evaluate our pool of 39 P-unit model cells, 172 P-units, and 30 ampullary afferents recorded in 80 specimen of \textit{Apteronotus leptorhynchus}. For comparison across cells we condensed the structure of the second-order susceptibilities into SI($r$) values, \eqnref{siindex}. In order to make the data comparable, both model and experimental SI($r$) estimates, \eqnref{siindex}, are based on 100 FFT segments.
+Now we are prepared to evaluate our pool of 39 P-unit model cells, 172 P-units, and 30 ampullary afferents recorded in 80 specimen of \textit{Apteronotus leptorhynchus}. For direct comparison across cells we condensed the structure of the second-order susceptibilities into SI($r$) values, \eqnref{siindex}. Both, model and experimental SI($r$) estimates, \eqnref{siindex}, are based on 100 FFT segments.
 
 In the P-unit models, each model cell contributed with three RAM stimulus presentations with contrasts of 1, 3, and 10\,\%, resulting in $n=117$ data points. 19 (16\,\%) had SI($r$) values larger than 1.8, indicating the expected ridges at the baseline firing rate in their second-order susceptibility. The lower the cell's baseline CV, i.e. the less intrinsic noise, the higher the SI($r$) (\figrefb{fig:dataoverview}\,\panel[i]{A}).
 
@@ -441,11 +446,11 @@ The effective stimulus strength also plays a role in predicting the SI($r$) valu
 
 In comparison to the experimentally measured P-unit recordings, the model cells are skewed to lower baseline CVs (Mann-Whitney $U=13986$, $p=3\times 10^{-9}$), because the models are not able to reproduce bursting, which we observe in many P-units and which leads to high CVs. Also the response modulation of the models is skewed to lower values (Mann-Whitney $U=15312$, $p=7\times 10^{-7}$), because in the measured cells, response modulation is positively correlated with baseline CV (Pearson $R=0.45$, $p=1\times 10^{-19}$), i.e. bursting cells are more sensitive. Median baseline firing rate in the models is by 53\,Hz smaller than in the experimental data (Mann-Whitney $U=17034$, $p=0.0002$).
 
-In the experimentally measured P-units, each of the $172$ units contributes on average with two RAM stimulus presentations, presented at contrasts ranging from 1 to 20\,\% to the 376 samples. Despite the mentioned differences between the P-unit models and the measured data, the SI($r$) values do not differ between models and data (median of 1.3, Mann-Whitney $U=19702$, $p=0.09$) and also 16\,\% of the samples from all presented stimulus contrasts exceed the threshold of 1.8. The SI($r$) values of the P-unit population correlate weakly with the CV of the baseline ISIs that range from 0.18 to 1.35 (median 0.49). Cells with lower baseline CVs tend to have more pronounced ridges in their second-order susceptibilities than those with higher baseline CVs (\figrefb{fig:dataoverview}\,\panel[i]{B}).
+In the experimentally measured P-units, each of the $172$ cells contributes on average with two RAM stimulus presentations, presented at contrasts ranging from 1 to 20\,\% to the 376 samples. Despite the mentioned differences between the P-unit models and the measured data, the SI($r$) values do not differ between models and data (median of 1.3, Mann-Whitney $U=19702$, $p=0.09$) and also 16\,\% of the samples from all presented stimulus contrasts exceed the threshold of 1.8. The SI($r$) values of the P-unit population correlate weakly with the CV of the baseline ISIs that range from 0.18 to 1.35 (median 0.49). Cells with lower baseline CVs tend to have more pronounced ridges in their second-order susceptibilities than those with higher baseline CVs (\figrefb{fig:dataoverview}\,\panel[i]{B}).
 
-Samples with weak responses to a stimulus, be it an insensitive P-unit or a weak stimulus, have higher SI($r$) values and thus a more pronounced ridge in the second-order susceptibility in comparison to strongly responding cells, most of them having flat second-order susceptibilities (\figrefb{fig:dataoverview}\,\panel[ii]{B}). P-units with low or high baseline firing rates can have large SI($r$) (\figrefb{fig:dataoverview}\,\panel[iii]{B}). How pronounced nonlinear response components are in P-units thus depends on the baseline CV (a proxy for the internal noise level), and the response strength during stimulation (effective output noise). 
+Samples with weak responses to a stimulus, due to low sensitivity or a weak stimulus, have higher SI($r$) values in comparison to strongly responding cells, most of them having flat second-order susceptibilities (\figrefb{fig:dataoverview}\,\panel[ii]{B}). P-units with low or high baseline firing rates can have large SI($r$) (\figrefb{fig:dataoverview}\,\panel[iii]{B}). How pronounced nonlinear response components are in P-units thus depends on the baseline CV (a proxy for the internal noise level), and the response strength during stimulation (effective output noise). 
 
-The population of ampullary cells is generally more homogeneous, with lower baseline CVs than P-units (Mann-Whitney $U=33464$, $p=9\times 10^{-49}$). Accordingly, SI($r$) values of ampullary cells (median 2.3) are indeed higher than in P-units (median 1.3, Mann-Whitney $U=6450$, $p=2\times 10^{-19}$). 52 samples (58\,\%) with SI($r$) values greater than 1.8 would have a triangular structure in their second-order susceptibilities. Ampullary cells also show a negative correlation with baseline CV, despite their narrow distribution of CVs ranging from 0.03 to 0.15 (median 0.09) (\figrefb{fig:dataoverview}\,\panel[i]{C}). Again, sensitive cells with stronger response modulations are at the bottom of the SI($r$) distribution with values close to one (\figrefb{fig:dataoverview}\,\panel[ii]{C}). As in P-units, baseline firing rate does not predict SI($r$) values (\figrefb{fig:dataoverview}\,\panel[iii]{C}).
+The population of ampullary cells is more homogeneous, with generally lower baseline CVs than P-units (Mann-Whitney $U=33464$, $p=9\times 10^{-49}$). Accordingly, SI($r$) values of ampullary cells (median 2.3) are indeed higher than in P-units (median 1.3, Mann-Whitney $U=6450$, $p=2\times 10^{-19}$). 52 samples (58\,\%) with SI($r$) values greater than 1.8 would have a triangular structure in their second-order susceptibilities. Ampullary cells also show a negative correlation with baseline CV, despite their narrow distribution of CVs ranging from 0.03 to 0.15 (median 0.09) (\figrefb{fig:dataoverview}\,\panel[i]{C}). Again, sensitive cells with stronger response modulations are at the bottom of the SI($r$) distribution with values close to one (\figrefb{fig:dataoverview}\,\panel[ii]{C}). Similar to P-units, the baseline firing rate does not predict SI($r$) values (\figrefb{fig:dataoverview}\,\panel[iii]{C}).
 
 \begin{figure*}[t]
   \includegraphics[width=\columnwidth]{model_full.pdf}
@@ -453,19 +458,19 @@ The population of ampullary cells is generally more homogeneous, with lower base
 \end{figure*}
 
 \subsection{Second-order susceptibility explains nonlinear responses to pure sinewave stimuli}
-Using the RAM stimulation we found pronounced nonlinear responses, in particular in the limit to weak stimulus amplitudes. How do these findings relate to the situation of pure sine wave stimuli with finite amplitudes that approximate the interference of EODs of real animals? For the P-units, the relevant signals are the beat frequencies $f_1$ and $f_2$ that arise from the interference of either of the two foreign EODs with the receiving fish's own EOD (\figref{fig:motivation}). In the introductory example, the response power spectrum showed peaks from nonlinear interactions at the sum of the two beat frequencies (orange marker, \subfigrefb{fig:motivation}{D}) and at the difference between the two beat frequencies (red marker, \subfigrefb{fig:motivation}{D}). In this example, $f_2$ was similar to the baseline firing rate, corresponding to the horizontal ridge of the second-order susceptibility (\subfigrefb{fig:lifresponse}{B}). The difference frequency is not covered by the so-far shown part of the second-order susceptibility, in which only the response at the sum of the two stimulus frequencies is addressed.
+Using the RAM stimulation we found pronounced nonlinear responses, in particular in the limit to weak stimulus amplitudes. How do these findings relate to the situation of pure sine wave stimuli \notejg{with finite amplitudes, needed?} that approximate the interference of EODs of real animals? For the P-units, the relevant signals are the beat frequencies $f_1$ and $f_2$ that arise from the interference of either of the two foreign EODs with the receiving fish's own EOD (\figref{fig:motivation}). In the introductory example, the response power spectrum showed peaks from nonlinear interactions at the sum of the two beat frequencies (orange marker, \subfigrefb{fig:motivation}{D}) and at the difference between the two beat frequencies (red marker, \subfigrefb{fig:motivation}{D}). In this example, $f_2$ was similar to the baseline firing rate, corresponding to the horizontal ridge of the second-order susceptibility (\subfigrefb{fig:lifresponse}{B}). The difference frequency is not covered by the so-far shown part of the second-order susceptibility, in which only the response at the sum of the two stimulus frequencies is addressed.
 
 However, the second-order susceptibility \eqnref{chi2} is a spectral measure that is based on Fourier transforms and thus is also defined for negative stimulus frequencies. The full $\chi_2$ matrix is symmetric with respect to the origin. In the upper-right and lower-left quadrants of $|\chi_2|(f_1, f_2)$, stimulus-evoked responses at the sum $f_1 + f_2$ are shown, whereas in the lower-right and upper-left quadrants responses at the difference $f_2 - f_1$ are shown (\figref{fig:model_full}). The vertical and horizontal lines at $f_1 = r$ and $f_2 = r$ are very pronounced in the upper-right quadrant of \subfigrefb{fig:model_full}{A} and extend into the upper-left quadrant, where they fade out towards more negative $f_1$ frequencies. The peak in the response power-spectrum at $f_2 - f_1$ 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}, \citealp{Schlungbaum2023}). 
 
 Is it possible to predict nonlinear responses in a three-fish setting based on second-order susceptibility matrices estimated from RAM stimulation in the noise-split configuration (\subfigrefb{fig:model_full}{A})? We test this by stimulating the same model with two weak beats (\subfigrefb{fig:model_full}{B--E}). Qualitatively, the second-order susceptibility predicts the presence and absence of peaks in the response spectrum at the sums and differences of the two beat frequencies. However, a quantitative prediction fails. This is because pure sine waves influence a nonlinear system in different ways than a white-noise stimulus. This will be addressed in detail in a future manuscript.
-
+\notejg{This section is not unimportant but opens many problems, e.g. the quantitative mismatch. Do we need it? Could move essential parts to the discussion, i.e. the f1-f2. have the full spectrum in figure 1 and ``gloss'' over the quality of the prediction.}
 
 \section{Discussion}
 
-Theoretical work \citep{Voronenko2017,Franzen2023} derived analytical expressions for weakly-nonlinear responses in LIF and theta model neurons driven by two sine waves with distinct frequencies. We here investigated such nonlinear responses in two types of electroreceptor afferents that differ in their intrinsic noise levels \citep{Grewe2017} using band-limited white-noise stimuli to estimate second-order susceptibilities. Following \citet{Voronenko2017} we expected to observe distinct ridges in the second-order susceptibility where either of the stimulus frequencies alone or their sum matche the baseline firing rate. We find traces of these nonlinear responses in the majority of the low-noise ampullary afferents and in only less than a fifth of the P-units that are characterized by low intrinsic noise levels and low output noise. Complementary model simulations demonstrate in the limit of high numbers of FFT segments, that the second order susceptibilities estimated from the electrophysiological data  are indeed indicative of the theoretically expected triangular structure of super.threshold weakly nonlinear responses. With this, we provide experimental evidence for nonlinear responses of a spike generator at low intrinsic noise levels or low stimulus amplitudes.
+Theoretical work \citep{Voronenko2017,Franzen2023} derived analytical expressions for weakly-nonlinear responses in LIF and theta model neurons driven by two sine waves with distinct frequencies. We here investigated such nonlinear responses in two types of electroreceptor afferents that differ in their intrinsic noise levels \citep{Grewe2017} using band-limited white-noise stimuli to estimate second-order susceptibilities. Following \citet{Voronenko2017} we expected to observe distinct ridges in the second-order susceptibility where either of the stimulus frequencies alone or their sum matches the baseline firing rate. We find traces of these nonlinear responses in the majority of ampullary afferents. In P-units, however, only a minority of the recorded cells, i.e. those characterized by low intrinsic noise levels and low output noise, show signs of such nonlinear responses. Complementary model simulations demonstrate in the limit of high numbers of FFT segments, that the estimates from the electrophysiological data are indeed indicative of the theoretically expected triangular structure of supra-threshold weakly nonlinear responses. With this, we provide experimental evidence for nonlinear responses of a spike generator at low intrinsic noise levels or low stimulus amplitudes. \notejg{not sure, if i like this last sentence}
 
 \subsection{Intrinsic noise limits nonlinear responses}
-The weakly nonlinear regime with its triangular pattern of elevated second-order susceptibility resides between the linear and a stochastic mode-locking regime . Too strong intrinsic noise linearizes the system and wipes out the triangular structure of the second-order susceptibility \citep{Voronenko2017}. Our electrophysiological recordings match this theoretical expectation. The lower the coefficient of variations of the P-units' baseline interspike intervals, the more cells show the expected nonlinearities (\subfigref{fig:dataoverview}{B}). Still, in only 18\,\% of the P-units analyzed in this study we find relevant nonlinear responses. On the other hand, the majority of the ampullary cells have generally lower CVs (median of 0.09) and indeed in the majority of ampullary afferents (73\,\%) we observe nonlinear responses (\subfigrefb{fig:dataoverview}{C}).
+The weakly nonlinear regime with its triangular pattern of elevated second-order susceptibility resides between the linear and a stochastic mode-locking regime. Too strong intrinsic noise linearizes the system and wipes out the triangular structure of the second-order susceptibility \citep{Voronenko2017}. Our electrophysiological recordings match this theoretical expectation. The lower the coefficient of variations of the P-units' baseline interspike intervals, the more cells show the expected nonlinearities (\subfigref{fig:dataoverview}{B})\notejg{There's something wrong here}. Still, only 18\,\% of the P-units analyzed in this study show relevant nonlinear responses. On the other hand, the majority of the ampullary cells have generally lower CVs (median of 0.09) and indeed we can observe nonlinear responses in the majority of ampullary afferents (74\,\%, \subfigrefb{fig:dataoverview}{C}).
 
 The CV is a proxy for the intrinsic noise in the cells (\citealp{Vilela2009}, note however the effect of coherence resonance for excitable systems close to a bifurcation, \citealp{Pikovsky1997, Lindner2004}). In both cell types, we observe a negative correlation between the second-order susceptibility at $f_1 + f_2 = r$ and the baseline CV (\figrefb{fig:dataoverview}), indicating that it is the level of intrinsic noise that shapes nonlinear responses. These findings are in line with previous theoretical and experimental studies showing the linearizing effects of noise in nonlinear systems \citep{Roddey2000, Chialvo1997, Voronenko2017}. Increased intrinsic noise has been demonstrated to increase the CV and to 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}. Only in cells with sufficiently low levels of intrinsic noise, weakly nonlinear responses can be observed.
 
@@ -482,14 +487,15 @@ Making assumptions about the nonlinearities in a system also reduces the amount
 \subsection{Nonlinear encoding in ampullary cells}
 The afferents of the passive electrosensory system, the ampullary cells, exhibit strong second-order susceptibilities (\figref{fig:dataoverview}). Ampullary cells more or less directly translate external low-frequency electric fields into afferent spikes, much like in the standard LIF and theta models used by Lindner and colleagues \citep{Voronenko2017,Franzen2023}. Indeed, we observe in the ampullary cells similar second-order nonlinearities as in LIF models.
 
-Ampullary stimuli originate from the muscle potentials induced by prey movement \citep{Kalmijn1974, Engelmann2010, Neiman2011fish}. For a single prey items such as \textit{Daphnia}, these potentials are often periodic but the simultaneous activity of a swarm of prey resembles Gaussian white noise \citep{Neiman2011fish}. Linear and nonlinear encoding in ampullary cells has been studied in great detail in the paddlefish \citep{Neiman2011fish}. The power spectrum of the baseline response shows two main peaks: One peak at the baseline firing frequency, a second one at the oscillation frequency of primary receptor cells in the epithelium, and interactions of both. Linear encoding in the paddlefish shows a gap at the epithelial oscillation frequency, instead, nonlinear responses are very pronounced there.
+Ampullary stimuli originate from the muscle potentials induced by prey movement \citep{Kalmijn1974, Engelmann2010, Neiman2011fish}. For a single prey items such as \textit{Daphnia}, these potentials are often periodic but the simultaneous activity of a swarm of prey resembles Gaussian white noise \citep{Neiman2011fish}. Linear and nonlinear encoding in ampullary cells has been studied in great detail in the paddlefish \citep{Neiman2011fish}. The power spectrum of the baseline response shows two main peaks: One peak at the baseline firing frequency, a second one at the oscillation frequency of primary receptor cells in the epithelium, plus interactions of both. Linear encoding in the paddlefish shows a gap at the epithelial oscillation frequency, instead, nonlinear responses are very pronounced there.
 
 Ampullary stimulus encoding is somewhat different in \lepto{}. The power spectrum of the spontaneous response is dominated by only the baseline firing rate and its harmonics, a second oscillator is not visible. The baseline firing 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 ampullary cells increases with stimulus intensity while it disappears in our case (\subfigrefb{fig:dataoverview}{F}) indicating that paddlefish data have been recorded above the weakly-nonlinear regime. 
 
 The population of ampullary cells is very homogeneous with respect to the baseline rate (131$\pm$29\,Hz) and stimulus encoding properties \citep{Grewe2017}. This implies that, if the stimulus contains the appropriate frequency components that sum up to the baseline rate, the resulting nonlinear response appears at the baseline rate that is similar in the full population of ampullary cells and that is outside the linear coding range. Postsynaptic cells integrating ampullary input might be able to extract such nonlinear responses. How such nonlinear effects might influence prey detection should be addressed in future studies.
 
 \subsection{Nonlinear encoding in P-units}
-Noise stimuli have the advantage that a range of frequencies can be measured with a single stimulus presentation and they have been successfully applied to characterize sensory coding in many systems \citep{French1973, Marmarelis1999, Borst1999, Chacron2005, Grewe2017}. In the real-world situation of wave-type electric fish, however, the natural stimuli encoded by P-units are periodic amplitude modulations of the self-generated electric field which arise from the superposition of the own and foreign EODs. Natural interactions usually occur between low numbers of close-by fish and thus the AMs are a mixture of a few distinct frequencies with substantial amplitudes \citep{Stamper2010,Fotowat2013, Henninger2020}. How informative are the second-order susceptibilities observed under noise stimulation for the encoding of distinct frequencies?
+\notejg{This section is a bit wild. Parts are repetitions of things that have been said earlier already. Other parts are nice but there is no real flow. We could/should merge is with the last section of the results section.}
+Noise stimuli have the advantage that a range of frequencies can be measured with a single stimulus presentation and they have been successfully applied to characterize sensory coding in many systems \citep{French1973, Marmarelis1999, Borst1999, Chacron2005, Grewe2017}. The natural stimuli encoded by P-units are, however, periodic amplitude modulations of the self-generated electric field which arise from the superposition of the own and foreign EODs. Such interactions usually occur between low numbers of close-by fish and thus the AMs are a mixture of a few distinct frequencies with substantial amplitudes \citep{Stamper2010,Fotowat2013, Henninger2020}. How informative are the second-order susceptibilities observed under noise stimulation for the encoding of distinct frequencies?
 
 We here applied broad-band noise amplitude-modulation stimuli to estimate the second-order susceptibility. Broadband noise stimuli introduce additional noise that linearizes the dynamics of the system. In contrast, pure sine wave stimulation is spectrally focused and drives the system on the background of the intrinsic noise. This explains why we can observe nonlinear interactions between sine wave stimuli with distinct frequencies (\figrefb{fig:motivation}) although they are barely visible when stimulating with strong broad-band noise (\figrefb{fig:noisesplit}).
 
@@ -501,7 +507,7 @@ The weakly nonlinear interactions in low-CV P-units could facilitate the detecta
 
 \subsection{Conclusions}
 
-We have demonstrated pronounced nonlinear responses in primary electrosensory afferents at weak stimulus amplitudes and sufficiently low intrinsic noise levels. The observed nonlinearities match the expectations from previous theoretical studies \citep{Voronenko2017,Franzen2023}. The resulting nonlinear components introduce spectral components not present in the original stimulus, but may provide an edge in the context of signal detection problems at stimulus amplitudes close to threshold \citep{Schlungbaum2023}. Electrosensory afferents share an evolutionary history with hair cells \citep{Baker2019} and share many response properties with mammalian auditory nerve fibers \citep{Barayeu2023, Joris2004}. Thus, we expect weakly nonlinear responses for near-threshold stimulation in auditory nerve fibers as well.
+We have demonstrated pronounced nonlinear responses in primary electrosensory afferents at weak stimulus amplitudes and sufficiently low intrinsic noise levels. The observed nonlinearities match the expectations from previous theoretical studies \citep{Voronenko2017,Franzen2023}. The resulting nonlinear components introduce spectral components not present in the original stimulus, that may provide an edge in the context of signal detection problems at stimulus amplitudes close to threshold \citep{Schlungbaum2023}. Electrosensory afferents share an evolutionary history with hair cells \citep{Baker2019} and share many response properties with mammalian auditory nerve fibers \citep{Barayeu2023, Joris2004}. Thus, we expect weakly nonlinear responses for near-threshold stimulation in auditory nerve fibers as well.
 
 
 \section{Acknowledgements}