updating model_full description and figure

model_and_data.py

@@ -186,7 +186,7 @@ def model_and_data2(eod_metrice=False, width=0.005, nffts=['whole'], powers=[1],
 
         axes = [ax_ams[0], axes[1], axes[2], ax_ams[1], axes[3], axes[4], ]
 
-        fig.tag([ax_data], xoffs=-3, yoffs=1.6)
+        fig.tag(ax_data, xoffs=-3, yoffs=1.6)
         fig.tag([axes[0:3]], xoffs=-3, yoffs=1.6)
         fig.tag([axes[3:6]], xoffs=-3, yoffs=1.6)
 

susceptibility1.tex

@@ -435,33 +435,40 @@ Field observations have shown that courting males were able to react to distant
 
 \section{Results}
 
+\begin{figure*}[t]
+  \includegraphics[width=\columnwidth]{motivation.pdf}
+  \caption{\label{fig:motivation} Nonlinearity in an electrophysiologically recorded P-unit of \lepto{} in a three-fish setting (cell identifier ``2021-08-03-ac"). Receiver with EOD frequency $\feod{} =664$\,Hz encounters fish with EOD frequencies $f_{1}=631$\,Hz and $f_{2}=797$\,Hz. Both encountered fish lead to a beat contrast of 10\,\%. Top: Scheme of a nonlinear system. Second row: Interference of the receiver EOD with the EODs of other fish. Third row: Spike trains of the P-unit. Fourth row: Firing rate, retrieved as the convolution of the spike trains with a Gaussian kernel ($\sigma = 1$\,ms). Bottom row: Power spectrum of the firing rate. \figitem{A} Baseline condition: Only the receiver is present. The baseline firing rate \fbase{} dominates the power spectrum of the firing rate. \figitem{B} The receiver and the fish with EOD frequency $f_{1}=631$\,Hz are present. \figitem{C} The receiver and the fish with EOD frequency $f_{2}=797$\,Hz are present. \figitem{D} All three fish with the EOD frequencies \feod{}, $f_{1}$ and $f_{2}$ are present. Nonlinear peaks occur at the sum and difference of the two beat frequencies in the power spectrum of the firing rate.
+  }
+\end{figure*}
+
 Theoretical work on leaky integrate-and-fire and conductance-based models suggests a distinct structure of the second-order response function for neurons with low levels of intrinsic noise driven in the super-threshold regime with low stimulus amplitudes (\figrefb{fig:lifresponse}, \citealp{Voronenko2017}). Here, we re-analyze 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} together with simulations of LIF-based models of P-unit spiking to search for such weakly nonlinear responses in real neurons. We start with a few example P-units to demonstrate the basic concepts.
 
 
 \subsection{Nonlinear responses in P-units stimulated with two beat frequencies}
 
+
+
 \begin{figure*}[t]
-  \includegraphics[width=\columnwidth]{motivation.pdf}
-  \caption{\label{fig:motivation} Nonlinearity in an electrophysiologically recorded P-unit of \lepto{} in a three-fish setting (cell identifier ``2021-08-03-ac"). Receiver with EOD frequency $\feod{} =664$\,Hz encounters fish with EOD frequencies $f_{1}=631$\,Hz and $f_{2}=797$\,Hz. Both encountered fish lead to a beat contrast of 10\,\%. Top: Scheme of a nonlinear system. Second row: Interference of the receiver EOD with the EODs of other fish. Third row: Spike trains of the P-unit. Fourth row: Firing rate, retrieved as the convolution of the spike trains with a Gaussian kernel ($\sigma = 1$\,ms). Bottom row: Power spectrum of the firing rate. \figitem{A} Baseline condition: Only the receiver is present. The baseline firing rate \fbase{} dominates the power spectrum of the firing rate. \figitem{B} The receiver and the fish with EOD frequency $f_{1}=631$\,Hz are present. \figitem{C} The receiver and the fish with EOD frequency $f_{2}=797$\,Hz are present. \figitem{D} All three fish with the EOD frequencies \feod{}, $f_{1}$ and $f_{2}$ are present. Nonlinear peaks occur at the sum and difference of the two beat frequencies in the power spectrum of the firing rate.
-  }
+  \includegraphics[width=\columnwidth]{nonlin_regime.pdf}
+  \caption{\label{fig:nonlin_regime} Nonlinear response of a model P-unit in response to increasing beat amplitudes in a three-fish or two-beat setting. The model used has the cell identifier 2013-01-08-aa (table~\ref{modelparams}). \figitem{A--D}. Top -- stimulus consisting of two beats. The beat frequencies are 30\,Hz (\bone{}) and 130\,Hz (\btwo{}). \btwo{} is equal to the baseline firing rate \fbase{}. The beat contrasts are equal in the respective example, and increase from \panel{A} to \panel{D}. Middle -- spike response of the model P-unit to the stimulus above. Bottom -- power spectrum of the firing rate of this P-unit. Nonlinear effects at \bsum{} (orange marker) and at \bdiff{} (red marker) increase for intermediate contrasts (\panel{B}), decrease for stronger contrasts (\panel{C}) and again increase for very strong contrasts (\panel{D}). \figitem{E} Amplitude of the linear (\bone{} and \btwo{}) and nonlinear (\bdiff{} and \bsum{}) responses of the model P-units plotted for increasing beat contrasts (contrasts equal for both beats). }
 \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 \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 chosen 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 the 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.
 
+
 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. 
 
+
 The response of a P-unit to varying beat amplitudes has been estimated by leaky-integrate-and-fire (LIF) models, fitted to the baseline firing properties of electrophysiologically measured P-units. In the chosen P-units model nonlinear peaks (red and orange markers) appear for intermediate beat stimuli (\subfigrefb{fig:motivation}{B}), decreases for stronger stimuli (\subfigrefb{fig:motivation}{C}) and again emerges for very strong stimuli (\subfigrefb{fig:motivation}{D}). Thus two regimes of nonlinearity can be observed for intermediate and strong beat amplitudes (\subfigrefb{fig:motivation}{E}).
 
 For this example, we have chosen two specific stimulus (beat) frequencies. For a full characterization of the nonlinear responses, we need to measure the response of the P-units to many different combinations of stimulus frequencies. 
 
-\begin{figure*}[t]
-  \includegraphics[width=\columnwidth]{nonlin_regime.pdf}
-  \caption{\label{fig:nonlin_regime} Nonlinear response of a model P-unit in response to increasing beat amplitudes in a three-fish or two-beat setting. The model used has the cell identifier 2013-01-08-aa (table~\ref{modelparams}). \figitem{A--D}. Top -- stimulus consisting of two beats. The beat frequencies are 30\,Hz (\bone{}) and 130\,Hz (\btwo{}). \btwo{} is equal to the baseline firing rate \fbase{}. The beat contrasts are equal in the respective example, and increase from \panel{A} to \panel{D}. Middle -- spike response of the model P-unit to the stimulus above. Bottom -- power spectrum of the firing rate of this P-unit. Nonlinear effects at \bsum{} (orange marker) and at \bdiff{} (red marker) increase for intermediate contrasts (\panel{B}), decrease for stronger contrasts (\panel{C}) and again increase for very strong contrasts (\panel{D}). \figitem{E} Amplitude of the linear (\bone{} and \btwo{}) and nonlinear (\bdiff{} and \bsum{}) responses of the model P-units plotted for increasing beat contrasts (contrasts equal for both beats). }
-\end{figure*}
+
+
 
 
 \subsection{Nonlinear signal transmission in low-CV P-units}
@@ -492,9 +499,8 @@ In contrast, a high-CV P-unit (CV$_{\text{base}}=0.4$) does not exhibit pronounc
 Irrespective of the CV, neither of the two example P-units shows the complete expected structure of nonlinear interactions. Electric fish possess an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogenous 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 its harmonics. 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 loses its distinct peak at \fsum{} and its overall 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-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}. 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-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}. 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.pdf}
@@ -505,11 +511,11 @@ In the example recordings shown above (\figsrefb{fig:cells_suscept} and \fref{fi
 
 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.
+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 mainly missing.
 
 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 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 \, 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 used in the experiments reveals key features of the nonlinear response.
+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^6$ (\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 used in the experiments reveals key features of the nonlinear response.
 
 With high levels of intrinsic noise, we would not expect the nonlinear response features to survive. Indeed, we do not find these features in a high-CV P-unit and its corresponding model (not shown).
 
@@ -522,8 +528,8 @@ We estimated the second-order susceptibility of P-unit responses using RAM stimu
   \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 2013-01-08-aa, 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. The baseline firing rate of this cell was at $\fbase=120$\,Hz. The position of the orange/red letters corresponds to the beat frequencies used for the stimulation with pure sine waves in the subsequent panels and indicates the sum/difference of those beat frequencies.  \figitem{B--E} 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 contrast of beat beats is 0.02.  Colored circles highlight the height of selected peaks in the power spectrum. 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{}.}
 \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}, \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}).
+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}, \citealp{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}{B--E}). 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}{B}), 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}{C}). If only one beat frequency is equal to \fbase{}, both a peak at the sum and at the difference frequency are present in the P-unit response (\subfigrefb{fig:model_full}{D}). 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}{E}).
 
 \begin{figure*}[tp]
   \includegraphics[width=\columnwidth]{data_overview_mod.pdf}

susceptibility1.tex.bak

@@ -435,26 +435,18 @@ Field observations have shown that courting males were able to react to distant
 
 \section{Results}
 
-Theoretical work on leaky integrate-and-fire and conductance-based models suggests a distinct structure of the second-order response function for neurons with low levels of intrinsic noise driven in the super-threshold regime with low stimulus amplitudes (\figrefb{fig:lifresponse}, \citealp{Voronenko2017}). Here, we re-analyze 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} together with simulations of LIF-based models of P-unit spiking to search for such weakly nonlinear responses in real neurons. We start with a few example P-units to demonstrate the basic concepts.
-
-
-\subsection{Nonlinear responses in P-units stimulated with two beat frequencies}
-
 \begin{figure*}[t]
   \includegraphics[width=\columnwidth]{motivation.pdf}
   \caption{\label{fig:motivation} Nonlinearity in an electrophysiologically recorded P-unit of \lepto{} in a three-fish setting (cell identifier ``2021-08-03-ac"). Receiver with EOD frequency $\feod{} =664$\,Hz encounters fish with EOD frequencies $f_{1}=631$\,Hz and $f_{2}=797$\,Hz. Both encountered fish lead to a beat contrast of 10\,\%. Top: Scheme of a nonlinear system. Second row: Interference of the receiver EOD with the EODs of other fish. Third row: Spike trains of the P-unit. Fourth row: Firing rate, retrieved as the convolution of the spike trains with a Gaussian kernel ($\sigma = 1$\,ms). Bottom row: Power spectrum of the firing rate. \figitem{A} Baseline condition: Only the receiver is present. The baseline firing rate \fbase{} dominates the power spectrum of the firing rate. \figitem{B} The receiver and the fish with EOD frequency $f_{1}=631$\,Hz are present. \figitem{C} The receiver and the fish with EOD frequency $f_{2}=797$\,Hz are present. \figitem{D} All three fish with the EOD frequencies \feod{}, $f_{1}$ and $f_{2}$ are present. Nonlinear peaks occur at the sum and difference of the two beat frequencies in the power spectrum of the firing rate.
   }
 \end{figure*}
 
+Theoretical work on leaky integrate-and-fire and conductance-based models suggests a distinct structure of the second-order response function for neurons with low levels of intrinsic noise driven in the super-threshold regime with low stimulus amplitudes (\figrefb{fig:lifresponse}, \citealp{Voronenko2017}). Here, we re-analyze 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} together with simulations of LIF-based models of P-unit spiking to search for such weakly nonlinear responses in real neurons. We start with a few example P-units to demonstrate the basic concepts.
 
 
-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 chosen 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 the 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.
+\subsection{Nonlinear responses in P-units stimulated with two beat frequencies}
 
-For this example, we have chosen two specific stimulus (beat) 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. 
 
-For a full characterization of the nonlinear responses, we need to measure the response of the P-units to many different combinations of stimulus frequencies. Such a P-unit response can be estimated by leaky-integrate-and-fire (LIF) models, fitted to the baseline firing properties of electrophysiologically measured P-units. In the chosen P-units model a nonlinear peak appears for intermediate beat stimuli (\subfigrefb{fig:motivation}{B}) decreases for stronger stimuli (\subfigrefb{fig:motivation}{C}) and again emerges for very strong stimuli (\subfigrefb{fig:motivation}{D}). 
 
 \begin{figure*}[t]
   \includegraphics[width=\columnwidth]{nonlin_regime.pdf}
@@ -462,6 +454,23 @@ For a full characterization of the nonlinear responses, we need to measure the r
 \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 \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 chosen 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 the 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.
+
+
+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. 
+
+
+The response of a P-unit to varying beat amplitudes has been estimated by leaky-integrate-and-fire (LIF) models, fitted to the baseline firing properties of electrophysiologically measured P-units. In the chosen P-units model nonlinear peaks (red and orange markers) appear for intermediate beat stimuli (\subfigrefb{fig:motivation}{B}), decreases for stronger stimuli (\subfigrefb{fig:motivation}{C}) and again emerges for very strong stimuli (\subfigrefb{fig:motivation}{D}). Thus two regimes of nonlinearity can be observed for intermediate and strong beat amplitudes (\subfigrefb{fig:motivation}{E}).
+
+For this example, we have chosen two specific stimulus (beat) frequencies. For a full characterization of the nonlinear responses, we need to measure the response of the P-units to many different combinations of stimulus frequencies. 
+
+
+
+
+
 \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 \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}).
 
@@ -490,9 +499,8 @@ In contrast, a high-CV P-unit (CV$_{\text{base}}=0.4$) does not exhibit pronounc
 Irrespective of the CV, neither of the two example P-units shows the complete expected structure of nonlinear interactions. Electric fish possess an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogenous 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 its harmonics. 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 loses its distinct peak at \fsum{} and its overall 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-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}. 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-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}. 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.pdf}
@@ -503,11 +511,11 @@ In the example recordings shown above (\figsrefb{fig:cells_suscept} and \fref{fi
 
 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.
+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 mainly missing.
 
 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 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 \, 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 used in the experiments reveals key features of the nonlinear response.
+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^6$ (\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 used in the experiments reveals key features of the nonlinear response.
 
 With high levels of intrinsic noise, we would not expect the nonlinear response features to survive. Indeed, we do not find these features in a high-CV P-unit and its corresponding model (not shown).