worked on Benjamins comments

.gitignore

@@ -1,4 +1,4 @@
-__pychache__
+__pycache__/
 *.blg
 *.aux
 *.bbl

model_full.py

@@ -1,4 +1,3 @@
-
 import sys
 import time
 
@@ -572,4 +571,4 @@ def plt_data_full_model(c1, chose_score, detections, devs, dfs, end, grid, mult_
 
 if __name__ == '__main__':
 
-    model_full()
+    model_full()

susceptibility1.tex

@@ -423,7 +423,7 @@ The transfer function used to describe linear properties of a system is the firs
 
 Noise in nonlinear systems, however, linearizes the system's response properties \citep{Yu1989, Chialvo1997}. Vice versa, in the limit to small stimuli, nonlinear systems can be well described by linear response theory \citep{Roddey2000, Doiron2004, Rocha2007, Sharafi2013}. With increasing stimulus amplitude, the contribution of the second-order kernel of the Volterra series becomes more relevant. For these weakly nonlinear responses analytical expressions for the second-order susceptibility have been derived for leaky-integrate-and-fire (LIF) \citep{Voronenko2017} and theta model neurons \citep{Franzen2023}. In the suprathreshold regime, where the LIF generates a baseline firing rate in the absence of an external stimulus, the linear response function has a peak at the baseline firing rate and its harmonics (\subfigrefb{fig:lifresponse}{A}) and the second-order susceptibility shows very distinct ridges of elevated nonlinear responses, exactly where two stimulus frequencies equal or add up to the neuron's baseline firing rate (\subfigrefb{fig:lifresponse}{B}). In experimental data, such structures in the second-order susceptibility have not been reported yet.
 
-Here we search for such weakly nonlinear responses in electroreceptors of the two electrosensory systems of the wave-type electric fish \textit{Apteronotus leptorhynchus}, i.e. the tuberous (active) and the ampullary (passive) electrosensory system. The p-type electroreceptors of the active system (P-units) are driven by the fish's high-frequency, quasi-sinusoidal electric organ discharges (EOD) and encode disturbances of it \citep{Bastian1981a}. The electroreceptors of the passive system are tuned to lower-frequency exogeneous electric fields such as caused by muscle activity of prey \citep{Kalmijn1974}. As different animals have different EOD-frequencies, being exposed to stimuli of multiple distinct frequencies is part of the animal's everyday life \citep{Benda2020,Henninger2020} and weakly non-linear interactions may occur in the electrosensory periphery. In communication contexts \citep{Walz2014, Henninger2018} the EODs of interacting fish superimpose and lead to periodic amplitude modulations (AMs or beats) of the receiver's EOD. Non-linear mechanisms in P-units, enable encoding of AMs in their time-dependent firing rates \citep{Bastian1981a, Walz2014, Middleton2006, Barayeu2023}. When multiple animals interact, the EOD interferences induce second-order amplitude modulations referred to as envelopes \citep{Yu2005, Fotowat2013, Stamper2012Envelope} and saturation non-linearities allow also for the encoding of these in the electrosensory periphery \citep{Savard2011}. Field observations have shown that courting males were able to react to the extremely weak signals of distant intruding males despite the strong foreground EOD of the nearby female \citep{Henninger2018}. Weakly nonlinear interactions at particular combinations of signals can be of immediate relevance in such settings as they could boost detectability of the faint signals \citep{Schlungbaum2023}.
+Here we search for such weakly nonlinear responses in electroreceptors of the two electrosensory systems of the wave-type electric fish \textit{Apteronotus leptorhynchus}, i.e. the tuberous (active) and the ampullary (passive) electrosensory system. The p-type electroreceptors of the active system (P-units) are driven by the fish's high-frequency, quasi-sinusoidal electric organ discharges (EOD) and encode disturbances of it \citep{Bastian1981a}. The electroreceptors of the passive system are tuned to lower-frequency exogeneous electric fields such as caused by muscle activity of prey \citep{Kalmijn1974}. As different animals have different EOD-frequencies, being exposed to stimuli of multiple distinct frequencies is part of the animal's everyday life \citep{Benda2020,Henninger2020} and weakly nonlinear interactions may occur in the electrosensory periphery. In communication contexts \citep{Walz2014, Henninger2018} the EODs of interacting fish superimpose and lead to periodic amplitude modulations (AMs or beats) of the receiver's EOD. Nonlinear mechanisms in P-units, enable encoding of AMs in their time-dependent firing rates \citep{Bastian1981a, Walz2014, Middleton2006, Barayeu2023}. When multiple animals interact, the EOD interferences induce second-order amplitude modulations referred to as envelopes \citep{Yu2005, Fotowat2013, Stamper2012Envelope} and saturation nonlinearities allow also for the encoding of these in the electrosensory periphery \citep{Savard2011}. Field observations have shown that courting males were able to react to the extremely weak signals of distant intruding males despite the strong foreground EOD of the nearby female \citep{Henninger2018}. Weakly nonlinear interactions at particular combinations of signals can be of immediate relevance in such settings as they could boost detectability of the faint signals \citep{Schlungbaum2023}.
 
 %The population of P-units is heterogeneous with respect to their baseline firing rates as well as in their intrinsic noise levels, as quantified by the coefficient of variation (CV) of the interspike intervals (ISI) \citep{Grewe2017, Hladnik2023}. Low-CV P-units have a less noisy firing pattern that is closer to pacemaker firing, whereas high-CV P-units show a more irregular firing pattern that is closer to a Poisson process. The population of ampullary cells of the passive electrosensory system, on the other hand, is homogeneous in their response properties and CVs are low \citep{Grewe2017}.
 
@@ -449,19 +449,18 @@ When stimulating with both foreign signals simultaneously, additional peaks appe
 \subsection{Linear and weakly nonlinear regimes}
 \begin{figure*}[tp]
   \includegraphics[width=\columnwidth]{nonlin_regime.pdf}
-  \caption{\label{fig:nonlin_regime} Linear and nonlinear responses of a model P-unit in a three-fish setting in dependence on increasing stimulus amplitudes. The model P-unit was stimulated with two sinewaves of equal amplitude (contrast) at difference frequencies $\bone=30$\,Hz and $\btwo=130$\,Hz relative the receiver's EOD frequency. \btwo{} was set to match the baseline firing rate \fbase{} of the P-unit. The model used has the cell identifier 2013-01-08-aa (table~\ref{modelparams}). \figitem{A--D} Top: the stimulus, an amplitude modulation of the receiver's EOD resulting from the stimulation with the two sine waves. The contrasts of both beats increase from \panel{A} to \panel{D} as indicated. Middle: Spike raster of the model P-unit response. Bottom: power spectrum of the firing rate estimated from the spike raster with a Gaussian kernel ($\sigma=1$\,ms). \figitem{A} At very low stimulus contrasts the response is linear. The only peaks in the response spectrum are at the two stimulating beat frequencies (green and purple marker). \figitem{B} At higher stimulus contrasts, nonlinear responses appear at the sum and the difference of the stimulus frequencies (orange and red marker). \figitem{C} At even stronger stimulation additional peaks appear in the power spectrum. \figitem{D} At a contrast of 10\,\% the response at the sum of the stimulus frequencies almost disappears. \figitem{E} Amplitude of the linear (\bone{}, \btwo{}) and nonlinear (\bdiff{}, \bsum{}) responses of the model P-unit as a function of beat contrast (thick lines). Thin lines indicate the initial linear and quadratic dependence on stimulus amplitude for the linear and non-linear responses, respectively.}
+  \caption{\label{fig:nonlin_regime} Linear and nonlinear responses of a model P-unit in a three-fish setting in dependence on increasing stimulus amplitudes. The model P-unit was stimulated with two sinewaves of equal amplitude (contrast) at difference frequencies $\bone=30$\,Hz and $\btwo=130$\,Hz relative the receiver's EOD frequency. \btwo{} was set to match the baseline firing rate \fbase{} of the P-unit. The model used has the cell identifier 2013-01-08-aa (table~\ref{modelparams}). \figitem{A--D} Top: the stimulus, an amplitude modulation of the receiver's EOD resulting from the stimulation with the two sine waves. The contrasts of both beats increase from \panel{A} to \panel{D} as indicated. Middle: Spike raster of the model P-unit response. Bottom: power spectrum of the firing rate estimated from the spike raster with a Gaussian kernel ($\sigma=1$\,ms). \figitem{A} At very low stimulus contrasts the response is linear. The only peaks in the response spectrum are at the two stimulating beat frequencies (green and purple marker). \figitem{B} At higher stimulus contrasts, nonlinear responses appear at the sum and the difference of the stimulus frequencies (orange and red marker). \figitem{C} At even stronger stimulation additional peaks appear in the power spectrum. \figitem{D} At a contrast of 10\,\% the response at the sum of the stimulus frequencies almost disappears. \figitem{E} Amplitude of the linear (\bone{}, \btwo{}) and nonlinear (\bdiff{}, \bsum{}) responses of the model P-unit as a function of beat contrast (thick lines). Thin lines indicate the initial linear and quadratic dependence on stimulus amplitude for the linear and nonlinear responses, respectively.}
 \end{figure*}
 
 The stimuli used in \figref{fig:motivation} had the same not-small amplitude. Whether this stimulus conditions falls into the weakly nonlinear regime as in \citet{Voronenko2017} is not clear. In order to illustrate how the responses to two beat frequencies develop over a range of amplitudes we use a stochastic leaky-integrate-and-fire (LIF) based P-unit model fitted to a specific electrophysiologically measured cell \citep{Barayeu2023}.
 
-At very low stimulus contrasts (less than approximately 0.5\,\% relative to the receiver's EOD amplitude) the spectrum has small peaks only at the beat frequencies (note, $\Delta f_2 = f_{base}$, \subfigref{fig:nonlin_regime}{A}, green and purple). The amplitudes of these peaks initially increase linearly with stimulus amplitude (\subfigref{fig:nonlin_regime}{E}, thin lines), an indication of the linear response at lowest stimulus amplitudes.
+At very low stimulus contrasts (in the example cell less than approximately 0.5\,\% relative to the receiver's EOD amplitude) the spectrum has small peaks only at the beat frequencies (note, $\Delta f_2 = f_{base}$, \subfigref{fig:nonlin_regime}{A}, green and purple). The amplitudes of these peaks initially increase linearly with stimulus amplitude (\subfigref{fig:nonlin_regime}{E}, thin lines), an indication of the linear response at lowest stimulus amplitudes.
 
-The linear regime is followed by the weakly non-linear regime. In addition to peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:nonlin_regime}{B}, orange and red). The amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:nonlin_regime}{E}, thin lines).
+This linear regime is followed by the weakly nonlinear regime. In addition to peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:nonlin_regime}{B}, orange and red). The amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:nonlin_regime}{E}, thin lines).
 
 At higher stimulus amplitudes (\subfigref{fig:nonlin_regime}{C \& D}) additional peaks appear in the response spectrum. The linear response and the weakly-nonlinear response start to deviate from their linear and quadratic dependence on amplitude (\subfigref{fig:nonlin_regime}{E}). The responses may even decrease for intermediate stimulus contrasts (\subfigref{fig:nonlin_regime}{D}). At high stimulus contrasts, additional nonlinearities in the system, in particular clipping of the firing rate, shape the responses.
 
-
+For this example, we have chosen two specific stimulus (beat) frequencies. One of these matching the P-unit's baseline firing rate. In the following, however, we are interested in how the nonlinear responses depend on different combinations of stimulus frequencies in the weakly nonlinear regime.
-For this example, we have chosen two specific stimulus (beat) frequencies. One of these matching the P-unit's baseline firing rate. In the following, however, we are interested in how the non-linear responses depend on different combinations of stimulus frequencies in the weakly nonlinear regime.
 
 
 \subsection{Nonlinear signal transmission in low-CV P-units}
@@ -475,7 +474,7 @@ Weakly nonlinear responses are expected in cells with sufficiently low intrinsic
 
 Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:punit}{C}, top trace, red line), are commonly used to characterize stimulus-driven responses of sensory neurons using transfer functions (first-order susceptibility), spike-triggered averages, or stimulus-response coherences. Here, we additionally estimate the second-order susceptibility to quantify nonlinear encoding. P-unit spikes align more or less clearly with fluctuations in the RAM stimulus. A higher stimulus intensity, here a higher contrast of the RAM relative to the EOD amplitude (see methods), entrains the P-unit response more clearly (light and dark purple for low and high contrast stimuli, \subfigrefb{fig:punit}{C}). Linear encoding, quantified by the transfer function, \eqnref{linearencoding_methods}, is similar for the two RAM contrasts in this low-CV P-unit (\subfigrefb{fig:punit}{D}), as expected for a linear system. The first-order susceptibility is low for low frequencies, peaks in the range below 100\,Hz and then falls off again \citep{Benda2005}.
 
-The second-order susceptibility, \eqnref{eq:susceptibility},  quantifies for each combination of two stimulus frequencies \fone{} and \ftwo{} the amplitude and phase of the stimulus-evoked response at the sum \fsum{} \notebl{but the  response spectrum is also multiplied with the stimulus spectra!}\notebl{Also add that this includes responses at f1-f2}. Large values of the second-order susceptibility indicate stimulus-evoked peaks in the response spectrum at the summed frequency that cannot be explained by linear response theory. Similar to the first-order susceptibility, the second-order susceptibility can be estimated directly from the response evoked by a RAM stimulus that stimulates the neuron with a whole range of frequencies simultaneously (\subfigsref{fig:punit}{E, F}). For LIF and theta neuron models driven in the supra-threshold regime, theory predicts nonlinear interactions between the two stimulus frequencies, when the two frequencies \fone{} and \ftwo{} or their sum \fsum{} exactly match the neuron's baseline firing rate \fbase{} \citep{Voronenko2017,Franzen2023}. Only then, additional stimulus-evoked peaks appear in the spectrum of the spiking response that would show up in the second-order susceptibility as a horizontal, a vertical, and an anti-diagonal line (\subfigrefb{fig:lifresponse}{B}, pink triangle in \subfigsref{fig:punit}{E, F}).
+The second-order susceptibility, \eqnref{eq:susceptibility},  quantifies for each combination of two stimulus frequencies \fone{} and \ftwo{} the amplitude and phase of the stimulus-evoked response at the sum \fsum{} (and also the difference, see below \subfigrefb{fig:model_full}{A}). Large values of the second-order susceptibility indicate stimulus-evoked peaks in the response spectrum at the summed frequency that cannot be explained by linear response theory. Similar to the first-order susceptibility, the second-order susceptibility can be estimated directly from the response evoked by a RAM stimulus that stimulates the neuron with a whole range of frequencies simultaneously (\subfigsref{fig:punit}{E, F}). For LIF and theta neuron models driven in the supra-threshold regime, theory predicts nonlinear interactions between the two stimulus frequencies, when the two frequencies \fone{} and \ftwo{} or their sum \fsum{} exactly match the neuron's baseline firing rate \fbase{} \citep{Voronenko2017,Franzen2023}. Only then, additional stimulus-evoked peaks appear in the spectrum of the spiking response that would show up in the second-order susceptibility as a horizontal, a vertical, and an anti-diagonal line (\subfigrefb{fig:lifresponse}{B}, pink triangle in \subfigsref{fig:punit}{E, F}).
 
 For the low-CV P-unit, we observe a band of slightly elevated second-order susceptibility for the low RAM contrast at \fsumb{} (yellowish anti-diagonal between pink edges, \subfigrefb{fig:punit}{E}). This structure vanishes for the stronger stimulus (\subfigref{fig:punit}{F}). Further, the overall level of the second-order susceptibility is reduced with increasing stimulus strength. To quantify the structural changes in the susceptibility matrices we projected the susceptibility values onto the diagonal by averaging over the anti-diagonals (\subfigrefb{fig:punit}{G}). At low RAM contrast this projected second-order susceptibility indeed has a small peak at \fbase{} (\subfigrefb{fig:punit}{G}, dot on top line). For the higher RAM contrast, however, this peak vanishes and the overall level of the second-order susceptibility is reduced (\subfigrefb{fig:punit}{G}). The reason behind this reduction is that a RAM with a higher contrast is not only a stimulus with an increased amplitude, but also increases the total noise in the system. Increased noise is known to linearize signal transmission \citep{Longtin1993, Chialvo1997,  Roddey2000, Voronenko2017} and thus the second-order susceptibility is expected to decrease.
 
@@ -555,12 +554,12 @@ Theoretical work \citep{Voronenko2017,Franzen2023} derived analytical expression
 % Weakly nonlinear responses versus saturation regime
 
 \subsection{Intrinsic noise limits nonlinear responses}
-The pattern of elevated  second-order susceptibility found in the experimental data matches the theoretical expectations only partially \notebl{too negative. The theory also predicts vanishing patterns for too much noise in the system}. Only P-units with low coefficients of variation (CV $<$ 0.25) of the interspike-interval distribution in their baseline response show the expected nonlinearities (\figref{fig:punit}, \figref{fig:model_full}, \subfigref{fig:dataoverview}{A}). Such low-CV cells are rare among the 221 P-units used in this study. On the other hand, the majority of the ampullary cells have generally lower CVs (median of 0.12) and have an approximately ten-fold higher level of second-order susceptibilities where \fsumb{} (\figref{fig:ampullary}, \subfigrefb{fig:dataoverview}{B}).
+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. Only P-units with low coefficients of variation (CV $<$ 0.25) of the interspike-interval distribution in their baseline response show the expected nonlinearities (\figref{fig:punit}, \figref{fig:model_full}, \subfigref{fig:dataoverview}{A}). Such low-CV cells are rare among the 221 P-units used in this study. On the other hand, the majority of the ampullary cells have generally lower CVs (median of 0.12) and have an approximately ten-fold higher level of second-order susceptibilities where \fsumb{} (\figref{fig:ampullary}, \subfigrefb{fig:dataoverview}{B}).
 
-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 \fsumb{} and the CV, 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.
+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 \fsumb{} and the 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.
 
 \subsection{Linearization by white-noise stimulation}
-Not only the intrinsic noise but also the stimulation with external white-noise linearizes the cells. This applies to both, the stimulation with AMs in P-units (\subfigrefb{fig:dataoverview}{E}) and direct stimulation in ampullary cells (\subfigrefb{fig:dataoverview}{F}). The stronger the effective stimulus, the less pronounced are the peaks in second-order susceptibility (see \subfigref{fig:punit}{E\&F} for a P-unit example and \subfigref{fig:ampullary}{E\&F} for an ampullary cell). This linearizing effect of noise stimuli limits the weakly nonlinear regime to small stimulus amplitudes. At higher stimulus amplitudes, however, other non-linearities of the system eventually show up in the second-order susceptibility.
+Not only the intrinsic noise but also the stimulation with external white-noise linearizes the cells. This applies to both, the stimulation with AMs in P-units (\subfigrefb{fig:dataoverview}{E}) and direct stimulation in ampullary cells (\subfigrefb{fig:dataoverview}{F}). The stronger the effective stimulus, the less pronounced are the peaks in second-order susceptibility (see \subfigref{fig:punit}{E\&F} for a P-unit example and \subfigref{fig:ampullary}{E\&F} for an ampullary cell). This linearizing effect of noise stimuli limits the weakly nonlinear regime to small stimulus amplitudes. At higher stimulus amplitudes, however, other nonlinearities of the system eventually show up in the second-order susceptibility.
 
 In order to characterize weakly nonlinear responses of the cells in the limit to vanishing stimulus amplitudes we utilized the Furutsu-Novikov theorem \citep{Novikov1965, Furutsu1963}. Following \citet{Lindner2022}, a substantial part of the intrinsic noise of a P-unit model \citep{Barayeu2023} is treated as signal. Performing this noise-split trick we can estimate the weakly nonlinear response without the linearizing effect of an additional external white noise stimulus. The model of a low-CV P-unit then shows the full nonlinear structure (\figref{model_and_data}) known from analytical derivations and simulations of basic LIF and theta models driven with pairs of sine-wave stimuli \citep{Voronenko2017,Franzen2023}.