fixed panel reference

Makefile

@@ -1,7 +1,7 @@
 TEXBASE=nonlinearbaseline
 BIBFILE=references.bib
-REBUTTALBASE=
-COVERBASE=cover2
+REBUTTALBASE=rebuttal3
+COVERBASE=cover3
 
 TEXFILE=$(TEXBASE).tex
 PDFFILE=$(TEXBASE).pdf
@@ -13,7 +13,7 @@ PT=$(wildcard *.py)
 PYTHONFILES=$(filter-out plotstyle.py spectral.py examplecells.py, $(PT))
 PYTHONPDFFILES=$(PYTHONFILES:.py=.pdf)
 
-REVISION=e3814a1be539f9424c17b7bd7ef45a8826a9f1e2
+REVISION=0a072e27a6736cc9384ed31efac16faabb0f3af3
 
 ifdef REBUTTALBASE
 REBUTTALTEXFILE=$(REBUTTALBASE).tex

data/2021-08-03-ac-invivo-1/info.dat

@@ -0,0 +1,58 @@
+# Recording:
+#     Recording quality: good
+#     Comment          : at end maldetection
+#     Experimenter     : Alexandra Rudnaya
+#     WaterTemperature : 26.5°C
+#     WaterConductivity: 300uS/cm
+#     Cell:
+#         CellType          : unkown
+#         Structure         : Nerve
+#         BrainRegion       : P-units
+#         BrainSubRegion    : ~
+#         Depth             : 1330um
+#         Lateral position  : 0mm
+#         Transverse section: 8
+#         Cell properties:
+#             Firing Rate1: 142Hz
+#             P-Value1    : 0.21
+#             X Position  : 0mm
+#             Y Position  : 0mm
+#     Subject:
+#         Species       : Apteronotus leptorhynchus
+#         Gender        : Female
+#         Size          : 15cm
+#         Weight        : 0.8g
+#         Identifier    : "2020lepto09"
+#         EOD Frequency : 667Hz
+#         Snout Position: 0mm
+#         Tail Position : 0mm
+#     Preparation:
+#         Type              : in vivo
+#         Anaesthesia       : true
+#         Anaesthetic       : MS 222
+#         AnaestheticDose   : 125mg/l
+#         LocalAnaesthesia  : true
+#         LocalAnaesthetic  : Lidocaine
+#         Immobilization    : true
+#         ImmobilizationDrug: Tubocurarin 5mg/L
+#         ImmobilizationDose: 75ul
+#     Name             : "2021-08-03-ac-invivo-1"
+#     Folder           : /home/efish/data/ephys/labrotation/2021-08-03-ac-invivo-1
+#     File             : [ trace-1.raw, trace-2.raw, trace-3.raw, trace-4.raw, trace-5.raw, stimulus-events.dat, restart-events.dat, recording-events.dat, spikes-1-events.dat, chirps-events.dat, localbeat-1-1-events.dat, localbeat-1-2-events.dat, stimuli.dat, stimulus-descriptions.dat ]
+#     Date             : "2021-08-03"
+#     Time             : "14:22:15.000"
+#     Recording duratio: 25.27min
+#     Mode             : Acquisition
+#     Software         : RELACS
+#     Software version : "0.9.8"
+# Setup:
+#     Identifier: invivo1
+#     Maintainer: Jan Grewe
+#     Creator   : Jan Grewe and Jan Benda
+#     Location  : "5A17"
+#     Lab       : Neuroethology Lab
+#     Institute : Institute for Neurobiology
+#     University: University Tuebingen
+#     Address   : Auf der Morgenstelle 28
+# Electrode:
+#     Resistance: 0MOhm

dataoverview.py

@@ -6,9 +6,6 @@ from scipy.stats import mannwhitneyu
 from thunderlab.tabledata import TableData
 from plotstyle import plot_style, lighter, significance_str
 
-
-data_path = Path('data')
-
 from noisesplit import model_cell as model_split_example
 from modelsusceptcontrasts import model_cells as model_contrast_examples
 from modelsusceptlown import model_cell as model_lown_example
@@ -18,6 +15,9 @@ from noisesplit import example_cell as punit_split_example
 from ampullaryexamplecell import example_cell as ampul_example
 from ampullaryexamplecell import example_cells as ampul_examples
 
+
+data_path = Path('data')
+
 model_examples = ([[model_lown_example, 0.01],
                    [model_lown_example, 0.03],
                    [model_lown_example, 0.1]],
@@ -26,6 +26,8 @@ model_examples = ([[model_lown_example, 0.01],
 punit_examples = (punit_example, [punit_split_example], punit_examples)
 ampul_examples = (ampul_example, [], ampul_examples)
 
+respmod = 'respmod1'
+
 
 def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
               si_thresh, example=[], split_example=[], examples=[]):
@@ -124,7 +126,7 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
     if 'cvbase' in xcol:
         ax.text(1, 0.64, f'$n={data.rows()}$', ha='right',
                 transform=ax.transAxes, fontsize='small')
-    print(f'    correlation {xcol:<8s} - {ycol}: r={r:5.2f}, p={p:.2e}')
+    print(f'    correlation {xcol:<11s} - {ycol:<11s}: r={r:5.2f}, p={p:.2e}')
     return cax
 
 
@@ -221,7 +223,7 @@ def plot_corr_contrast(ax, data, xcol, ycol, xpdfmax, color,
     if 'cvbase' in xcol:
         ax.text(1, 0.64, f'$n={data.rows()}$', ha='right',
                 transform=ax.transAxes, fontsize='small')
-    print(f'    correlation {xcol:<8s} - {ycol}: r={r:5.2f}, p={p:.2e}')
+    print(f'    correlation {xcol:<11s} - {ycol:<11s}: r={r:5.2f}, p={p:.2e}')
 
 
 def si_stats(title, data, sicol, si_thresh, nsegscol):
@@ -266,7 +268,7 @@ def si_stats(title, data, sicol, si_thresh, nsegscol):
     print(f'    total duration:     {np.min(duration):.1f}s - {np.max(duration):.1f}s, median={np.median(duration):.1f}s, mean={np.mean(duration):.1f}s, std={np.std(duration):.1f}s')
     duration = data['durationbase']
     print(f'    baseline duration:  {np.min(duration):.1f}s - {np.max(duration):.1f}s, median={np.median(duration):.1f}s, mean={np.mean(duration):.1f}s, std={np.std(duration):.1f}s')
-    cols = ['cvbase', 'respmod2', 'ratebase', 'vsbase', 'serialcorr1', 'burstfrac', 'ratestim', 'cvstim']
+    cols = ['cvbase', respmod, 'ratebase', 'vsbase', 'serialcorr1', 'burstfrac', 'ratestim', 'cvstim']
     for i in range(len(cols)):
         for j in range(i + 1, len(cols)):
             xcol = cols[i]
@@ -285,7 +287,7 @@ def plot_cvbase_si_punit(ax, data, ycol, si_thresh, color):
     ax.set_ylim(0, 6.5)
     ax.set_yticks_delta(2)
     examples = punit_examples if 'stimindex' in data else model_examples
-    #cax = plot_corr(ax, data, 'cvbase', ycol, 'respmod2', 0, 250, 3,
+    #cax = plot_corr(ax, data, 'cvbase', ycol, respmod, 0, 250, 3,
     #                'coolwarm', color, si_thresh, *examples)
     #cax.set_ylabel('Response mod.', 'Hz')
     plot_corr_contrast(ax, data, 'cvbase', ycol, 3, color,
@@ -300,7 +302,7 @@ def plot_cvstim_si_punit(ax, data, ycol, si_thresh, color):
     ax.set_ylim(0, 6.5)
     ax.set_yticks_delta(2)
     examples = punit_examples if 'stimindex' in data else model_examples
-    #cax = plot_corr(ax, data, 'cvstim', ycol, 'respmod2', 0, 250, 2,
+    #cax = plot_corr(ax, data, 'cvstim', ycol, respmod, 0, 250, 2,
     #                'coolwarm', color, si_thresh, *examples)
     #cax.set_ylabel('Response mod.', 'Hz')
     #cax = plot_corr(ax, data, 'cvstim', ycol, 'cvbase', 0, 1.5, 2,
@@ -318,16 +320,16 @@ def plot_cvstim_si_punit(ax, data, ycol, si_thresh, color):
     
 def plot_rmod_si_punit(ax, data, ycol, si_thresh, color):
     ax.set_xlabel('Response modulation', 'Hz')
-    ax.set_xlim(0, 250)
+    ax.set_xlim(0, 350)
     ax.set_xticks_delta(100)
     ax.set_ylabel('SI($r$)')
     ax.set_ylim(0, 6.5)
     ax.set_yticks_delta(2)
     examples = punit_examples if 'stimindex' in data else model_examples
-    #cax = plot_corr(ax, data, 'respmod2', ycol, 'cvbase', 0, 1.5, 0.016,
+    #cax = plot_corr(ax, data, respmod, ycol, 'cvbase', 0, 1.5, 0.016,
     #                'coolwarm', color, si_thresh, *examples)
     #cax.set_ylabel('CV$_{\\rm base}$')
-    plot_corr_contrast(ax, data, 'respmod2', ycol, 0.016, color,
+    plot_corr_contrast(ax, data, respmod, ycol, 0.016, color,
                        si_thresh, *examples)
 
     
@@ -355,7 +357,7 @@ def plot_cvbase_si_ampul(ax, data, ycol, si_thresh, color):
     ax.set_ylabel('SI($r$)')
     ax.set_ylim(0, 10)
     ax.set_yticks_delta(2)
-    #cax = plot_corr(ax, data, 'cvbase', ycol, 'respmod2', 0, 80, 20,
+    #cax = plot_corr(ax, data, 'cvbase', ycol, respmod, 0, 80, 20,
     #                'coolwarm', color, si_thresh, *ampul_examples)
     #cax.set_ylabel('Response mod.', 'Hz')
     plot_corr_contrast(ax, data, 'cvbase', ycol, 20, color,
@@ -369,7 +371,7 @@ def plot_cvstim_si_ampul(ax, data, ycol, si_thresh, color):
     ax.set_ylabel('SI($r$)')
     ax.set_ylim(0, 10)
     ax.set_yticks_delta(2)
-    #cax = plot_corr(ax, data, 'cvstim', ycol, 'respmod2', 0, 80, 6,
+    #cax = plot_corr(ax, data, 'cvstim', ycol, respmod, 0, 80, 6,
     #                'coolwarm', color, si_thresh, *ampul_examples)
     #cax.set_ylabel('Response mod.', 'Hz')
     #cax = plot_corr(ax, data, 'cvstim', ycol, 'cvbase', 0, 0.2, 6,
@@ -386,16 +388,16 @@ def plot_cvstim_si_ampul(ax, data, ycol, si_thresh, color):
     
 def plot_rmod_si_ampul(ax, data, ycol, si_thresh, color):
     ax.set_xlabel('Response modulation', 'Hz')
-    ax.set_xlim(0, 80)
-    ax.set_xticks_delta(20)
+    ax.set_xlim(0, 100)
+    ax.set_xticks_delta(30)
     ax.set_ylabel('SI($r$)')
     ax.set_ylim(0, 10)
     ax.set_yticks_delta(2)
-    #cax = plot_corr(ax, data, 'respmod2', ycol, 'cvbase', 0, 0.2, 0.06,
+    #cax = plot_corr(ax, data, respmod, ycol, 'cvbase', 0, 0.2, 0.06,
     #                'coolwarm', color, si_thresh, *ampul_examples)
     #cax.set_ylabel('CV$_{\\rm base}$')
     #cax.set_yticks_delta(0.1)
-    plot_corr_contrast(ax, data, 'respmod2', ycol, 0.06, color,
+    plot_corr_contrast(ax, data, respmod, ycol, 0.06, color,
                        si_thresh, *ampul_examples)
 
         
@@ -412,7 +414,7 @@ def plot_rate_si_ampul(ax, data, ycol, si_thresh, color):
     #cax.set_yticks_delta(0.1)
     plot_corr_contrast(ax, data, 'ratebase', ycol, 0.06, color,
                        si_thresh, *ampul_examples)
-    ax.legend(title='contrast', loc='upper right', bbox_to_anchor=(1.95, 1),
+    ax.legend(title='contrast', loc='upper right', bbox_to_anchor=(1.9, 1),
               markerfirst=False, handletextpad=0.5)
 
     
@@ -453,8 +455,8 @@ if __name__ == '__main__':
     print(f'    baseline rate model: min={np.min(ratemodel):3.0f}Hz max={np.max(ratemodel):3.0f}Hz median={np.median(ratemodel):3.0f}Hz')
     print(f'    baseline rate data:  min={np.min(ratedata):3.0f}Hz max={np.max(ratedata):3.0f}Hz median={np.median(ratedata):3.0f}Hz')
     print()
-    rmmodel = punit_model['respmod2']
-    rmdata = punit_data['respmod2']
+    rmmodel = punit_model[respmod]
+    rmdata = punit_data[respmod]
     u, p = mannwhitneyu(rmmodel, rmdata)
     print('Response modulation differs between P-unit models and data:')
     print(f'    U={u:g}, p={p:.2g}')
@@ -489,7 +491,7 @@ if __name__ == '__main__':
     s = plot_style()
     fig, axs = plt.subplots(3, 3, cmsize=(s.plot_width, 0.75*s.plot_width),
                             height_ratios=[1, 0, 1, 0.3, 1])
-    fig.subplots_adjust(leftm=6.5, rightm=18, topm=4.5, bottomm=4,
+    fig.subplots_adjust(leftm=6.5, rightm=17.5, topm=4.5, bottomm=4,
                         wspace=0.6, hspace=0.6)
     
     si_stats('P-unit model:', punit_model, 'dsinorm100', si_thresh,

modelsusceptlown.py

@@ -58,11 +58,12 @@ def plot_si_diags(ax, s, data, alphax, alphay, xthresh, ythresh, cell_name):
     l = linregress(six, siy)
     x = np.linspace(0, 10, 10)
     ax.set_visible(True)
+    ax.set_aspect('equal')
     ax.set_title(f'$c$={100*alphay:g}\\,\\%', fontsize='medium')
     ax.plot(x, x, **s.lsLine)
     ax.plot(x, l.slope*x + l.intercept, **s.lsGrid)
     ax.axhline(ythresh, **s.lsLine)
-    ax.axvline(xthresh, 0, 0.5, **s.lsLine)
+    ax.axvline(xthresh, 0, 0.65, **s.lsLine)
     if alphax == 0:
         mask = datax['triangle'] > 0.5
         ax.plot(six[mask], siy[mask], zorder=30, label='strong', **s.psA1m)
@@ -137,8 +138,8 @@ if __name__ == '__main__':
     xthresh = 1.2
     ythresh = 1.8
     s = plot_style()
-    fig, axs = plt.subplots(3, 4, cmsize=(s.plot_width, 0.6*s.plot_width),
-                            height_ratios=[3, 3, 0, 2])
+    fig, axs = plt.subplots(3, 4, cmsize=(s.plot_width, 0.62*s.plot_width),
+                            height_ratios=[3, 3, 0, 3])
     fig.subplots_adjust(leftm=7, rightm=9, topm=2, bottomm=4,
                         wspace=1, hspace=0.6)
     for ax in axs.flat:

nonlinearbaseline.tex

@@ -2,8 +2,6 @@
 
 \title{Spike generation in electroreceptor afferents introduces additional spectral response components by weakly nonlinear interactions}
 
-%\title{Weakly nonlinear interactions of the spike generator introduce additional spectral response components in electroreceptor afferents}
-
 \author{Alexandra Barayeu\textsuperscript{1},
   Maria Schlungbaum\textsuperscript{2,3},
   Benjamin Lindner\textsuperscript{2,3},
@@ -22,11 +20,14 @@
 %%%%% overall style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % number the lines:
-\newcommand{\setlinenumbers}{true}
+\newcommand{\setlinenumbers}{false}
 
 % use double spacing:
 \newcommand{\setdoublespacing}{false}
 
+% hide secrets (names and locations, for double-blind review):
+\newcommand{\hidesecrets}{false}
+
 % show authors:
 \newcommand{\showauthors}{true}
 
@@ -42,6 +43,11 @@
 % remove images from figures:
 \newcommand{\nofigs}{false}
 
+% for diff file:
+%\providecommand{\DIFaddtex}[1]{{\textbf{#1}}} %DIF PREAMBLE
+%\providecommand{\DIFdeltex}[1]{}              %DIF PREAMBLE
+
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%% line numbering, spacing, authors, keywords %%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -56,6 +62,8 @@
 \ifthenelse{\boolean{\setdoublespacing}}{%
 \usepackage{setspace}\doublespacing}{}
 
+\newcommand{\secret}[1]{\ifthenelse{\boolean{\hidesecrets}}{\colorbox{black}{\rule{2em}{1.5ex}}}{#1}}
+
 \ifthenelse{\boolean{\showauthors}}{}{\author{}\date{}}
 
 \usepackage[inline]{enumitem}
@@ -74,12 +82,12 @@
 \pagestyle{myheadings}
 \markright{\runningtitle}
 
-%%%%% language %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\usepackage[english]{babel}
-
 %%%%% fonts %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %\usepackage{pslatex}   % nice font for pdf file
 
+%%%%% language %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%\usepackage[english]{babel}
+
 %%%%% section style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \usepackage[sf,bf,it,big,clearempty]{titlesec}
 \usepackage{titling}
@@ -89,11 +97,8 @@
 
 %%%%% math %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \usepackage{amsmath}
-%\usepackage{lualatex-math}
 
 \usepackage{iftex}
-
-% LuaTeX-spezifischen Code nur ausführen, wenn LuaTeX verwendet wird
 \ifluatex
   \usepackage{lualatex-math}
 \fi
@@ -214,12 +219,6 @@
 \newcommand{\eqnref}[1]{Eq.~\eqref{#1}}
 \newcommand{\eqnsref}[1]{Eqs.~\eqref{#1}}
 
-%%%%% species names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\newcommand{\Lepto}{\emph{Apteronotus leptorhynchus}}
-\newcommand{\lepto}{\emph{A. leptorhynchus}}
-\newcommand{\Eigen}{\emph{Eigenmannia virescens}}
-\newcommand{\eigen}{\emph{E. virescens}}
-
 %%%%% notes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \newcommand{\note}[2][]{\textbf{[#1: #2]}}
 %\newcommand{\note}[2][]{}
@@ -236,27 +235,14 @@
 
 \maketitle
 
-\thispagestyle{empty}
+% 250 words
+\section{Abstract}
+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 one of 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} of either sex. In our combined experimental and modeling approach we identify these predicted nonlinear responses in those 31 out of 172 P-units that are characterized by low intrinsic noise. In contrast, the majority (22 out of 30) ampullary cells show nonlinear responses. 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.
 
-\noindent
-\begin{tabular}{@{}lr}
-  Number of figures & 10 \\
-  Number of tables & 0 \\
-  Number of multimedia & 0 \\
-  Number of words in abstract & 194 \\
-  Number of words in introduction & 659 \\
-  Number of words in discussion & 1891
-\end{tabular}
-
-\paragraph{Acknowledgements}
-We thank Tim Hladnik, Henriette Walz, Franziska Kuempfbeck, Fabian Sinz, Laura Seidler, Eva Vennemann, and Ibrahim Tunc for the data they recorded over the years in our lab.
-
-\paragraph{Conflict of interest:}The authors declare no conflict of interest.
-
-%\paragraph{Author Contributions:} All authors designed the study and discussed the results. AB performed the data analyses and modeling, AB and JG drafted the paper, all authors discussed and revised the manuscript.
-
-\paragraph{Founding sources}
-Supported by SPP 2205 ``Evolutionary optimisation of neuronal processing'' by the DFG, project number 430157666.
+% max 120 words
+\section{Significance statement}
+The generation of action potentials involves a strong threshold nonlinearity. Nevertheless, the encoding of stimuli with small amplitudes by neurons with sufficient intrinsic noise can be well described as a linear system. As the stimulus amplitude is increased, new spectral components start to appear in the so called weakly nonlinear regime. Theory predicts nonlinear interactions whenever one or the sum of two stimulus frequencies matches the neuron's baseline firing rate. Indeed, we find these interactions in a large set of electrophysiological recordings from primary electroreceptive afferents of a weakly electric fish. The nonlinear response components could boost sensory responses to weak signals emitted, for example, by distant conspecifics.
   
 \begin{keywords}
 \item Volterra series
@@ -266,75 +252,55 @@ Supported by SPP 2205 ``Evolutionary optimisation of neuronal processing'' by th
 \end{keywords}
 
 
-\newpage
-\setcounter{page}{1}
-
-\begin{center}
-  \sffamily\bfseries\LARGE Spike generation in electroreceptor afferents
-  introduces\\ additional spectral response components by\\ weakly
-  nonlinear interactions
-\end{center}
-
-% 250 words
-\section{Abstract}
-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.
-
-% max 120 words
-\section{Significance statement}
-The generation of action potentials involves a strong threshold nonlinearity. Nevertheless, the encoding of stimuli with small amplitudes by neurons with sufficient intrinsic noise can be well described as a linear system. As the stimulus amplitude is increased, new spectral components start to appear in the so called weakly nonlinear regime. Theory predicts nonlinear interactions whenever one or the sum of two stimulus frequencies matches the neuron's baseline firing rate. Indeed, we find these interactions in a large set of electrophysiological recordings from primary electroreceptive afferents of a weakly electric fish. The nonlinear response components could boost sensory responses to weak signals emitted, for example, by distant conspecifics.
-
-
 \section{Introduction}
 
-\begin{figure*}[t]
-  \includegraphics[width=\columnwidth]{lifsuscept}
-  \caption{\label{fig:lifsuscept} 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 (nonlinear) response function $|\chi_2(f_1, f_2)|$ quantifies the response at the sum of two stimulus frequencies. Because frequencies can also be negative, the sum is actually a difference in the upper-left and lower-right quadrants, since $f_2 + f_1 = f_2 - |f_1|$ for $f_1 < 0$. 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 \subfigref{fig:lifsuscept}{A}, 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.
 
 Noise in nonlinear systems, however, linearizes the system's response properties \citep{Yu1989, Chialvo1997}. Also, 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, in which 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:lifsuscept}{A}) and the second-order susceptibility shows very distinct ridges of elevated nonlinear responses, exactly where one of two stimulus frequencies equals or both frequencies add up to the neuron's baseline firing rate (\subfigrefb{fig:lifsuscept}{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 electroreceptor afferents 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}.
+\begin{figure*}[t]
+  \includegraphics[width=\columnwidth]{lifsuscept}
+  \caption{\label{fig:lifsuscept} 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. $r$ and $2r$ mark the neuron's baseline firing rate and its harmonic (dashed vertical lines). \figitem{B} Magnitude of the second-order (nonlinear) response function $|\chi_2(f_1, f_2)|$ quantifies the response at the sum of two stimulus frequencies. Because Fourier spectra have positive and negative frequencies, the sum is actually a difference in the upper-left and lower-right quadrants, since $f_2 + f_1 = f_2 - |f_1|$ for $f_1 < 0$. 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*}
+
+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 electroreceptor afferents 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 superpositions of the EODs 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}.
 
 
 \section{Materials and Methods}
 
 \subsection{Experimental subjects and procedures}
 
-Within this project, we re-evaluated datasets that were recorded between 2010 and 2023 at the Ludwig Maximilian University (LMU) M\"unchen and the Eberhard-Karls University T\"ubingen. All experimental protocols complied with national and European law and were approved by the respective Ethics Committees of the Ludwig-Maximilians Universität München (permit no. 55.2-1-54-2531-135-09) and the Eberhard-Karls Unversität Tübingen (permit no. ZP 1/13 and ZP 1/16).
-The final sample consisted of 172 P-units and 30 ampullary electroreceptor afferents recorded in 80 weakly electric fish of both sexes of the species \lepto{}. Fish were obtained from a commercial supplier for tropical fish (Aquarium Glaser GmbH, Rodgau, Germany) and kept in tanks with a water temperature of $25\pm1\,^\circ$C and a conductivity of around $270\,\micro\siemens\per\centi\meter$ under a 12\,h:12\,h light-dark cycle. 
+Within this project, we re-evaluated datasets that were recorded between 2010 and 2023 at the \secret{Ludwig Maximilian University (LMU) M\"unchen} and the \secret{Eberhard-Karls University T\"ubingen}. All experimental protocols complied with national and \secret{European} law and were approved by the respective Ethics Committees of the \secret{Ludwig-Maximilians Universität München} (permit no. \secret{55.2-1-54-2531-135-09}) and the \secret{Eberhard-Karls Unversität Tübingen} (permit no. \secret{ZP 1/13} and \secret{ZP 1/16}).
+The final sample consisted of 172 P-units and 30 ampullary electroreceptor afferents recorded in 80 weakly electric fish of both sexes of the species \textit{Apteronotus leptorhynchus}. Fish were obtained from a commercial supplier for tropical fish (\secret{Aquarium Glaser GmbH, Rodgau, Germany}) and kept in tanks with a water temperature of $25\pm1\,^\circ$C and a conductivity of around $270\,\micro\siemens\per\centi\meter$ under a 12\,h:12\,h light-dark cycle. 
 
 Before surgery, the animals were deeply anesthetized via bath application of a solution of MS222 (120\,mg/l, PharmaQ, Fordingbridge, UK) buffered with Sodium Bicarbonate (120\,mg/l). The posterior anterior lateral line nerve (pALLN) was exposed by making a small cut into the skin covering the nerve. The cut was placed dorsal of the operculum just before the nerve descends towards the anterior lateral line ganglion (ALLNG). Those parts of the skin that were to be cut were locally anesthetized by cutaneous application of liquid lidocaine hydrochloride (20\,mg/ml, bela-pharm GmbH). During the surgery, water supply was ensured by a mouthpiece to maintain anesthesia with a solution of MS222 (100\,mg/l) buffered with Sodium Bicarbonate (100\,mg/l). After surgery, fish were immobilized by intramuscular injection of from 25\,$\micro$l to 50\,$\micro$l of tubocurarine (5\,mg/ml dissolved in fish saline; Sigma-Aldrich).
 Respiration was then switched to normal tank water and the fish was transferred to the experimental tank.
 
 \subsection{Electrophysiological recordings} 
-For the recordings fish were positioned centrally in the experimental tank, with the major parts of their body submerged into the water. Those body parts that were above the water surface were covered with paper tissue to avoid drying of the skin. Local analgesia was refreshed in intervals of two hours by cutaneous application of Lidocaine (2\,\%; bela-pharm, Vechta, Germany) around the surgical wounds. Electrodes (borosilicate; 1.5\,mm outer diameter; GB150F-8P; Science Products, Hofheim, Germany) were pulled to a resistance of 50--100\,\mega\ohm{} (model P-97; Sutter Instrument, Novato, CA) and filled with 1\,M KCl solution. Electrodes were fixed in a microdrive (Luigs-Neumann, Ratingen, Germany) and lowered into the nerve. Recordings of electroreceptor afferents were amplified and lowpass filtered at 10\,kHz (SEC-05, npi-electronics, Tamm, Germany, operated in bridge mode). All signals, neuronal recordings, recorded EOD, and the generated stimulus, were digitized with sampling rates of 20 or 40\,kHz (PCI-6229, National Instruments, Austin, TX). RELACS (\url{https://github.com/relacs/relacs}) running on a Linux computer was used for online spike and EOD detection, stimulus generation, and calibration. Recorded data was then stored on the hard drive for offline analysis.
+For the recordings fish were positioned centrally in the experimental tank, with the major parts of their body submerged into the water. Those body parts that were above the water surface were covered with paper tissue to avoid drying of the skin. Local analgesia was refreshed in intervals of two hours by cutaneous application of Lidocaine (2\,\%; \secret{bela-pharm, Vechta, Germany}) around the surgical wounds. Electrodes (borosilicate; 1.5\,mm outer diameter; GB150F-8P; \secret{Science Products, Hofheim, Germany}) were pulled to a resistance of 50--100\,\mega\ohm{} (model P-97; Sutter Instrument, Novato, CA) and filled with 1\,M KCl solution. Electrodes were fixed in a microdrive (\secret{Luigs-Neumann, Ratingen, Germany}) and lowered into the nerve. Recordings of electroreceptor afferents were amplified and lowpass filtered at 10\,kHz (\secret{SEC-05, npi-electronics, Tamm, Germany}, operated in bridge mode). All signals, neuronal recordings, recorded EOD, and the generated stimulus, were digitized with sampling rates of 20 or 40\,kHz (PCI-6229, National Instruments, Austin, TX). RELACS (\url{https://github.com/relacs/relacs}) running on a Linux computer was used for online spike and EOD detection, stimulus generation, and calibration. Recorded data was then stored on the hard drive for offline analysis.
 
 \subsection{Identification of P-units and ampullary cells}
-Recordings were classified as P-units if baseline action potentials phase locked to the EOD with vectors strengths between 0.7 and 0.95, a baseline firing rate larger than 30\,Hz, a serial correlation of subsequent interspike intervals below zero, a coefficient of variation of baseline interspike intervals below 1.5 und during stimulation below 2. As ampullary cells we classified recordings with vector strengths below 0.15, baseline firing rate above 10\,Hz, baseline CV below 0.18, CV during stimulation below 1.0, and a response modulation during stimulation below 80\,Hz \citep{Grewe2017}. We here selected only those cells of which the neuron's baseline activity as well as the responses to band-limited white noise stimuli were recorded.
+Recordings were classified as P-units if baseline action potentials phase locked to the EOD with vectors strengths between 0.7 and 0.95, a baseline firing rate larger than 30\,Hz, a serial correlation of subsequent interspike intervals below zero, a coefficient of variation of baseline interspike intervals below 1.5 and during stimulation below 2. P-units are clearly distinguished from T-type electroreceptors, that we did not analyze here, by having firing rates much lower than the EOD frequency of the fish (no 1:1 locking to the EOD). As ampullary cells we classified recordings with vector strengths below 0.15, baseline firing rate above 10\,Hz, baseline CV below 0.18, CV during stimulation below 1.0, and a response modulation during stimulation below 80\,Hz \citep{Grewe2017}. We here selected only those cells of which the neuron's baseline activity as well as the responses to band-limited white noise stimuli were recorded.
 
 \subsection{Electric field recordings}
-For monitoring the EOD without the stimulus, two vertical carbon rods ($11\,\centi\meter$ long, 8\,mm diameter) in a head-tail configuration were placed isopotential to the stimulus. Their signal was differentially amplified with a gain factor between 100 and 500 (depending on the recorded animal) and band-pass filtered (3 to 1500\,Hz pass-band, DPA2-FX; npi electronics, Tamm, Germany). For an estimate of the transdermal potential that drives the electroreceptors, two silver wires spaced by 1\,cm were located next to the left gill of the fish and orthogonal to the fish's longitudinal body axis (amplification 100 to 500 times, band-pass filtered with 3 to 1\,500\,Hz pass-band, DPA2-FX; npi-electronics, Tamm, Germany). This local EOD measurement recorded the combination of the fish's own EOD and the applied stimulus.
+For monitoring the EOD without the stimulus, two vertical carbon rods ($11\,\centi\meter$ long, 8\,mm diameter) in a head-tail configuration were placed isopotential to the stimulus. Their signal was differentially amplified with a gain factor between 100 and 500 (depending on the recorded animal) and band-pass filtered (3 to 1500\,Hz pass-band, DPA2-FX; \secret{npi electronics, Tamm, Germany}). For an estimate of the transdermal potential that drives the electroreceptors, two silver wires spaced by 1\,cm were located next to the left gill of the fish and orthogonal to the fish's longitudinal body axis (amplification 100 to 500 times, band-pass filtered with 3 to 1\,500\,Hz pass-band, DPA2-FX; \secret{npi-electronics, Tamm, Germany}). This local EOD measurement recorded the combination of the fish's own EOD and the applied stimulus.
 
 
 \subsection{Stimulation}\label{rammethods}
-Electric stimuli were attenuated (ATN-01M, npi-electronics, Tamm, Germany), isolated from ground (ISO-02V, npi-electronics, Tamm, Germany) and delivered via two horizontal carbon rods (30 cm length, 8 mm diameter) located $15\,\centi\meter$ laterally to the fish.
-The fish were stimulated with band-limited white noise stimuli with a cut-off frequency of 150\,Hz (ampullary), 300 or 400\,Hz (P-units). The stimulus intensity is given as a contrast, i.e. the standard deviation of the noise stimulus relative to the fish's EOD amplitude. The contrast varied between 1 and 20\,\% (median 5\,\%) for P-units and 2.5 and 20\,\% (median 5\,\%) for ampullary cells. Only noise stimuli with a duration of at least 2\,s (maximum of 50\,s, median 10\,s) and enough repetitions to results in at least 100 FFT segments (see below, P-units: 100--1520, median 313, ampullary cells: 105 -- 3648, median 722) were included into the analysis. When ampullary cells were recorded, the white noise was directly applied as the stimulus. To create random amplitude modulations (RAM) for P-unit recordings, the noise stimulus was first multiplied with the EOD of the fish (MXS-01M; npi electronics).
+Electric stimuli were attenuated (ATN-01M, \secret{npi-electronics, Tamm, Germany}), isolated from ground (ISO-02V, \secret{npi-electronics, Tamm, Germany}) and delivered via two horizontal carbon rods (30 cm length, 8 mm diameter) located $15\,\centi\meter$ parallel to each side of the fish. The fish were stimulated with band-limited Gaussian white noise stimuli, i.e. signals with equal power at all frequencies up to a cut-off frequency and with a stationary Gaussian probability density. For the ampullary cells we chose a cut-off frequency of 150\,Hz, whereas for the P-units we used either 300 or 400\,Hz. The stimuli were generated by drawing normally distributed real and imaginary numbers for all frequencies up to the desired cut-off frequency in the Fourier domain and then applying an inverse Fourier transform \citep{Billah1990,Skorjanc2023}. The stimulus intensity is given as a contrast, i.e. the standard deviation of the resulting amplitude modulation relative to the fish's EOD amplitude. The contrast varied between 1 and 20\,\% (median 5\,\%) for P-units and 2.5 and 20\,\% (median 5\,\%) for ampullary cells. Only recordings with noise stimuli with a duration of at least 2\,s (maximum of 50\,s, median 10\,s) and enough repetitions to results in at least 100 FFT segments (see below, P-units: 100--1520, median 313, ampullary cells: 105 -- 3648, median 722) were included into the analysis. When ampullary cells were recorded, the noise stimuli $s(t)$ were directly applied as the stimulus and thus were simply added to the fish's own EOD: $s(t) + EOD(t)$. To create random amplitude modulations (RAM) for P-unit recordings, the noise stimulus was first multiplied with the EOD of the fish (MXS-01M; npi electronics) and then played back through the stimulation electrodes: $EOD(t) + s(t)EOD(t) = (1 + s(t))EOD(t)$.
 
-\subsection{Data analysis} Data analysis was done in Python 3 using the packages matplotlib \citep{Hunter2007}, numpy \citep{Walt2011}, scipy \citep{scipy2020}, pandas \citep{Mckinney2010}, nixio \citep{Stoewer2014}, and thunderlab (\url{https://github.com/bendalab/thunderlab}).
+\subsection{Data analysis} Data analysis was done in Python 3 using the packages matplotlib \citep{Hunter2007}, numpy \citep{Walt2011}, scipy \citep{scipy2020}, nixio \citep{Stoewer2014}, and thunderlab (\secret{\url{https://github.com/bendalab/thunderlab}}).
 
-\paragraph{Code accessibility}
-The P-unit model parameters and spectral analysis algorithms are available at \url{https://github.com/bendalab/punitmodel/tree/v1}.
+\paragraph{Code and data accessibility}
+The P-unit model parameters and spectral analysis algorithms are available at \secret{\url{https://github.com/bendalab/punitmodel/tree/v1}}. The electrophysiological data of all the P-units and Ampullary cells analyzed here, the analysis scripts, and the scripts simulating the P-unit model responses are accessible at \secret{\url{https://doi.org/10.12751/g-node.4w8ebq}}.
 
 \paragraph{Baseline analysis}\label{baselinemethods}
 The baseline firing rate $r$ was calculated as the number of spikes divided by the duration of the baseline recording (median 32\,s). The coefficient of variation (CV) of the interspike intervals (ISI) is their standard deviation relative to their mean: $\rm{CV}_{\rm base} = \sqrt{\langle (ISI- \langle ISI \rangle) ^2 \rangle} / \langle ISI \rangle$. If the baseline was recorded several times in a recording, the measures from the longest recording were taken.
 
 \paragraph{White noise analysis} \label{response_modulation}
-When stimulated with band-limited white noise stimuli, neuronal activity is modulated around the average firing rate that is similar to the baseline firing rate and in that way encodes the time-course of the stimulus. For an estimate of the time-dependent firing rate $r(t)$ we convolved each spike train with normalized Gaussian kernels with standard deviation of 2\,ms, if not mentioned otherwise, and averaged the resulting single-trail firing rates over trials. The response modulation quantifies the variation of $r(t)$ computed as the standard deviation in time $\sigma_{s} = \sqrt{\langle (r(t)-\langle r(t) \rangle_t )^2\rangle_t}$, where $\langle \cdot \rangle_t$ denotes averaging over time.
+When stimulated with band-limited white noise stimuli, neuronal activity is modulated around the average firing rate that is similar to the baseline firing rate and in that way encodes the time-course of the stimulus. For an estimate of the time-dependent firing rate $r(t)$ we convolved each spike train with normalized Gaussian kernels with standard deviation of 1\,ms and averaged the resulting single-trail firing rates over trials. The response modulation quantifies the variation of $r(t)$ computed as the standard deviation in time $\sigma_{s} = \sqrt{\langle (r(t)-\langle r(t) \rangle_t )^2\rangle_t}$, where $\langle \cdot \rangle_t$ denotes averaging over time.
 
 \paragraph{Spectral analysis}\label{susceptibility_methods}
 To characterize the relation between the spiking response evoked by white-noise  stimuli, we estimated the first- and second-order susceptibilities in the frequency domain. For this we converted spike times into binary vectors $x_k$ with $\Delta t = 0.5$\,ms wide bins that are set to 2\,kHz where a spike occurred and zero otherwise. Fast Fourier transforms (FFT) $S(\omega)$ and $X(\omega)$ of the stimulus $s_k$ (also down-sampled to a sampling rate of 2\,kHz) and $x_k$, respectively, were computed numerically according to
@@ -342,7 +308,7 @@ To characterize the relation between the spiking response evoked by white-noise
   \label{fourier}
   X(\omega) = \sum_{k=0}^{N-1} \, x_k e^{- i \omega k}
 \end{equation}
-for $N=512$ long segments of $T=N \Delta t = 256$\,ms duration with no overlap, resulting in a spectral resolution of about 4\,Hz. Note, that for a real Fourier integral a factor $\Delta t$ is missing. For simplicity we use angular frequencies $\omega=2\pi f$ instead of frequencies $f$.
+for $N=512$ long segments of $T=N \Delta t = 256$\,ms duration with no overlap, resulting in a spectral resolution of about 4\,Hz. In \figref{fig:regimes} we used $N=2^{18}$ and $\Delta t=0.1$\,ms (26.2\,s). Note, that for a real Fourier integral a factor $\Delta t$ is missing. For simplicity we use angular frequencies $\omega=2\pi f$ instead of frequencies $f$. Also note, that the Fourier transform of a signal results in both positive as well as negative frequencies. This is because the signal is decomposed into complex exponentials, not sine or cosine waves. To get back a real-valued sine or cosine wave, the components at positive and negative frequencies need to be combined. For example, the sum $e^{-\omega t} + e^{+\omega t}$ of two complex exponentials at frequencies $\pm\omega$ is $2\cos(\omega t)$.   
 
 In the experimental data the duration of the noise stimuli varied and they were presented once or repeatedly (frozen noise). For the analysis we discarded the responses within the initial 200\,ms of stimulation in each trial. To make the recordings comparable we always used the first 100 segments from as many trials as needed for the following analysis.
 
@@ -402,7 +368,7 @@ Values larger than one indicate a sharp ridge in the susceptibility matrix close
 
 
 \begin{figure*}[t]
-  \includegraphics[width=\columnwidth]{flowchart.pdf}
+  \includegraphics[width=\columnwidth]{flowchart}
   \caption{\label{flowchart}
     Architecture of the P-unit model. Each row illustrates subsequent processing steps for three different stimulation regimes: (i) baseline activity without external stimulus, only the fish's self-generated EOD (the carrier, \eqnref{eq:eod}) is present. (ii) RAM stimulation, \eqnref{eq:ram_equation}. The amplitude of the EOD carrier is modulated with a weak (2\,\% contrast) band-limited white-noise stimulus. (iii) Noise split, \eqnsref{eq:ram_split}--\eqref{eq:Noise_split_intrinsic}, where 90\,\% of the intrinsic noise is replaced by a RAM stimulus, whose amplitude is scaled to maintain the mean firing rate and the CV of the ISIs of the model's baseline activity. As an example, simulations of the model for cell ``2012-07-03-ak'' are shown. \figitem{A} The stimuli are thresholded, \eqnref{eq:threshold2}, by setting all negative values to zero. \figitem{B} Subsequent low-pass filtering, \eqnref{eq:dendrite}, attenuates the carrier and carves out the AM signal. \figitem{C} Intrinsic Gaussian white-noise is added to the signals shown in \panel{B}. Note the reduced internal noise amplitude in the noise split (iii) condition. \figitem{D} Spiking output of the LIF model, \eqnsref{eq:LIF}--\eqref{spikethresh}, in response to the sum of \panel{B} and \panel{C}. \figitem{E} Power spectra of the LIF neuron's spiking activity. Both, baseline activity (\panel[i]{E}) and noise split (\panel[iii]{E}), have the same peaks in the response spectrum at $r$, $f_{EOD} - r$, $f_{EOD}$, and $f_{EOD} + r$. With RAM stimulation (\panel[ii]{E}), the peak at the baseline firing rate, $r$, is washed out.}
 \end{figure*}
@@ -423,7 +389,7 @@ To mimic the interaction with other fish, the EODs of a second or third fish wit
 \end{equation}
 For two fish, $c_2 = 0$.
 
-Random amplitude modulations (RAMs) were simulated by first generating the AM as a band-limited white noise stimulus $s(t)$. For this, random real and imaginary numbers were drawn from Gaussian distributions for each frequency component in the range from 0 to 300\,Hz in the Fourier domain \citep{Billah1990,Skorjanc2023}. By means of the inverse Fourier transform, the time course of the RAM stimulus, $s(t)$,  was generated. The input to the model was then
+Random amplitude modulations (RAMs) were simulated by first generating the AM as a band-limited white noise stimulus $s(t)$. For this, random real and imaginary numbers were drawn from normal distributions for each frequency component in the range from 0 to 300\,Hz in the Fourier domain \citep{Billah1990,Skorjanc2023}. By means of the inverse Fourier transform, the time course of the RAM stimulus, $s(t)$,  was generated. The input to the model was then
 \begin{equation}
   \label{eq:ram_equation}
   y(t) = (1+ s(t)) \cos(2\pi f_{EOD} t)
@@ -435,12 +401,14 @@ First, the input $y(t)$ is thresholded by setting negative values to zero:
   \label{eq:threshold2}
   \lfloor y(t) \rfloor_0 = \left\{ \begin{array}{rcl} y(t) & ; & y(t) \ge 0 \\ 0 & ; & y(t) < 0 \end{array} \right.
 \end{equation}
-(\subfigrefb{flowchart}{A}). This thresholds models the transfer function of the synapses between the primary receptor cells and the afferent. Together with a low-pass filter
+(\subfigrefb{flowchart}{A}). This threshold models the transfer function of the synapses between the primary receptor cells and the afferent \citep{Moser2016} and may also include nonlinear properties introduced by low-threshold Kv1 channels known to be present in the afference \citep{TroySmith2006,Nogueira2013}. Together with a low-pass filter
 \begin{equation}
   \label{eq:dendrite}
   \tau_{d} \frac{d V_{d}}{d t} = -V_{d}+  \lfloor y(t) \rfloor_{0}
 \end{equation}
-the threshold operation is required for extracting the amplitude modulation from the input \citep{Barayeu2024}. The low-pass filter models passive signal conduction in the afferent's dendrite (\subfigrefb{flowchart}{B}) and $\tau_{d}$ is the membrane time constant of the dendrite. Dendritic low-pass filtering was also necessary to reproduce the loose coupling of P-unit spikes to the EOD while maintaining high sensitivity at small amplitude modulations.
+the threshold operation is required for extracting the amplitude modulation from the input \citep{Barayeu2024}. As detailed in the discussion, the intricacies of the threshold nonlinearity \eqnref{eq:threshold2} and the low-pass filter \eqnref{eq:dendrite} do not matter here, as all cross-spectra were computed between the spiking response and the amplitude-modulation $s(t)$ --- and not the full input signal $y(t)$, \eqnref{eq:ram_equation}.
+
+The low-pass filter models passive signal conduction in the afferent's dendrite (\subfigrefb{flowchart}{B}) from the synapse to the spike initiation zone. $\tau_{d}$ is the effective membrane time constant of this part of the dendrite. Dendritic low-pass filtering was also necessary to reproduce the loose coupling of P-unit spikes to the EOD while maintaining high sensitivity at small amplitude modulations.
 
 The dendritic voltage $V_d(t)$ is then fed into a stochastic leaky integrate-and-fire (LIF) model with adaptation,
 \begin{equation}
@@ -464,7 +432,7 @@ Whenever the membrane voltage $V_m(t)$ crosses the spiking threshold $\theta=1$,
 
 The P-unit models were integrated by the Euler forward method with a time-step of $\Delta t = 0.05$\,ms. For each trial of a simulation, $V_{m}$ was drawn from a uniform distribution between 0 and 1 and the initial value of $A$ was jittered by adding a random number drawn from a normal distribution with standard deviation of 2\,\% of its initial value. Then the first 500\,ms of any simulation were discarded to remove remaining transients.
 
-The eight free parameters of the P-unit model $\beta$, $\tau_m$, $\mu$, $D$, $\tau_A$, $\Delta_A$, $\tau_d$, and $t_{ref}$, were fitted to both the baseline activity (baseline firing rate, CV of ISIs, serial correlation of ISIs at lag one, and vector strength of spike coupling to EOD) and the responses to step increases and decreases in EOD amplitude (onset and steady-state responses, effective adaptation time constant, \citealp{Benda2005}) of recorded P-units. Model parameters of all 39 cells are summarized in file \texttt{models.csv} of our \texttt{punitmodel} repository at \url{https://github.com/bendalab/punitmodel/tree/v1}.
+The eight free parameters of the P-unit model $\beta$, $\tau_m$, $\mu$, $D$, $\tau_A$, $\Delta_A$, $\tau_d$, and $t_{ref}$, were fitted to both the baseline activity (baseline firing rate, CV of ISIs, serial correlation of ISIs at lag one, and vector strength of spike coupling to EOD) and the responses to step increases and decreases in EOD amplitude (onset and steady-state responses, effective adaptation time constant, \citealp{Benda2005}) of recorded P-units. Model parameters of all 39 cells are summarized in file \texttt{models.csv} of our \texttt{punitmodel} repository at \secret{\url{https://github.com/bendalab/punitmodel/tree/v1}}.
 
 
 \subsection{Noise split}
@@ -484,36 +452,34 @@ Both, the reduced intrinsic noise and the RAM stimulus, need to replace the orig
 \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 $f_{\rm EOD} =664$\,Hz encounters fish with EOD frequencies $f_{1}=631$\,Hz and $f_{2}=797$\,Hz. Both foreign signals have the same strength relative to the own field amplitude (10\,\% contrast). Top row: Sketch of signal processing in the nonlinear system (black box). Second row: Interference of the receiver EOD with the EODs of other fish, bold line highlights the amplitude modulation. Third row: Respective spike trains of the recorded P-unit. Fourth row: Firing rate, estimated by convolution of the spike trains with a Gaussian kernel ($\sigma = 1$\,ms). Bottom row: Power spectrum of the firing rate. \figitem{A} Baseline condition: The cell is driven by the self-generated field alone. The baseline firing rate $r$ dominates the power spectrum of the firing rate ($f_{base} = 139$\,Hz). \figitem{B} The receiver's EOD and a foreign fish with an EOD frequency $f_{1}=631$\,Hz are present. EOD interference induces an amplitude modulation, referred to as beat. \figitem{C} The receiver and a fish with an EOD frequency $f_{2}=797$\,Hz are present. The resulting beat is faster as the difference between the individual frequencies is larger. \figitem{D} All three fish with the EOD frequencies $f_{\rm EOD}$, $f_1$ and $f_2$ are present. A second-order amplitude modulation occurs, commonly referred to as envelope. Nonlinear peaks occur at the sum and difference of the two beat frequencies in the power spectrum of the firing rate.
-  }
+  \includegraphics[width=\columnwidth]{twobeats}
+  \caption{\label{fig:twobeats} Nonlinearity in an electrophysiologically recorded P-unit of \textit{A. leptorhynchus} in a three-fish setting (cell identifier ``2021-08-03-ac"). Receiver with EOD frequency $f_{\rm EOD} =664$\,Hz encounters fish with EOD frequencies $f_{1}=631$\,Hz and $f_{2}=797$\,Hz. Both foreign signals have the same strength relative to the own field amplitude (10\,\% contrast). Top row: Sketch of signal processing in the nonlinear system (black box). Second row: Superposition of the receiver EOD with the EODs of other fish, colored line highlights the amplitude modulation. Third row: Three trials of spike trains of the recorded P-unit. Fourth row: Firing rate, estimated by convolution of the spike trains with a Gaussian kernel. Bottom row: Power spectrum of the spike trains. \figitem{A} Baseline condition: The cell is driven by the self-generated field alone. The baseline firing rate $r = 139$\,Hz dominates the power spectrum (blue circle). \figitem{B} The receiver's EOD and a single conspecific with an EOD frequency $f_{1}=631$\,Hz are present. Superposition of the two EODs induces a periodic amplitude modulation, referred to as beat, with beat frequency $\Delta f_1=33$\,Hz. The P-unit strongly responds to this beat (purple). \figitem{C} The receiver and a fish with an EOD frequency $f_{2}=797$\,Hz are present. The resulting beat $\Delta f_2=133$\,Hz is faster as the difference between the EOD frequencies is larger. The P-unit response to this faster beat is weaker (green). \figitem{D} All three fish with EOD frequencies $f_{\rm EOD}$, $f_1$, and $f_2$ are present. Additional peaks occur in the power spectrum of the spike response at the sum (orange) and difference (red) of the two beat frequencies, indicating nonlinear interactions between the two frequencies in the P-unit. Note, the spectrum of the raw signal (top row, gray) has power only at the three EOD frequencies $f_{EOD}$, $f_1$, and $f_2$.}
 \end{figure*}
 
 We explored a large set of electrophysiological data from primary afferents of the active and passive electrosensory system, P-units and ampullary cells \citep{Grewe2017, Hladnik2023}, that were recorded in the brown ghost knifefish \textit{Apteronotus leptorhynchus}. We re-analyzed this dataset to search for weakly nonlinear responses that have been predicted in previous theoretical work \citep{Voronenko2017}. Additional simulations of LIF-based models of P-unit spiking help to interpret the experimental findings in this theoretical framework. 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}). 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.
+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:twobeats}{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:twobeats}{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:twobeats}{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.
+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:twobeats}{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.
 
 
 \subsection{Linear and weakly nonlinear regimes}
 
 \begin{figure*}[p]
   \includegraphics[width=\columnwidth]{regimes}
-  \caption{\label{fig:regimes} Linear and nonlinear responses of a model P-unit in a three-fish setting in dependence on stimulus amplitudes. The model P-unit (identifier ``2018-05-08-ad'') was stimulated with two sine waves of equal amplitude (contrast) at difference frequencies $\Delta f_1=40$\,Hz and $\Delta f_2=228$\,Hz relative the receiver's EOD frequency. $\Delta f_2$ was set to match the baseline firing rate $r$ of the P-unit. \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. \figitem{A} At 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 moderately higher stimulus contrast, the peaks in the response spectrum at the two beat frequencies are larger. \figitem{C} At intermediate stimulus contrasts, nonlinear responses start to appear at the sum and the difference of the stimulus frequencies (orange and red marker). \figitem{D} At higher stimulus contrasts additional peaks appear in the power spectrum. \figitem{E} Amplitude of the linear (at  $\Delta f_1$ and  $\Delta f_2$) and nonlinear (at $\Delta f_2 - \Delta f_1$ and $\Delta f_1 + \Delta f_2$) 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. In the linear regime, below a stimulus contrast of about 1.2\,\% (left vertical line), the only peaks in the response spectrum are at the stimulus frequencies. In the weakly nonlinear regime up to a contrast of about 3.5\,\% peaks arise at the sum and the difference of the two stimulus frequencies. At stronger stimulation the amplitudes of these nonlinear responses deviate from the quadratic dependency on stimulus contrast.}
+  \caption{\label{fig:regimes} Linear and nonlinear responses of a model P-unit in a three-fish setting in dependence on stimulus amplitudes. The model P-unit (identifier ``2018-05-08-ad'') was stimulated with two sine waves of equal amplitude (contrast) at difference frequencies $\Delta f_1=40$\,Hz and $\Delta f_2=228$\,Hz relative the receiver's EOD frequency $f_{EOD}=656$\,Hz. $\Delta f_2$ was set to match the baseline firing rate $r$ of the P-unit. \figitem{A--D} Top row: 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. Second row: Spike raster of the model P-unit response. Third row: average firing rate and, bottom row: power spectrum estimated from 1000 trials of spike train responses. \figitem{A} At low stimulus contrasts the response to the two beat frequencies is linear. These frequencies are present in the response spectrum (green and purple circles), the peak at $\Delta f_2$ (purple) enhances the peak at baseline firing rate (blue). The largest peak, however, is always the one at the EOD frequency of the receiver (black), reflecting the locking of P-unit spikes to the receiver's own EOD. The peak at $f_{EOD}$ is also flanked by harmonics of $\Delta f_1$ (gray triangles), where $f_{EOD} + \Delta f_1 = f_1$ is the EOD frequency of the first fish.  \figitem{B} At moderately higher stimulus contrast, the peaks in the response spectrum at the two beat frequencies become larger. In addition, two peaks at $f_{EOD} - \Delta f_2$ and  $f_{EOD} + \Delta f_2$ appear (blue diamonds), the latter is the EOD frequency of the second fish, $f_2$. These peaks indicate nonlinear processes acting on the three EOD frequencies from which the beat frequencies are generated. \figitem{C} At intermediate stimulus contrasts, nonlinear responses start to appear at the sum and the difference of the beat frequencies (orange and red circles). Similarly, side peaks appear at $\pm \Delta f_1$ around $f_{EOD} \pm \Delta f_2$ (purple squares) as well as harmonics of $\Delta f_1$ (small green circles). \figitem{D} At higher stimulus contrasts many more peaks appear in the power spectrum. Most of them are interactions of $f_{EOD}$ and $\Delta f_2$ with harmonics of $\Delta f_1$ ($k=1, 2, ...$), since the cell responds strongest to $\Delta f_1$. \figitem{E} Amplitude of the linear (at  $\Delta f_1$ and  $\Delta f_2$) and nonlinear (at $\Delta f_2 - \Delta f_1$ and $\Delta f_1 + \Delta f_2$) 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. In the linear regime, below a stimulus contrast of about 1.2\,\% (left vertical line), the only peaks in the response spectrum are at the stimulus frequencies. In the weakly nonlinear regime up to a contrast of about 3.5\,\% peaks arise at the sum and the difference of the two stimulus frequencies. At stronger stimulation the amplitudes of these nonlinear responses deviate from the quadratic dependency on stimulus contrast.}
 \end{figure*}
 
-The stimuli used in \figref{fig:motivation} had the same not-small amplitude. Whether this stimulus condition 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}.
+The stimuli used in \figref{fig:twobeats} had the same AM contrasts of 10\,\% that evoke rather strong modulations in a P-unit's firing rate response. Whether these definitely not-small stimulus amplitudes fall 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 (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.
+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 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. The largest peak at the receiver's EOD frequency indicates the stochastic locking of the P-unit spikes to the EOD, and its side-peaks are remnants of the two stimulating EOD frequencies \citep{Sinz2020}.
 
-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.
+This linear regime is followed by the weakly nonlinear regime (in the example cell between approximately 1.2\,\% and 2.6\,\% stimulus contrast). In addition to the peaks at the stimulating beat frequencies, peaks at the sum and the difference of the beat 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 $r$ 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.
 
-For this example, we chose very 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 the sake of simplicity we will drop the $\Delta$ notation even though P-unit stimuli are beats.
+For this example, we chose very specific stimulus (beat) frequencies. 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 the sake of simplicity we will drop the $\Delta$ notation even though P-unit stimuli are beats.
 
 
 \subsection{Nonlinear signal transmission in P-units}
@@ -522,12 +488,12 @@ P-units fire action potentials probabilistically phase-locked to the self-genera
 
 \begin{figure*}[p]
   \includegraphics[width=\columnwidth]{punitexamplecell} 
-  \caption{\label{fig:punit} Linear and nonlinear stimulus encoding in example P-units. \figitem{A} Interspike interval (ISI) distribution of a cell's baseline activity, i.e. the cell is driven only by the unperturbed own electric field (cell identifier ``2020-10-27-ag''). This cell has a rather high baseline firing rate $r=405$\,Hz and an intermediate CV$_{\text{base}}=0.49$ of its interspike intervals. \figitem{B} Power spectral density of the cell's baseline response with marked peaks at the cell's baseline firing rate $r$ and the fish's EOD frequency $f_{\text{EOD}}$. \figitem{C} Random amplitude modulation (RAM) stimulus (top, red, with cutoff frequency of 300\,Hz) and evoked responses (spike raster, bottom) of the same P-unit for two different stimulus contrasts (right). The stimulus contrast quantifies the standard deviation of the RAM relative to the fish's EOD amplitude. \figitem{D} Gain of the transfer function (first-order susceptibility), \eqnref{linearencoding_methods}, computed from the responses to 10\,\% (light blue) and 20\,\% contrast (dark blue) RAM stimulation of 5\,s duration. \figitem{E} Absolute value of the second-order susceptibility, \eqnref{chi2}, for both the low and high stimulus contrast. At the lower stimulus contrast an anti-diagonal where the sum of the two stimulus frequencies equals the neuron's baseline frequency clearly sticks out of the noise floor. \figitem{F} At the higher contrast, the anti-diagonal is much weaker. \figitem{G} Second-order susceptibilities projected onto the diagonal (averages over all anti-diagonals of the matrices shown in \panel{E, F}). The anti-diagonals from \panel{E} and \panel{F} show up as a peak close to the cell's baseline firing rate $r$. The susceptibility index, SI($r$) \eqnref{siindex}, quantifies the height of this peak relative to the values in the vicinity. \figitem{H} ISI distributions (top) and second-order susceptibilities (bottom) of two more example P-units (``2021-06-18-ae'', ``2017-07-18-ai'') showing an anti-diagonal, but not the full expected triangular structure. \figitem{I} Most P-units, however, have a flat second-order susceptibility and consequently their SI($r$) values are close to one (cell identifiers ``2018-08-24-ak'', ``2018-08-14-ac'').}
+  \caption{\label{fig:punit} Linear and nonlinear stimulus encoding in example P-units. \figitem{A} Interspike interval (ISI) distribution of a cell's baseline activity, i.e. the cell is driven only by the unperturbed own electric field (cell identifier ``2020-10-27-ag''). This cell has a rather high baseline firing rate $r=405$\,Hz and an intermediate CV$_{\text{base}}=0.49$ of its interspike intervals. \figitem{B} Power spectral density of the cell's baseline response with marked peaks at the cell's baseline firing rate $r$ and the fish's EOD frequency $f_{\text{EOD}}$. \figitem{C} Random amplitude modulation (RAM) stimulus (top, red, with cutoff frequency of 300\,Hz) and evoked responses (spike raster, bottom) of the same P-unit for two different stimulus contrasts (right). The stimulus contrast quantifies the standard deviation of the RAM relative to the fish's EOD amplitude. \figitem{D} Gain of the transfer function (first-order susceptibility), \eqnref{linearencoding_methods}, computed from the responses to 10\,\% (light blue) and 20\,\% contrast (dark blue) RAM stimulation of 5\,s duration. \figitem{E} Absolute value of the second-order susceptibility, \eqnref{chi2}, for the low stimulus contrast. An anti-diagonal where the sum of the two stimulus frequencies equals the neuron's baseline frequency clearly sticks out of the noise floor. \figitem{F} At the higher contrast, the anti-diagonal in the absolute value of the second-order susceptibility is much weaker. \figitem{G} Second-order susceptibilities projected onto the diagonal (averages over all anti-diagonals of the matrices shown in \panel{E, F}). The anti-diagonals from \panel{E} and \panel{F} show up as a peak close to the cell's baseline firing rate $r$. The susceptibility index, SI($r$) \eqnref{siindex}, quantifies the height of this peak relative to the values in the vicinity. \figitem{H} ISI distributions (top) and second-order susceptibilities (bottom) of two more example P-units (``2021-06-18-ae'', ``2017-07-18-ai'') showing an anti-diagonal, but not the full expected triangular structure. \figitem{I} Most P-units, however, have a flat second-order susceptibility and consequently their SI($r$) values are close to one (cell identifiers ``2018-08-24-ak'', ``2018-08-14-ac'').}
 \end{figure*}
 
 Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:punit}{C}, top trace, red line), have been 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 from existing recordings 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 blue for low and high contrast stimuli, respectively, \subfigrefb{fig:punit}{C}). Linear encoding, quantified by the first-order susceptibility or 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{chi2},  quantifies for each combination of two stimulus frequencies $f_1$ and $f_2$ the amplitude and phase of the stimulus-evoked response at the sum $f_1 + f_2$ (and also the difference, \subfigrefb{fig:lifsuscept}{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 $f_1$ and $f_2$ or their sum $f_1 + f_2$ exactly match the neuron's baseline firing rate $r$ \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:lifsuscept}{B}).
+The second-order susceptibility, \eqnref{chi2},  quantifies for each combination of two stimulus frequencies $f_1$ and $f_2$ the amplitude and phase of the stimulus-evoked response at the sum $f_1 + f_2$ (and also the difference, \subfigrefb{fig:lifsuscept}{B}). 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 $f_1$ and $f_2$ or their sum $f_1 + f_2$ exactly match the neuron's baseline firing rate $r$ \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:lifsuscept}{B}).
 
 For the example P-unit, we observe a ridge of elevated second-order susceptibility for the low RAM contrast at $f_1 + f_2 = r$ (yellowish anti-diagonal, \subfigrefb{fig:punit}{E}). This structure is less prominent 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 (white dashed line) by averaging over the anti-diagonals (\subfigrefb{fig:punit}{G}). At low RAM contrast this projection indeed has a distinct peak close to the neuron's baseline firing rate (\subfigrefb{fig:punit}{G}, dot on top line). For the higher RAM contrast this peak is much smaller 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.
 
@@ -545,7 +511,7 @@ 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 where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, but not where either og the frequencies alone matches the baseline rate. In the following, we investigate this discrepancy to the theoretical expectations \citep{Voronenko2017,Franzen2023}.
+In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampullary}), we only observe nonlinear responses where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, but not where either of the frequencies alone matches the baseline rate. In the following, we investigate this discrepancy to the theoretical expectations (\figrefb{fig:lifsuscept}, \citealp{Voronenko2017,Franzen2023}).
 
 \begin{figure*}[p]
   \includegraphics[width=\columnwidth]{noisesplit}
@@ -554,7 +520,7 @@ In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampull
 
 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_1 + 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}).
 
@@ -564,7 +530,7 @@ In the model, however, we know the time course of the intrinsic noise and can us
 
 \begin{figure}[p]
   \includegraphics[width=\columnwidth]{modelsusceptcontrasts}
-  \caption{\label{fig:modelsusceptcontrasts}Dependence of second order susceptibility on stimulus contrast. \figitem{A} Second-order susceptibilities estimated for increasing stimulus contrasts of $c=0, 1, 3$ and $10$\,\% as indicated ($N=10^7$ FFT segments for $c=1$\,\%, $N=10^6$ segments for all other contrasts). $c=0$\,\% refers to the noise-split configuration (limit to vanishing external RAM signal, see \subfigrefb{fig:noisesplit}{D}). Shown are simulations of the P-unit model cell ``2017-07-18-ai'') with a baseline firing rate of $82$\,Hz and CV$_{\text{base}}=0.23$. The cell shows a clear triangular pattern in its second-order susceptibility even up to a contrast of $10$\,\%. Note, that for $c=1$\,\% (\panel[ii]{D}), the estimate did not converge yet. \figitem{B} Cell ``2012-12-13-ao'' (baseline firing rate of $146$\,Hz, CV$=0.23$) also has strong interactions at its baseline firing rate that survive up to a stimulus contrast of $3$\,\%. \figitem{C} Model cell ``2012-12-20-ac'' (baseline firing rate of $212$\,Hz, CV$=0.26$) shows a weak triangular structure in the second-order susceptibility that vanishes for stimulus contrasts larger than $1$\,\%. \figitem{D} Cell ``2013-01-08-ab'' (baseline firing rate of $218$\,Hz, CV$=0.55$) does not show any triangular pattern in its second-order susceptibility. Nevertheless, interactions between low stimulus frequencies become substantial at higher contrasts. \figitem{E} The presence of an elevated second-order susceptibility where the stimulus frequency add up to the neuron's baseline frequency, can be identified by the susceptibility index (SI($r$), \eqnref{siindex}) greater than one (horizontal black line). The SI($r$) (density to the right) is plotted as a function of the model neuron's baseline CV for all $39$ model cells. Model cells have been visually categorized based on the strong (11 cells) or weak (5 cells) presence of a triangular pattern in their second-order susceptibility estimated in the noise-split configuration (legend). The cells from \panel{A--D} are marked by black circles. Pearson's correlation coefficients $r$, the corresponding significance level $p$ and regression line (dashed gray line) are indicated. The higher the stimulus contrast, the less cells show weakly nonlinear interactions as expressed by the triangular structure in the second-order susceptibility.}
+  \caption{\label{fig:modelsusceptcontrasts}Dependence of second order susceptibility on stimulus contrast. \figitem{A} Second-order susceptibilities estimated for increasing stimulus contrasts of $c=0, 1, 3$ and $10$\,\% as indicated ($N=10^7$ FFT segments for $c=1$\,\%, $N=10^6$ segments for all other contrasts). $c=0$\,\% refers to the noise-split configuration (limit to vanishing external RAM signal, see \subfigrefb{fig:noisesplit}{D}). Shown are simulations of the P-unit model cell ``2017-07-18-ai'') with a baseline firing rate of $82$\,Hz and CV$_{\text{base}}=0.23$. The cell shows a clear triangular pattern in its second-order susceptibility even up to a contrast of $10$\,\%. Note, that for $c=1$\,\% (\panel[ii]{D}), the estimate did not converge yet. \figitem{B} Cell ``2012-12-13-ao'' (baseline firing rate of $146$\,Hz, CV$=0.23$) also has strong interactions at its baseline firing rate that survive up to a stimulus contrast of $3$\,\%. \figitem{C} Model cell ``2012-12-20-ac'' (baseline firing rate of $212$\,Hz, CV$=0.26$) shows a weak triangular structure in the second-order susceptibility that vanishes for stimulus contrasts larger than $1$\,\%. \figitem{D} Cell ``2013-01-08-ab'' (baseline firing rate of $218$\,Hz, CV$=0.55$) does not show any triangular pattern in its second-order susceptibility. Nevertheless, interactions between low stimulus frequencies become substantial at higher contrasts. \figitem{E} The presence of an elevated second-order susceptibility where the stimulus frequency add up to the neuron's baseline frequency, can be identified by the susceptibility index (SI($r$), \eqnref{siindex}) greater than one (horizontal black line). The SI($r$) (density to the right) is plotted as a function of the model neuron's baseline CV for all $39$ model cells. Model cells have been visually categorized based on the strong (11 cells) or weak (5 cells) presence of a triangular pattern in their second-order susceptibility estimated in the noise-split configuration (legend). The cells from \panel{A--D} are marked by black circles. Pearson's correlation coefficients $r$, the corresponding significance level $p$ and regression line (dashed gray line) are indicated. The higher the stimulus contrast, the fewer cells show weakly nonlinear interactions as expressed by the triangular structure in the second-order susceptibility.}
 \end{figure}
 
 \subsection{Weakly nonlinear interactions in many model cells}
@@ -578,13 +544,15 @@ The SI($r$) correlates with the CVs of the cell's baseline interspike intervals
 
 
 \subsection{Weakly nonlinear interactions vanish for higher stimulus contrasts}
-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.
+As pointed out above, noise stimuli act as a linearizing background noise and thus may shift the neural signal transmission into a regime free of nonlinearities (see discussion). Sine-wave stimuli, as for example used in the introductory figure~\ref{fig:regimes}, do not have such a linearizing effect. At which amplitudes a noise stimulus effectively linearizes the system, also depends on cellular properties of the neuron, or on the parameters in our leaky integrate-and-fire neuron. This is what we explore now by varying the contrast of RAM stimuli driving our model cells.
+
+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\,\%.
 
 \begin{figure}[tp]
   \includegraphics[width=\columnwidth]{modelsusceptlown}
-  \caption{\label{fig:modelsusceptlown}Inferring the triangular structure of the second-order susceptibility from limited data. \figitem{A} Reliably estimating the structure of the second-order susceptibility requires a high number of FFT segments $N$ in the order of one or even ten millions. As an example, susceptibilities of the model cell ``2012-12-21-ak-invivo-1'' (baseline firing rate of 157\,Hz, CV=0.15) are shown for the noise-split configuration ($c=0$\,\%) and RAM stimulus contrasts of $c=1$, $3$, and $10$\,\% as indicated. For contrasts below $10$\,\% this cell shows a nice triangular pattern in its susceptibilities, quite similar to the introductory example of a LIF in \figrefb{fig:lifsuscept}. \figitem{B} However, with limited data of $N=100$ trials the susceptibility estimates are noisy and show much less structure, except for the anti-diagonal at the cell's baseline firing rate. The SI($r$) quantifies the height of this ridge where the two stimulus frequencies add up to the neuron's baseline firing rate. \figitem{C} Correlations between the estimates of SI($r$) based on 100 FFT segments (density to the right) with the converged ones based on one or ten million segments at a given stimulus contrast for all $n=39$ model cells. The black circle marks the model cell shown in \panel{A} and \panel{B}. The black diagonal line is the identity line and the dashed line is a linear regression. The correlation coefficient and corresponding significance level are indicated in the top left corner. The thin vertical line is a threshold at 1.2, the thin horizontal line a threshold at 1.8. The number of cells within each of the resulting four quadrants denote the false positives (top left), true positives (top right), true negatives (bottom left), and false negatives (bottom right) for predicting a triangular structure in the converged susceptibility estimate from the estimates based on only 100 segments.}
+  \caption{\label{fig:modelsusceptlown}Inferring the triangular structure of the second-order susceptibility from limited data. \figitem{A} Reliably estimating the structure of the second-order susceptibility requires a high number of FFT segments $N$ in the order of one or even ten millions. As an example, susceptibilities of the model cell ``2012-12-21-ak'' (baseline firing rate of 157\,Hz, CV=0.15) are shown for the noise-split configuration ($c=0$\,\%) and RAM stimulus contrasts of $c=1$, $3$, and $10$\,\% as indicated. For contrasts below $10$\,\% this cell shows a nice triangular pattern in its susceptibilities, quite similar to the introductory example of a LIF in \figrefb{fig:lifsuscept}. \figitem{B} However, with limited data of $N=100$ trials the susceptibility estimates are noisy and show much less structure, except for the anti-diagonal at the cell's baseline firing rate. The SI($r$) quantifies the height of this ridge where the two stimulus frequencies add up to the neuron's baseline firing rate. \figitem{C} Correlations between the estimates of SI($r$) based on 100 FFT segments (density to the right) with the converged ones based on one or ten million segments at a given stimulus contrast for all $n=39$ model cells. The black circle marks the model cell shown in \panel{A} and \panel{B}. The black diagonal line is the identity line and the dashed line is a linear regression. The correlation coefficient and corresponding significance level are indicated in the top left corner. The thin vertical line is a threshold at 1.2, the thin horizontal line a threshold at 1.8. The number of cells within each of the resulting four quadrants denote the false positives (top left), true positives (top right), true negatives (bottom left), and false negatives (bottom right) for predicting a triangular structure in the converged susceptibility estimate from the estimates based on only 100 segments.}
 \end{figure}
 
 \subsection{Weakly nonlinear interactions can be deduced from limited data}
@@ -648,7 +616,7 @@ In the P-unit models, each model cell contributed with three RAM stimulus presen
 
 The effective stimulus strength also plays a role in predicting the SI($r$) values. We quantify the effect of stimulus strength on a cell's response by the response modulation --- the standard deviation of a cell's firing rate in response to a RAM stimulus. The lower the response modulation, i.e. the weaker the effective stimulus, the higher the S($r$) (\figrefb{fig:dataoverview}\,\panel[ii]{A}). Although there is a tendency of low stimulus contrasts to evoke lower response modulations, response modulations evoked by each of the three contrasts overlap substantially, emphasizing the strong heterogeneity of the P-units' sensitivity \citep{Grewe2017}. Cells with high SI($r$) values are the ones with baseline firing rate below 200\,Hz (\figrefb{fig:dataoverview}\,\panel[iii]{A}).
 
-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 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=14051$, $p=4\times 10^{-9}$), because in the measured cells, response modulation is positively correlated with baseline CV (Pearson $R=0.34$, $p=1\times 10^{-4}$), 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$ 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}).
 
@@ -658,7 +626,7 @@ The population of ampullary cells is more homogeneous, with generally lower base
 
 \section{Discussion}
 
-Theoretical work \citep{Voronenko2017,Franzen2023} derived analytical expressions for weakly-nonlinear responses of spike generating 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 evidence for weakly-nonlinear responses of a spike generator at low intrinsic noise levels or low stimulus amplitudes in real sensory neurons.
+Theoretical work \citep{Voronenko2017,Franzen2023} studied analytically and numerically the weakly-nonlinear responses of spike generating 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 evidence for weakly-nonlinear responses of a spike generator at low intrinsic noise levels or low stimulus amplitudes in real sensory neurons.
 
 \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 structure of the second-order susceptibility (\citealp{Voronenko2017}, \subfigref{fig:lifsuscept}{B}). The CV of the baseline interspike interval 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. Still, only 18\,\% of the P-units analyzed in this study show relevant nonlinear responses. On the other hand, the majority (74\,\%) of the ampullary cells show nonlinear responses as they have generally lower CVs (median of 0.09).
@@ -679,9 +647,9 @@ Making assumptions about the nonlinearities in a system also reduces the amount
 
 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, 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 stimuli originate from the muscle potentials induced by prey movement \citep{Kalmijn1974, Engelmann2010, Neiman2011fish}. For a single prey item 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 (\figrefb{fig:dataoverview}~\panel[ii]{C}) indicating that paddlefish data have been recorded above the weakly-nonlinear regime. 
+Ampullary stimulus encoding is somewhat different in \textit{A. leptorhynchus}. 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 (\figrefb{fig:dataoverview}~\panel[ii]{C}) 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.
 
@@ -689,15 +657,24 @@ The population of ampullary cells is very homogeneous with respect to the baseli
 
 In contrast to the ampullary cells, P-units respond to the amplitude modulation of the self-generated EOD. Extracting the AM requires a (threshold) nonlinearity \citep{Middleton2006, Stamper2012Envelope, Savard2011, Barayeu2023}. This nonlinearity, however, does not show up in our estimates of the susceptibilities, because in our analysis we directly relate the AM waveform to the recorded cellular responses. Encoding the time-course of the AM, however, has been shown to be linear over a wide range of AM amplitudes and frequencies \citep{Xu1996, Benda2005, Gussin2007, Grewe2017, Savard2011}. In contrast, we here have demonstrated nonlinear interactions originating from the spike generator for broad-band noise stimuli with small amplitudes and for stimulation with two distinct frequencies. Both settings have not been studied yet.
 
-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? 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 and substantial power (\figrefb{fig:motivation}) although these interactions vanish when stimulating with noise stimuli of similar contrast (\figrefb{fig:modelsusceptcontrasts}).
+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? Broadband noise stimuli introduce additional noise that linearizes the dynamics of the system. In contrast, a 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 and substantial power (\figrefb{fig:twobeats}) although these interactions vanish when stimulating with noise stimuli of similar contrast (\figrefb{fig:modelsusceptcontrasts}).
 
-The encoding of secondary AMs or social envelopes that arise from relative movement or the interaction of more than two animals \citep{Stamper2012Envelope} requires an another nonlinearity in addition to the one needed for extracting the AM. Initially, this nonlinearity was attributed to downstream processing \citep{Middleton2006, Middleton2007}. Later studies showed that already the electroreceptors can encode such information whenever the firing rate saturates at zero or the maximum rate at the EOD frequency \citep{Savard2011}. Based on our work, we predict that P-units with low CVs encode the social envelopes even under weak stimulation, whenever the resulting beat frequencies match or add up to the baseline firing rate. Then difference frequencies show up in the response spectrum that characterize the slow envelope.
+The encoding of secondary AMs or social envelopes that arise from relative movement or the interaction of more than two animals \citep{Stamper2012Envelope} requires another nonlinearity in addition to the one needed for extracting the AM. Initially, this nonlinearity was attributed to downstream processing \citep{Middleton2006, Middleton2007}. Later studies showed that already the electroreceptors can encode such information whenever the firing rate saturates at zero or the maximum rate at the EOD frequency \citep{Savard2011}. Based on our work, we predict that P-units with low CVs encode the social envelopes even under weak stimulation, whenever the resulting beat frequencies match or add up to the baseline firing rate. Then difference frequencies show up in the response spectrum that characterize the slow envelope.
 
 The weakly nonlinear interactions in low-CV P-units could facilitate the detectability of faint signals during three-animal interactions as found in courtship contexts in the wild \citep{Henninger2018}. The detection of a faint, distant intruder male could be improved by the presence of a nearby strong female stimulus, because of the nonlinear interaction terms \citep{Schlungbaum2023}. This boosting effect is, however, very specific with respect to the stimulus frequencies and a given P-unit's baseline frequency. The population of P-units is very heterogeneous in their baseline firing rates and CVs (50--450\,Hz and 0.1--1.4, respectively, \subfigref{fig:dataoverview}{B}, \citealp{Grewe2017, Hladnik2023}). The range of baseline firing rates thus covers substantial parts of the beat frequencies that may occur during animal interactions \citep{Henninger2018, Henninger2020}, while the number of P-units showing weakly nonlinear responses is small. Whether and how this information is specifically maintained and read out by pyramidal cells in the electrosensory lateral line lobe (ELL) in the hindbrain onto which P-units converge \citep{Krahe2014, Maler2009a} is an open question.
 
+Electric fish are able to slowly modulate their EOD frequency, as for example during the so-called jamming-avoidance-response \citep{Fortune2020}. Such behaviors modify the resulting beat frequency by a few Hertz. This could in principle increase the chance that the now slowly changing beat frequency matches at some point the baseline firing rate of a P-unit, where the weakly nonlinear responses then enhance the detectability of another conspecific \citep{Schlungbaum2023}. Furthermore, transient changes in EOD frequency on timescales of tens of milliseconds up to a few seconds are known as chirps and rises, respectively, and are involved in courtship and aggression behaviors \citep{Henninger2018, Raab2021}. How the encoding of such transient frequency modulations is affected by the nonlinearities described here is another open question, since the presented analysis focuses on stationary signals.
+
 \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, 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.
+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. These could boost or distort responses to two simultaneously presented tones and thus might play a role in forming perception of music.
+
+\paragraph{Acknowledgements}
+We thank Tim Hladnik, Henriette Walz, Franziska Kuempfbeck, Fabian Sinz, Laura Seidler, Eva Vennemann, and Ibrahim Tunc for the data they recorded over the years in our lab.
+
+\paragraph{Founding sources}
+Supported by SPP 2205 ``Evolutionary optimisation of neuronal processing'' by the DFG, project number 430157666,
+and by the Open Access Publication Fund of the University of T\"ubingen.
 
 
 %\bibliographystyle{apalike}%alpha}%}%alpha}%apalike}
@@ -718,4 +695,4 @@ We have demonstrated pronounced nonlinear responses in primary electrosensory af
 % LocalWords:  Furutsu Novikov interspike durations Apteronotus
 % LocalWords:  leptorhynchus nonlinearity linearizing spectrally
 % LocalWords:  electroreceptors lowpass differentially isopotential
-% LocalWords:  transdermal convolved
+% LocalWords:  transdermal convolved electrophysiologically

plotstyle.py

@@ -125,7 +125,7 @@ def plot_style():
     pt.make_line_styles(ns, 'ls', 'AMsplit', '', palette['orange'], '-', lwmid)
     pt.make_line_styles(ns, 'ls', 'Noise', '', palette['gray'], '-', lwmid)
     pt.make_line_styles(ns, 'ls', 'Median', '', palette['black'], '-', lwthick)
-    pt.make_line_styles(ns, 'ls', 'Max', '', palette['black'], '-', lwmid)
+    pt.make_line_styles(ns, 'ls', 'Max', '', palette['black'], '--', lwthin)
 
     ns.psC1 = dict(color=palette['red'], marker='o', linestyle='none', markersize=3, mec='none', mew=0)
     ns.psC3 = dict(color=palette['orange'], marker='o', linestyle='none', markersize=3, mec='none', mew=0)
@@ -135,24 +135,31 @@ def plot_style():
 
     ns.lsStim = dict(color=palette['gray'], lw=ns.lwmid)
     ns.lsRaster = dict(color=palette['black'], lw=ns.lwthin)
+    ns.lsRate = dict(color=palette['blue'], lw=ns.lwmid)
     ns.lsPower = dict(color=palette['gray'], lw=ns.lwmid)
-    ns.lsF0 = dict(color='blue', lw=ns.lwthick)
-    ns.lsF01 = dict(color='green', lw=ns.lwthick)
-    ns.lsF02 = dict(color='purple', lw=ns.lwthick)
-    ns.lsF012 = dict(color='orange', lw=ns.lwthick)
-    ns.lsF01_2 = dict(color='red', lw=ns.lwthick)
-    ns.lsF0m = dict(color=lighter('blue', 0.5), lw=ns.lwthin)
-    ns.lsF01m = dict(color=lighter('green', 0.6), lw=ns.lwthin)
-    ns.lsF02m = dict(color=lighter('purple', 0.5), lw=ns.lwthin)
-    ns.lsF012m = dict(color=darker('orange', 0.9), lw=ns.lwthin)
-    ns.lsF01_2m = dict(color=darker('red', 0.9), lw=ns.lwthin)
+    ns.lsF0 = dict(color=palette['blue'], lw=ns.lwthick)
+    ns.lsF01 = dict(color=palette['green'], lw=ns.lwthick)
+    ns.lsF02 = dict(color=palette['purple'], lw=ns.lwthick)
+    ns.lsF012 = dict(color=palette['orange'], lw=ns.lwthick)
+    ns.lsF01_2 = dict(color=palette['red'], lw=ns.lwthick)
+    ns.lsF0m = dict(color=lighter(palette['blue'], 0.5), lw=ns.lwthin)
+    ns.lsF01m = dict(color=lighter(palette['green'], 0.6), lw=ns.lwthin)
+    ns.lsF02m = dict(color=lighter(palette['purple'], 0.5), lw=ns.lwthin)
+    ns.lsF012m = dict(color=darker(palette['orange'], 0.9), lw=ns.lwthin)
+    ns.lsF01_2m = dict(color=darker(palette['red'], 0.9), lw=ns.lwthin)
 
-    ns.psFEOD = dict(color='black', marker='o', linestyle='none', markersize=5, mec='none', mew=0)
-    ns.psF0 = dict(color='blue', marker='o', linestyle='none', markersize=5, mec='none', mew=0)
-    ns.psF01 = dict(color='green', marker='o', linestyle='none', markersize=5, mec='none', mew=0)
-    ns.psF02 = dict(color='purple', marker='o', linestyle='none', markersize=5, mec='none', mew=0)
-    ns.psF012 = dict(color='orange', marker='o', linestyle='none', markersize=5, mec='none', mew=0)
-    ns.psF01_2 = dict(color='red', marker='o', linestyle='none', markersize=5, mec='none', mew=0)
+    ns.psFEOD = dict(color=palette['black'], marker='o', linestyle='none', markersize=4, mec='none', mew=0)
+    ns.psF0 = dict(color=palette['blue'], marker='o', linestyle='none', markersize=4, mec='none', mew=0)
+    ns.psF01 = dict(color=palette['green'], marker='o', linestyle='none', markersize=4, mec='none', mew=0)
+    ns.psF02 = dict(color=palette['purple'], marker='o', linestyle='none', markersize=4, mec='none', mew=0)
+    ns.psF012 = dict(color=palette['orange'], marker='o', linestyle='none', markersize=4, mec='none', mew=0)
+    ns.psF01_2 = dict(color=palette['red'], marker='o', linestyle='none', markersize=4, mec='none', mew=0)
+    ns.psFEODm = dict(color=lighter(palette['black'], 0.5), marker='^', linestyle='none', markersize=3, mec='none', mew=0)
+    ns.psF0m = dict(color=lighter(palette['blue'], 0.6), marker='D', linestyle='none', markersize=2.5, mec='none', mew=0)
+    ns.psF01m = dict(color=lighter(palette['green'], 0.6), marker='o', linestyle='none', markersize=3, mec='none', mew=0)
+    ns.psF02m = dict(color=lighter(palette['purple'], 0.6), marker='s', linestyle='none', markersize=2.4, mec='none', mew=0)
+    ns.psF012m = dict(color=lighter(palette['orange'], 0.6), marker='p', linestyle='none', markersize=3, mec='none', mew=0)
+    ns.psF01_2m = dict(color=palette['red'], marker='o', linestyle='none', markersize=3, mec='none', mew=0)
     
     ns.model_color1 = palette['purple']
     ns.model_color2 = lighter(ns.model_color1, 0.6)

punitexamplecell.py

@@ -141,11 +141,13 @@ def plot_response(ax, s, eodf, time1, stimulus1, contrast1, spikes1,
         eod = np.sin(2*np.pi*eodf*time1) * (1 + stim)
     else:
         eod = np.sin(2*np.pi*eodf*time1) + stim
-    ax.plot(time1, 4*eod + 7, **s.lsEOD)
-    ax.plot(time1, 4*(1 + stim) + 7, **s.lsAM)
+    ax.plot(time1, 4*eod + maxtrials - 1, **s.lsEOD)
+    ax.plot(time1, 4*(1 + stim) + maxtrials - 1, **s.lsAM)
     ax.set_xlim(t0, t1)
     ax.set_ylim(-2*maxtrials - 0.5, 14)
     ax.xscalebar(1, -0.05, 0.01, None, '10\\,ms', ha='right')
+    ax.text(t1 + 0.003, maxtrials - 1 + 4, f'${100*contrast1:.0f}$\\,\\%',
+            va='center', color=s.lsAM['color'])
     ax.text(t1 + 0.003, -0.5*maxtrials, f'${100*contrast1:.0f}$\\,\\%',
             va='center', color=s.cell_color1)
     ax.text(t1 + 0.003, -1.55*maxtrials, f'${100*contrast2:.0f}$\\,\\%',
@@ -179,7 +181,7 @@ def plot_diagonals(ax, s, fbase, contrast1, freqs1, chi21,
         sis.append(sinorm)
         sips.append(sip)
         sifs.append(sif)
-        print(f'    SI at {100*contrast:.1f}% contrast: {sinorm:.2f}')
+        print(f'    SI at {100*contrast:4.1f}% contrast: {sinorm:.2f}')
     ax.plot(diags[1][0], 1e-4*diags[1][1], **s.lsC2)
     ax.plot(diags[0][0], 1e-4*diags[0][1], **s.lsC1)
     offs = 0.05*ymax

rebuttal2.tex

@@ -0,0 +1,406 @@
+\documentclass[11pt]{article}
+
+\usepackage[utf8]{inputenc}
+\usepackage{textcomp}
+\usepackage{xcolor}
+\usepackage{graphicx}
+
+\usepackage[ngerman,english]{babel}
+
+\usepackage[left=25mm, right=25mm, top=20mm, bottom=20mm]{geometry}
+\setlength{\parskip}{2ex}
+
+\usepackage[mediumqspace,Gray,squaren]{SIunits}      % \ohm, \micro
+
+\usepackage{natbib}
+%\bibliographystyle{jneurosci}
+
+\usepackage[breaklinks=true,bookmarks=true,bookmarksopen=true,pdfpagemode=UseNone,pdfstartview=FitH,colorlinks=true,citecolor=blue,urlcolor=blue]{hyperref}
+
+\newcommand{\issue}[1]{\textbf{#1}}
+\newcommand{\issueg}[1]{\foreignlanguage{ngerman}{\textbf{#1}}}
+\newcounter{responsecounter}
+\newcommand{\response}[1]{\refstepcounter{responsecounter}\begin{quote}\arabic{responsecounter}. #1\end{quote}}
+\newcommand{\note}[2][]{\textcolor{red!80!black}{\textbf{[#1: #2]}}}
+\newcommand{\notejb}[1]{\note[JB]{#1}}
+\newcommand{\notejg}[1]{\note[JG]{#1}}
+\newcommand{\notesr}[1]{\note[SR]{#1}}
+
+\newcommand{\changed}[1]{\textcolor{blue!50!black}{#1}}
+
+\setlength{\parindent}{0em}
+
+
+\begin{document}
+
+We would like to thank both reviewers for their valuable
+feedback. Note that line numbers mentioned in our following responses refer to
+the new version of the manuscript, not the redlined one.
+
+\issue{\large Reviewer \#1}
+
+\issue{The manuscript "Spike generation in electroreceptor afferents
+  introduces additional spectral response components by weakly
+  nonlinear interactions" submitted to eNeuro represents a noteworthy
+  advancement in the field, as it elucidates that, under an often
+  naturally occurring scenario the non linear responses of the pair
+  electroreceptor-primary afferent ensemble may intervene in signal
+  encoding. The manuscript shows that during the reception of signals
+  originating distantly from multiple individual conspecifics,
+  electroreceptor primary afferents may exhibit nonlinear responses
+  allowing the fish to guess the presence of more than one
+  individual. This is articulated with clarity through a
+  straightforward leaky-integrate-and-fire model of electroreceptor
+  responsiveness. The illustrations are both lucid and enhance the
+  comprehension of the results section. The discussion is
+  sound. However, the intrinsic value of the manuscript would likely
+  be obscure without a more "biologist-friendly" approach. I would
+  like to offer several suggestions that may serve to either enhance
+  the manuscript or inspire future research endeavors.}
+
+\response{Thank you for trying to make our manuscript more
+  biologist-friendly! And yes, some of your comments indeed inspired
+  our thinking for future projects.}
+
+\issue{First, I should point out that beyond the presence of a
+  threshold-induced nonlinearity, the complex structure of the
+  axon-like dendritic innervation receptor cell terminals within the
+  electroreceptor organ. This analogical nonlinear response may have
+  its origin in the branched anatomy of the dendrite-axon terminals,
+  easily verifiable by anatomical studies, and the presence, hardly
+  demonstrable but plausible, of ion channel diversity; see, for
+  example, Trigo, F. F. (2019) Antidromic analog signaling. Frontiers
+  in Cellular Neuroscience, 13, 354. for a discussion of the general
+  case and the study by Troy Smith, Unguez and Weber (2006, Fig. 3) in
+  which receptor cells of tuberous electroreceptor organs and their
+  afferents from Apteronotus leptorhynchus were labeled to varying
+  degrees by six anti-Kv1 antibodies. Kv1.1 and Kv1.4 immunoreactivity
+  was intense in the afferent axons of electroreceptor organs. It is
+  noteworthy that Kv1 are low-threshold channels and, in some cases,
+  exhibit a prolonged refractory period (Nogueira and Caputi,
+  2013). These sources of nonlinearity could be mentioned to
+  strengthen the links between well-written theoretical analysis and
+  the practical field of experimental physiology.}
+
+\response{You are right, there are more nonlinear mechanisms
+  potentially contributing to the threshold nonlinearity. We now
+  mention this in the methods when introducing the threshold
+  nonlinearity (after eq. 13) and cite the corresponding
+  articles.}
+
+\issue{Second, and along the same lines, the discussion could be
+  improved by mentioning the effects and significance of these
+  nonlinearities when the recipient fish makes changes in its EOD
+  frequency in at least two cases: a) sustained changes, as in
+  interference avoidance responses, and b) transient changes, as in
+  chirps.}
+
+\response{We added a paragraph addressing JARs, chirps, and rises to
+  the discussion (lines 697--705).}
+
+\issue{Finally, the precise description of the methods could be
+  expanded for reaching a broader biology audience; in particular, the
+  purpose of some procedures should be explained in some way. While
+  the meaning seems clear as the reader scrolls through the results, a
+  first reading of the methods, although accurate, does not offer the
+  biology reader a quick and intuitive approach to the study.}
+
+\issue{Next, I list some minor more detailed comments that may clarify
+  the design and methods and facilitate their understanding by a
+  broader audience.}
+
+\issue{In general, you referenced P receptors and ampullary
+  structures; however, what about T receptors? How can one distinguish
+  between T and P in the recordings? Might it be possible that the
+  negative results observed in certain receptors are attributable to
+  the type of receptor (P or T)? Did you postulate, as suggested by
+  Viancour (1979), that there exists a continuum of responsiveness
+  between the extreme profiles of P (signal amplitude) and T (signal
+  slope)?}
+
+\response{T-units are characterized by 1:1 locking to the EOD, i.e. by
+  having a baseline firing rate matching the EOD frequency. We
+  definitely have no T-units in our data set, since our P-unit firing
+  rates are well below the EOD frequencies. This we explain now in the
+  ``Identification of P-units and ampullary cells'' section in the
+  methods.}
+
+\issue{In line 147, rather than using the term
+  "laterally," I believe it would enhance clarity to state "parallel
+  to each side of the fish," as the orientation of the electrodes may
+  otherwise remain ambiguous.}
+
+\response{Done.}
+
+\issue{Furthermore, no commentary or discussion is provided regarding
+  the fact that the stimulation procedure, which is transverse to the
+  main axis of the body, neglects to account for the effects on the
+  field foveal perioral region where the majority of receptors are
+  located.}
+
+\response{As stated in ``Experimental subjects and procedures'', all
+  recordings were done in the posterior lateral line nerve. So we did
+  not record from the foveal perioral region, and hence this problem
+  is not relevant.}
+
+\issue{Line 148, the phrase "band limited white noise" lacks
+  clarity. Upon my initial reading, I assumed that the cutoff limit
+  you referenced pertained to a low pass filter applicable to both
+  ampullary and P-type tuberous receptors; however, it could indeed be
+  interpreted as the opposite. In a strict sense, all "white noise
+  stimuli" are band-passed. The duration of the stimulus establishes a
+  lower cutoff for the band pass in one instance, while the
+  responsiveness of the stimulation apparatus delineates the upper
+  cutoff limits in another. Nevertheless, once one comprehends the
+  objective of the experiment, the implicit significance of the white
+  noise filtering becomes exceedingly apparent. Thus, this description
+  could benefit from greater clarity to avoid the need to explore the
+  results first in order to understand well.}
+
+\response{We are sorry for the confusion. The cutoff frequencies
+  stated are pure stimulus parameters and not related to the filtering
+  performed by the respective neurons. ``White noise'' refers to a
+  time series that has equal power at all frequencies (like white
+  light) --- this choice of signal is agnostic with respect to the
+  preferred time scales of the system because all frequencies (or,
+  timescales) appear equally on the stimulus side. Bandpass-limited
+  white noise has equal power at all frequencies up to a cutoff
+  frequency that the experimenters choose in order to distribute the
+  total power over a reasonable frequency range in which they expect a
+  measurable response of the system under investigation. The choice
+  was different for ampullary receptors and P-units as stated in the
+  manuscript, but the stated values are not related with the actual
+  bandpass filtering that the neurons perform on the input
+  stimulus. The latter are quantified in the paper when we look at the
+  linear and nonlinear response functions of the cells. We completely
+  rewrote the description of the white-noise stimuli in the methods
+  sections (lines 155--160).}
+
+\issue{Line 154. This procedure elicits a modulation of the envelope
+  of the reafferent signal. To achieve this, you adopted distinct
+  approaches for the ampullary and P receptors: a) in the case of
+  ampullary receptors, you presented white noise and incrementally
+  elevated its amplitude (variance) until the mean amplitude of the
+  averaged sine wave recorded via local electrodes adjacent to the
+  gills exhibited an increase of 1 to 5\%, is this correct?}
+
+\response{We increased the amplitude of the white noise until the
+  standard deviation (not the mean) of the resulting modulation of the
+  EOD reached 1 to 5\,\%. We rephrased the description of the
+  stimulation and hope that this is clearer now (lines 166--169).}
+
+\issue{b) with regard to P receptors, you multiplied the head-to-tail
+  ongoing signal by a white noise signal and played the resultant
+  output, adjusting the amplitude until the local signal experienced
+  an enhancement of 1 to 5\% in average, is this interpretation
+  accurate?  Since the head to tail EOD and the local signals over the
+  body are out of phase this process induces both amplitude and phase
+  modulation of the stimulus signal, which will be contingent upon the
+  phase lag of the local EOD at the receptor site in relation to the
+  head-to-tail EOD. This phase lag, as reported in the literature,
+  exhibits a shift ranging from pi to 2pi between a receptor situated
+  at the head and another at the tail. (I posit that this may not
+  significantly impact individual receptors response; however, how
+  does this influence the relative timing among distinct receptors,
+  and what is its correlation with jamming avoidance mechanisms?)
+  Furthermore, does this form of noise modulation exert a comparable
+  effect on the flanks (i.e., the apex of the derivative) as it does
+  on the peaks of the signal themselves? How does this affect the
+  recruitment of P and T receptors?}
+
+\response{You are right about the phase shifts and that this does not
+  ``significantly impact individual receptors response''. This is a
+  standard stimulation procedure for characterizing receptor responses
+  that are located mainly on the sides of a fish's flat body as the
+  recording site is on the posterior branch of the lateral line
+  nerve. See, for example, Hladnik and Grewe, 2023. And yes, this will
+  probably impact relative spike timing in distinct receptors and thus
+  may also impact the JAR mechanisms. However, this manuscript is
+  about single receptor responses and not about T-units, and we feel
+  it is already complicated enough. Therefore we would rather prefer
+  to not open up all these issues, since they are not relevant for the
+  results we present.}
+
+\issue{Line 238. Are you referring to the terminal non-myelinated
+  branches that connect receptor cells to the initial Ranvier node?
+  The peripheral afferent constitutes a myelinated and active dendrite
+  whose distal branches receive synapses from receptor
+  cells. Consequently, there exists a summation occurring at some
+  juncture, likely at the first node, that facilitates the generation
+  of an action potential. Otherwise, the signal would not be
+  effectively propagated from the receptors to the ganglia where the
+  somata reside. Receptors across various species exhibit notable
+  differences; some are myelinated within the electroreceptor organ,
+  while others display the first node external to the electroreceptor
+  organ. Could you discuss this aspect, considering the anatomical
+  structure of the receptor in your species?}
+
+\response{Exactly. We slightly expanded our description to make clear
+  that we talk about the signal transduction until it reaches the
+  spike initiation zone (lines 260--261).}
+
+
+\issue{\large Reviewer \#2}
+
+\issue{This work is a nice contribution to our general understanding
+  of nonlinearities in sensory coding, and to our detailed
+  understanding of behaviorally relevant information processing in the
+  electrosensory systems of weakly electric fish. I have several
+  suggestions for the authors.}
+
+\issue{(1) Abstract, line 29. "...if these frequencies or their sum
+  match the neuron's baseline firing rate" is not quite accurate
+  because "these frequencies" implies BOTH input frequencies must
+  match the baseline firing rate. I think you mean to state, "...if
+  one of these frequencies or their sum match the neuron's baseline
+  firing rate."}
+
+\response{Your are right! We changed the sentence as suggested.}
+
+\issue{(2) Abstract, line 33. The wording here is unclear,
+  specifically what you mean by "much stronger." Much stronger what
+  exactly? I think you mean to refer to the fact that these nonlinear
+  responses were more common and stronger in ampullary units than
+  P-units, but "much stronger" does not clearly convey this,
+  especially in the abstract.}
+
+\response{We changed the sentence to ``...  identify these predicted
+  nonlinear responses only in individual low-noise P-units, but in
+  more than half of the ampullary cells.''}
+
+\issue{(3) Figure 1A. "r" needs to be clearly defined here. Based on
+  the text, it seems to be the baseline firing rate of the neuron, but
+  this needs to be made clear in the figure legend.}
+
+\response{We added a brief explanatory sentence to the caption.}
+
+\issue{(4) Figure 1B. "Because frequencies can also be negative..."
+  This is unclear and needs more explanation, especially because there
+  are no negative frequencies in your actual data. How can frequencies
+  be negative?}
+
+\response{We added a few sentences following equation (1) to motivate
+  the existence of negative frequencies in Fourier transforms. And we
+  added a hint in the caption of figure 1B.}
+
+\issue{(5) Figure 3 and 4. Why are the power spectra clipped at such
+  low frequencies? This makes it impossible to see peaks due to
+  potential df2 harmonics and fEOD. Figure 5 extends to higher
+  frequencies to illustrate these and it is not clear why these are
+  clipped in these two figures.}
+
+\response{You are right. In figure 4 we show now the spectrum up to
+  950\,Hz, such that $f_{EOD}$ and its interactions with $f_1$ and
+  $f_2$ are included. We labeled the additional peaks and expanded the
+  figure caption accordingly. In figure 3 we stay with the small
+  range, because we have so little data for this special setting where
+  one of the beat frequencies approximately matches the P-units
+  baseline firing rate (only three trials of 500ms duration). This is
+  why the power spectra are very noisy. Also, for an introductory
+  figure we prefer to only show the few peaks that are relevant for
+  the rest of the manuscript, to not overwhelm the reader right at the
+  start.}
+
+\issue{(6) Figure 3. Why are these example firing rates based on
+  convolution with a 1 ms Gaussian kernel if the analyses were based
+  on convolution with a 2 ms Gaussian kernel (line 169)? It seems that
+  example data should effectively illustrate how the data were
+  actually analyzed. More fundamentally, why would a 2-fold difference
+  in kernel width be appropriate for presentation vs. analysis?}
+
+\response{Thank you for addressing this inconsistency. This was for
+  ``historical'' reasons. We now decided to use the 1\,ms kernel for
+  all figures and analysis. We changed the sentence in the methods
+  accordingly (line 183). In doing so we also added panels showing
+  firing rates in addition to the response spectra in figure 4. Using
+  the more narrow kernel better reveals the details of the time course
+  of the firing rate and this way improves the connection between the
+  firing rate and the response spectra. In figure 10, middle column,
+  the range of possible values of the response modulations is a bit
+  enlarged by using the 1\,ms kernel, but the correlations and their
+  significance did not change a lot either.}
+
+\issue{(7) Figure 3D legend. The relationship between 2nd order AM
+  (envelope) and the two nonlinear peaks should be made clear. I
+  believe the envelope is represented by both peaks, correct?}
+
+\response{No, what is shown is the power spectrum of the spike
+  response, not the one of the amplitude modulation or envelope of the
+  stimulus. We added a sentence to the end of the figure caption to
+  make this clear.\\
+  If it were the power spectrum of the signal after it passed
+  a non-linearity (rectification or thresholding at zero), then there
+  could be also peaks at the sum and difference of the beat
+  frequencies. However, since they are close to the higher one of the
+  two beat frequencies they do not show up in the AM as obviously as
+  for the settings used in the social envelope papers by Eric Fortune
+  and Andre Longtin and colleges (we guess this is what you had in
+  mind).}
+
+\issue{(8) Line 302. "not-small amplitude" is arbitrary and
+  vague. Please be clearer and more precise.}
+
+\response{We rephrased to two sentences in lines 325--327.}
+
+\issue{(9) Figures 5C and 6C. For the stimuli with the red RAM
+  waveforms, please make it clear which contrast is being represented
+  by these traces, as responses to two different contrasts are shown.}
+
+\response{We added the contrast values to both figures.}
+
+\issue{(10) Figure 5E, F. The legend states that second-order
+  susceptibility for both the low and high stimulus contrasts are
+  shown in E, but E shows the low contrast and F shows the high
+  contrast.}
+
+\response{Good catch! Fixed.}
+
+\issue{(11) Lines 453-465, Figure 8. This section was confusing to
+  me. Why does second-order susceptibility decrease as stimulus
+  contrast increases, when theory predicts that higher signal-to-noise
+  ratios should result in larger nonlinearities?}
+
+\response{Yes, in figure 4 increasing stimulus contrast results in
+  stronger nonlinearities. There, the stimuli are narrow-band sine
+  waves. However, as pointed out in the context of figure 7, when
+  using a broad-band noise stimulus instead, this stimulus by itself
+  adds background noise to the system that linearizes the response. In
+  this context, it is crucial to realize that the (linear and
+  nonlinear) transfer of a nonlinear system like a neuron depends on
+  the background noise. A Gaussian noise stimulus acts here both (i)
+  as a signal that evokes a response (linear and nonlinear) and (ii)
+  as an additional background noise linearizing the (linear and
+  nonlinear) response. In the context of our study it implies the
+  susceptibilities estimated from noise stimuli decrease for higher
+  stimulus contrasts.\\
+  We added a whole paragraph at the beginning of this section to make
+  this clear (lines 479--484).}
+
+\issue{(12) Lines 655-675. This was a very nice end to the discussion,
+  but I would like to see more. I would like the broader significance
+  of this study to be expanded upon with respect to (1) behavioral
+  relevance for signal detection in weakly electric fish, and (2)
+  comparative relevance for other modalities and species. Speculation
+  is fine so long as it is clearly indicated as such. It might work
+  best to expand upon and distribute the information in lines 655-666
+  throughout the discussion at relevant points, rather than as an
+  afterthought. The conclusion section in lines 667-675 could then
+  reiterate these points briefly and delve into more detail on
+  comparative considerations.}
+
+\response{We also like to see more on this, but we feel that we
+  already speculated enough. Without further studies on the readout of
+  the receptor responses, we cannot make any convincing claim about
+  whether and how weakly nonlinear interactions are actually utilized
+  in a neural system. The problem is that a match of one of the
+  stimulating frequencies or their sum with the neuron's baseline
+  firing rate is required. This is all addressed in the (now second
+  final) paragraph of the ``Nonlinear encoding in P-units'' section.
+  Nevertheless, in response to reviewer \#1, we added another
+  paragraph discussing various behaviors that modulate the EOD
+  frequency and how these may exploit the weakly nonlinear
+  interactions.  However, we agree that the comparative aspect of the
+  conclusion could be expanded. We therefore added one more final
+  speculative sentence to the conclusion (lines 715--716).}
+
+\end{document}

rebuttal3.tex

@@ -0,0 +1,52 @@
+\documentclass[11pt]{article}
+
+\usepackage[utf8]{inputenc}
+\usepackage{textcomp}
+\usepackage{xcolor}
+\usepackage{graphicx}
+
+\usepackage[ngerman,english]{babel}
+
+\usepackage[left=25mm, right=25mm, top=20mm, bottom=20mm]{geometry}
+\setlength{\parskip}{2ex}
+
+\usepackage[mediumqspace,Gray,squaren]{SIunits}      % \ohm, \micro
+
+\usepackage{natbib}
+%\bibliographystyle{jneurosci}
+
+\usepackage[breaklinks=true,bookmarks=true,bookmarksopen=true,pdfpagemode=UseNone,pdfstartview=FitH,colorlinks=true,citecolor=blue,urlcolor=blue]{hyperref}
+
+\newcommand{\issue}[1]{\textbf{#1}}
+\newcommand{\issueg}[1]{\foreignlanguage{ngerman}{\textbf{#1}}}
+\newcounter{responsecounter}
+\newcommand{\response}[1]{\refstepcounter{responsecounter}\begin{quote}\arabic{responsecounter}. #1\end{quote}}
+\newcommand{\note}[2][]{\textcolor{red!80!black}{\textbf{[#1: #2]}}}
+\newcommand{\notejb}[1]{\note[JB]{#1}}
+\newcommand{\notejg}[1]{\note[JG]{#1}}
+\newcommand{\notesr}[1]{\note[SR]{#1}}
+
+\newcommand{\changed}[1]{\textcolor{blue!50!black}{#1}}
+
+\setlength{\parindent}{0em}
+
+
+\begin{document}
+
+Thank you for your final positive assessment of our manuscript.
+
+\issue{\large Reviewer \#1}
+
+\issue{Line 18 In the abstract the sentence "...only in individual low-noise P-units, but in more than half of the ampullary cells". is not quite clear. Please rephrase.}
+
+\issue{\large Reviewer \#2}
+
+\issue{The authors have done an excellent job of responding to my comments. I will note that the line numbers referenced in the response letter are incorrect, which made it challenging to review these changes. However, I have just one remaining minor comment. In the abstract, line 18, the wording is still unclear. "...more than half of the ampullary cells" is clear, but "...only in individual low-noise P-units" is not. Why not precisely state that these predicted nonlinear responses were observed in X out of Y P-units?}
+
+\response{We have rephrased this sentence as requested: ``In our
+  combined experimental and modeling approach we identify these
+  predicted nonlinear responses in those 31 out of 172 P-units that
+  are characterized by low intrinsic noise. In contrast, the majority
+  (22 out of 30) ampullary cells show nonlinear responses.''}
+
+\end{document}

references.bib

@@ -4178,6 +4178,17 @@ We collected weakly electric gymnotoid fish in the vicinity of Manaus, Amazonas,
   Timestamp                = {2020.02.11}
 }
 
+
+@article{Fortune2020,
+  title={Spooky interaction at a distance in cave and surface dwelling electric fishes.},
+  author={Eric S. Fortune and Nicole Andanar and Manu Madhav and Ravikrishnan P. Jayakumar and Noah J. Cowan and Maria Elina Bichuette and Daphne Soares},
+  journal={Front Integr Neurosci},
+  volume={14},
+  pages={561524},
+  year={2020}
+}
+
+
 @article{Maier2008,
   title={Integration of bimodal looming signals through neuronal coherence in the temporal lobe},
   author={Maier, Joost X and Chandrasekaran, Chandramouli and Ghazanfar, Asif A},
@@ -5158,7 +5169,7 @@ and Keller, Clifford H.},
   publisher={American Psychological Association}
 }
 
-@Article{Raab2019,
+@Article{Raa2019,
   Title                    = {Dominance in Habitat Preference and Diurnal Explorative Behavior of the Weakly Electric Fish \textit{Apteronotus leptorhynchus}},
   Author                   = {Raab, Till and Linhart, Laura and Wurm, Anna and Benda, Jan},
   Journal                  = {Frontiers in Integrative Neuroscience},
@@ -5173,6 +5184,15 @@ and Keller, Clifford H.},
   Url                      = {https://www.frontiersin.org/article/10.3389/fnint.2019.00021}
 }
 
+@article{Raab2021,
+  title={Electrocommunication signals indicate motivation to compete during dyadic interactions of an electric fish.},
+  author={Till Raab and Sercan Bayezit and Saskia Erdle and Jan Benda},
+  journal={J Exp Biol},
+  volume={224},
+  pages={jeb242905},
+  year={2021}
+}
+
 @Article{Rao1999,
   Title                    = {Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects.},
   Author                   = {Rajesh P.N. Rao and Dana H. Ballard},
@@ -7113,3 +7133,30 @@ microelectrode recordings from visual cortex and functional implications.},
   pages={145--156},
   year={1981},
 }
+
+@article{Nogueira2013,
+  title={From the intrinsic properties to the functional role of a neuron phenotype: an example from electric fish during signal trade-off.},
+  author={Javier Nogueira and Angel A. Caputi},
+  journal={J Exp Biol},
+  volume={216},
+  pages={2380--2392},
+  year={2013},
+}
+
+@article{Nogueira2014,
+  title={Pharmacological study of the one spike spherical neuron phenotype in {\textit{Gymnotus omarorum}}.},
+  author={Javier Nogueira and Angel A. Caputi},
+  journal={Neuroscience},
+  volume={258},
+  pages={347--354},
+  year={2014},
+}
+
+@article{TroySmith2006,
+  title={Distribution of {Kv1}-like potassium channels in the electromotor and electrosensory systems of the weakly electric fish {\textit{Apteronotus leptorhynchus}}.},
+  author={G. Troy Smith and Graciela A. Unguez and Christopher M. Weber},
+  journal={J Neurobiol},
+  volume={66},
+  pages={1011--1031},
+  year={2006},
+}

regimes.py

@@ -1,16 +1,39 @@
 import numpy as np
 import matplotlib.pyplot as plt
+
 from pathlib import Path
 from scipy.stats import linregress
 from numba import jit
+
+from spectral import rate
 from plotstyle import plot_style, lighter, darker
 
 
 model_cell = '2018-05-08-ad-invivo-1'      # 228Hz, CV=0.67
+#alphas = [0.002, 0.01, 0.03, 0.06]
+alphas = [0.002, 0.01, 0.025, 0.05]
+rmax = 500
+amax = 60
+cthresh1 = 1.2
+cthresh2 = 2.6
+
+#model_cell = '2018-05-08-ab-invivo-1'      # 116, CV=0.68
+#alphas = [0.002, 0.008, 0.025, 0.05]
+#rmax = 400
+#amax = 50
+#cthresh1 = 1.2
+#cthresh2 = 3.0
 
 data_path = Path('data')
 sims_path = data_path / 'simulations'
 
+trials = 1000
+spec_trials = 1000   # set to zero to only recompute firng rates
+sigma = 0.001
+nfft = 2**18
+
+recompute = False
+
 
 def load_data(file_path):
     data = np.load(file_path)
@@ -127,13 +150,13 @@ def punit_spikes(parameter, alpha, beatf1, beatf2, tmax, trials):
     return spikes
 
     
-def plot_am(ax, s, alpha, beatf1, beatf2, tmax):
-    time = np.arange(0, tmax, 0.0001)
+def plot_am(ax, s, alpha, beatf1, beatf2, tmin, tmax):
+    time = np.arange(tmin, tmax, 0.0001)
     am = alpha*np.sin(2*np.pi*beatf1*time)
     am += alpha*np.sin(2*np.pi*beatf2*time)
     ax.show_spines('l')
-    ax.plot(1000*time, -100*am, **s.lsAM)
-    ax.set_xlim(0, 1000*tmax)
+    ax.plot(1000*(time - tmin), -100*am, **s.lsAM)
+    ax.set_xlim(0, 1000*(tmax - tmin))
     ax.set_ylim(-13, 13)
     ax.set_yticks_delta(10)
     #ax.set_xlabel('Time', 'ms')
@@ -142,22 +165,47 @@ def plot_am(ax, s, alpha, beatf1, beatf2, tmax):
             transform=ax.transAxes, ha='right')
 
     
-def plot_raster(ax, s, spikes, tmax):
-    spikes_ms = [1000*s[s<tmax] for s in spikes[:16]]
+def plot_raster(ax, s, spikes, tmin, tmax):
+    spikes_ms = [1000*(s[(s >= tmin) & (s <= tmax)] - tmin)
+                 for s in spikes[:16]]
     ax.show_spines('')
     ax.eventplot(spikes_ms, linelengths=0.9, **s.lsRaster)
-    ax.set_xlim(0, 1000*tmax)
+    ax.set_xlim(0, 1000*(tmax - tmin))
     #ax.set_xlabel('Time', 'ms')
     #ax.set_ylabel('Trials')
 
+    
+def plot_rate(ax, s, path, spikes, tmin, tmax, sigma=0.002):
+    if recompute or not path.is_file():
+        print('    compute firing rate')
+        time = np.arange(0, tmin + tmax, sigma/4)
+        r, rsd = rate(time, spikes, sigma)
+        np.savez(path, time=time, rate=r, ratesd=rsd,
+                 sigma=sigma, trials=len(spikes))
+    else:
+        print(f'    load firing rate from {path}')
+        data = np.load(path)
+        time = data['time']
+        r = data['rate']
+        rsd = data['ratesd']
+    mask = (time >= tmin) &(time <= tmax)
+    time = time[mask] - tmin
+    r = r[mask]
+    ax.show_spines('l')
+    ax.plot(1000*time, r, clip_on=False, **s.lsRate)
+    ax.set_xlim(0, 1000*(tmax - tmin))
+    ax.set_ylim(0, rmax)
+    ax.set_ylabel('Rate', 'Hz')
+    ax.set_yticks_delta(200)
 
-def compute_power(path, contrast, spikes, nfft, dt):
-    if not path.exists():
-        print(f'Compute power spectrum for contrast = {100*contrast:4.1f}%')
+
+def compute_power(path, spikes, nfft, dt):
+    if spec_trials > 0 and (recompute or not path.is_file()):
+        print('    compute power spectrum')
         psds = []
         time = np.arange(nfft + 1)*dt
         tmax = nfft*dt
-        for s in spikes:
+        for s in spikes[:spec_trials]:
             b, _ = np.histogram(s, time)
             b = b / dt
             fourier = np.fft.rfft(b - np.mean(b))
@@ -165,9 +213,9 @@ def compute_power(path, contrast, spikes, nfft, dt):
             freqs = np.fft.rfftfreq(nfft, dt)
             prr = np.mean(psds, 0)*dt/nfft
         np.savez(path, nfft=nfft, deltat=dt, nsegs=len(spikes),
-                 freqs=freqs, prr=prr)
+                 freqs=freqs, prr=prr, trials=len(spikes))
     else:
-        print(f'Load power spectrum for contrast = {100*contrast:4.1f}%')
+        print(f'    load power spectrum from {path}')
         data = np.load(path)
         freqs = data['freqs']
         prr = data['prr']
@@ -175,45 +223,105 @@ def compute_power(path, contrast, spikes, nfft, dt):
 
 
 def decibel(x):
-    return 10*np.log10(x/1e8)
+    return 10*np.log10(x/1e8 + 1e-12)
 
 
-def plot_psd(ax, s, path, contrast, spikes, nfft, dt, beatf1, beatf2):
-    offs = 4
-    freqs, psd = compute_power(path, contrast, spikes, nfft, dt)
+def peak_ampl(freqs, psd, f, df=2):
+    if f < 0:
+        f = 5
+    psd_snippet = psd[(freqs > f - df) & (freqs < f + df)]
+    return np.max(psd_snippet)
+
+
+def plot_psd(ax, s, path, contrast, spikes, nfft, dt, beatf1, beatf2, eodf):
+    offs = 5
+    offsm = 3
+    freqs, psd = compute_power(path, spikes, nfft, dt)
     psd /= freqs[1]
     ax.plot(freqs, decibel(psd), **s.lsPower)
+    # mark frequencies:
+    ax.plot(eodf, decibel(peak_ampl(freqs, psd, eodf)) + offs,
+            label=r'$f_{EOD}$', clip_on=False, zorder=50, **s.psFEOD)
+    if contrast < alphas[1]:
+        ax.text(eodf + 40, -2, '$f_{EOD}$', va='center', color=s.psFEOD['color'])
+        ax.text(eodf + beatf1 + 5, -33, '$f_1$', color=s.psFEOD['color'])
+        #ax.text(beatf1 + 10, -30, r'$\Delta f_1$', ha='center', color=s.psF01['color'])
+        ax.text(beatf2 + 10, -24, '$r$', ha='center', color=s.psF0['color'])
+        #ax.text(beatf2 + 10, -16, r'$\Delta f_2$', ha='center', color=s.psF02['color'])
     ax.plot(beatf2, decibel(peak_ampl(freqs, psd, beatf2)) + offs,
-            label=r'$r$', clip_on=False, **s.psF0)
+            label=r'$r$', clip_on=False, zorder=50, **s.psF0)
     ax.plot(beatf1, decibel(peak_ampl(freqs, psd, beatf1)) + offs,
-            label=r'$\Delta f_1$', clip_on=False, **s.psF01)
-    ax.plot(beatf2, decibel(peak_ampl(freqs, psd, beatf2)) + 2*offs + 3,
-            label=r'$\Delta f_2$', clip_on=False, **s.psF02)
-    ax.plot(beatf2 - beatf1, decibel(peak_ampl(freqs, psd, beatf2 - beatf1)) + offs,
-            label=r'$\Delta f_2 - \Delta f_1$', clip_on=False, **s.psF01_2)
-    ax.plot(beatf1 + beatf2, decibel(peak_ampl(freqs, psd, beatf1 + beatf2)) + offs,
-            label=r'$\Delta f_1 + \Delta f_2$', clip_on=False, **s.psF012)
-    ax.set_xlim(0, 300)
+            label=r'$\Delta f_1$', clip_on=False, zorder=50, **s.psF01)
+    ax.plot(beatf2, decibel(peak_ampl(freqs, psd, beatf2)) + 2*offs + 2,
+            label=r'$\Delta f_2$', clip_on=False, zorder=50, **s.psF02)
+    ax.plot(eodf + beatf1, decibel(peak_ampl(freqs, psd, eodf + beatf1)) + offsm,
+            label=r'$f_{EOD} \pm k \Delta f_1$', zorder=40, **s.psFEODm)
+    ax.plot(eodf - beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf1)) + offsm, **s.psFEODm)
+    if contrast >= alphas[1]:
+        ax.plot(eodf - beatf2, decibel(peak_ampl(freqs, psd, eodf - beatf2)) + offsm,
+                label=r'$f_{EOD} \pm \Delta f_2$', zorder=40, **s.psF0m)
+        ax.plot(eodf + beatf2, decibel(peak_ampl(freqs, psd, eodf - beatf2)) + offsm,
+                zorder=40, **s.psF0m)
+        if contrast < alphas[2]:
+            ax.text(eodf + beatf2 + 10, -25, '$f_2$', ha='center', color=s.psF0['color'])
+    if contrast >= alphas[2]:
+        ax.plot(beatf2 - beatf1, decibel(peak_ampl(freqs, psd, beatf2 - beatf1)) + offs,
+                label=r'$\Delta f_2 - \Delta f_1$', clip_on=False, zorder=50, **s.psF01_2)
+        ax.plot(beatf1 + beatf2, decibel(peak_ampl(freqs, psd, beatf1 + beatf2)) + offs,
+                label=r'$\Delta f_1 + \Delta f_2$', clip_on=False, zorder=50, **s.psF012)
+        ax.plot(2*beatf1, decibel(peak_ampl(freqs, psd, 2*beatf1)) + offsm,
+                label=r'$k\Delta f_1$', zorder=40, **s.psF01m)
+        ax.plot(eodf + 2*beatf1, decibel(peak_ampl(freqs, psd, eodf + 2*beatf1)) + offsm, zorder=40, **s.psFEODm)
+        ax.plot(eodf - beatf2 + beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 + beatf1)) + offsm,
+                label=r'$f_{EOD} \pm \Delta f_2 \pm \Delta f_1$', zorder=40, **s.psF02m)
+        ax.plot(eodf - beatf2 - beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 - beatf1)) + offsm,
+                zorder=40, **s.psF02m)
+        ax.plot(eodf + beatf2 + beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 - beatf1)) + offsm,
+                zorder=40, **s.psF02m)
+        ax.plot(eodf + beatf2 - beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 - beatf1)) + offsm,
+                zorder=40, **s.psF02m)
+    if contrast >= alphas[3]:
+        ax.plot(beatf2 + 2*beatf1, decibel(peak_ampl(freqs, psd, beatf2 + 2*beatf1)) + offsm,
+                label=r'$\Delta f_2 \pm k\Delta f_1$', zorder=40, **s.psF012m)
+        ax.plot(beatf2 + 3*beatf1, decibel(peak_ampl(freqs, psd, beatf2 + 3*beatf1)) + offsm, zorder=40, **s.psF012m)
+        ax.plot(beatf2 - 2*beatf1, decibel(peak_ampl(freqs, psd, beatf2 - 2*beatf1)) + offsm, zorder=40, **s.psF012m)
+        ax.plot(beatf2 - 3*beatf1, decibel(peak_ampl(freqs, psd, beatf2 - 3*beatf1)) + offsm, zorder=40, **s.psF012m)
+        ax.plot(3*beatf1, decibel(peak_ampl(freqs, psd, 3*beatf1)) + offsm, zorder=40, **s.psF01m)
+        ax.plot(4*beatf1, decibel(peak_ampl(freqs, psd, 4*beatf1)) + offsm, zorder=40, **s.psF01m)
+        ax.plot(eodf - 2*beatf1, decibel(peak_ampl(freqs, psd, eodf - 2*beatf1)) + offsm, zorder=40, **s.psFEODm)
+        ax.plot(eodf - 3*beatf1, decibel(peak_ampl(freqs, psd, eodf - 3*beatf1)) + offsm, zorder=40, **s.psFEODm)
+        #ax.plot(eodf - 4*beatf1, decibel(peak_ampl(freqs, psd, eodf - 4*beatf1)) + offsm, zorder=40, **s.psFEODm)
+        #ax.plot(eodf - 2*beatf2, decibel(peak_ampl(freqs, psd, eodf - 2*beatf2)) + offsm, zorder=40, **s.psF0m)
+        #ax.plot(eodf - beatf2 + 2*beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 + 2*beatf1)) + offsm,
+        #        zorder=40, **s.psF02m)
+        #ax.plot(eodf - beatf2 + 3*beatf1, decibel(peak_ampl(freqs, psd, eodf - beatf2 + 2*beatf1)) + offsm,
+        #        zorder=40, **s.psF02m)
+    ax.set_xlim(0, 950)
     ax.set_ylim(-60, 0)
+    ax.set_xticks_delta(200)
+    ax.set_yticks_delta(20)
     ax.set_xlabel('Frequency', 'Hz')
     ax.set_ylabel('Power [dB]')
 
 
-def plot_example(axs, axr, axp, s, path, cell, alpha, beatf1, beatf2,
+def plot_example(axs, axr, axf, axp, s, path, cell, alpha, beatf1, beatf2, eodf,
                  nfft, trials):
+    spec_path = path.with_name(path.stem + f'-contrastspectrum-c{1000*alpha:03.0f}.npz')
+    rate_path = path.with_name(path.stem + f'-contrastrates-c{1000*alpha:03.0f}-s{1000*sigma:.0f}.npz')
     dt = 0.0001
     tmax = nfft*dt
-    t1 = 0.1
+    t0 = 0.112
+    t1 = 0.212
+    if not recompute and spec_path.is_file() and rate_path.is_file():
+        tmax = t0 + t1
+        trials = 20
+    else:
+        print('    compute spike response')
     spikes = punit_spikes(cell, alpha, beatf1, beatf2, tmax, trials)
-    plot_am(axs, s, alpha, beatf1, beatf2, t1)
-    plot_raster(axr, s, spikes, t1)
-    plot_psd(axp, s, path, alpha, spikes, nfft, dt, beatf1, beatf2)
-
-
-def peak_ampl(freqs, psd, f):
-    df = 2
-    psd_snippet = psd[(freqs > f - df) & (freqs < f + df)]
-    return np.max(psd_snippet)
+    plot_am(axs, s, alpha, beatf1, beatf2, t0, t1)
+    plot_raster(axr, s, spikes, t0, t1)
+    plot_rate(axf, s, rate_path, spikes, t0, t1, sigma)
+    plot_psd(axp, s, spec_path, alpha, spikes, nfft, dt, beatf1, beatf2, eodf)
 
 
 def amplitude(power):
@@ -243,16 +351,15 @@ def amplitude_squarefit(contrast, power, max_contrast):
 
 
 def plot_peaks(ax, s, alphas, contrasts, powerf1, powerf2, powerfsum,
-               powerfdiff):
-    cmax = 10
+               powerfdiff, cmax=8):
     contrasts *= 100
-    ax.plot(contrasts, amplitude_linearfit(contrasts, powerf1, 4),
+    ax.plot(contrasts, amplitude_linearfit(contrasts, powerf1, cthresh1),
             **s.lsF01m)
-    ax.plot(contrasts, amplitude_linearfit(contrasts, powerf2, 2),
+    ax.plot(contrasts, amplitude_linearfit(contrasts, powerf2, cthresh1),
             **s.lsF02m)
-    ax.plot(contrasts, amplitude_squarefit(contrasts, powerfsum, 4),
+    ax.plot(contrasts, amplitude_squarefit(contrasts, powerfsum, cthresh2),
             **s.lsF012m)
-    ax.plot(contrasts, amplitude_squarefit(contrasts, powerfdiff, 4),
+    ax.plot(contrasts, amplitude_squarefit(contrasts, powerfdiff, cthresh2),
             **s.lsF01_2m)
     ax.plot(contrasts, amplitude(powerf1), **s.lsF01)
     ax.plot(contrasts, amplitude(powerf2), **s.lsF02)
@@ -261,28 +368,29 @@ def plot_peaks(ax, s, alphas, contrasts, powerf1, powerf2, powerfsum,
             clip_on=False, **s.lsF012)
     ax.plot(contrasts[mask], amplitude(powerfdiff)[mask],
             clip_on=False, **s.lsF01_2)
-    ymax = 60
     for alpha, tag in zip(alphas, ['A', 'B', 'C', 'D']):
-        ax.plot(100*alpha, ymax*0.95, 'vk', ms=4, clip_on=False)
-        ax.text(100*alpha, ymax, tag, ha='center')
+        ax.plot(100*alpha, 1.05*amax, 'vk', ms=4, clip_on=False)
+        ax.text(100*alpha, 1.1*amax, tag, ha='center')
         #ax.axvline(contrast, **s.lsGrid)
         #ax.text(contrast, 630, tag, ha='center')
-    cthresh1 = 1.2
-    cthresh2 = 3.5
     print(f'Linear            regime ends at a contrast of about {cthresh1:4.1f}%')
     print(f'Weakly non-linear regime ends at a contrast of about {cthresh2:4.1f}%')
     ax.axvline(cthresh1, **s.lsLine)
     ax.axvline(cthresh2, **s.lsLine)
-    yoffs = 35
+    yoffs = 35 if amax == 60 else 31
     ax.text(cthresh1/2, yoffs, 'linear\nregime',
             ha='center', va='center')
-    ax.text((cthresh1 + cthresh2)/2, yoffs, 'weakly\nnonlinear\nregime',
-            ha='center', va='center')
-    ax.text(5.5, yoffs, 'strongly\nnonlinear\nregime',
+    ax.text(2.2, yoffs, 'weakly\nnonlinear\nregime',
             ha='center', va='center')
+    if amax == 60:
+        ax.text(5.5, yoffs, 'strongly\nnonlinear\nregime',
+                ha='center', va='center')
+    else:
+        ax.text(5.5, 6, 'strongly\nnonlinear\nregime',
+                ha='center', va='bottom')
     ax.set_xlim(0, cmax)
-    ax.set_ylim(0, ymax)
-    ax.set_xticks_delta(2)
+    ax.set_ylim(0, amax)
+    ax.set_xticks_delta(1)
     ax.set_yticks_delta(20)
     ax.set_xlabel('Contrast', r'\%')
     ax.set_ylabel('Amplitude', 'Hz')
@@ -292,35 +400,38 @@ if __name__ == '__main__':
     ratebase, cvbase, beatf1, beatf2, \
         contrasts, powerf1, powerf2, powerfsum, powerfdiff = \
             load_data(sims_path / f'{model_cell}-contrastpeaks.npz')
-    alphas = [0.002, 0.01, 0.03, 0.06]
 
     parameters = load_models(data_path / 'punitmodels.csv')
     cell = cell_parameters(parameters, model_cell)
-    nfft = 2**18
+    eodf = cell['EODf']
     
     print(f'Loaded data for cell {model_cell}: ')
     print(f'    baseline rate = {ratebase:.0f}Hz, CV = {cvbase:.2f}')
-    print(f'    f1 = {beatf1:.0f}Hz, f2 = {beatf2:.0f}Hz')
+    print(f'     f1 = {eodf + beatf1:.0f}Hz, f2  = {eodf + beatf2:.0f}Hz, eodf = {eodf:.0f}Hz')
+    print(f'    Df1 = {beatf1:.0f}Hz, Df2 = {beatf2:.0f}Hz')
     print()
 
     s = plot_style()
-    fig, (axes, axa) = plt.subplots(2, 1, height_ratios=[4, 3],
-                                    cmsize=(s.plot_width, 0.6*s.plot_width))
+    fig, (axes, axa) = plt.subplots(2, 1, height_ratios=[5, 3],
+                                    cmsize=(s.plot_width, 0.7*s.plot_width))
     fig.subplots_adjust(leftm=8, rightm=2, topm=2, bottomm=3.5, hspace=0.6)
-    axe = axes.subplots(3, 4, wspace=0.4, hspace=0.2,
-                        height_ratios=[1, 2, 3])
+    axe = axes.subplots(4, 4, wspace=0.4, hspace=0.3,
+                        height_ratios=[1, 2, 2, 0.2, 3])
     fig.show_spines('lb')
     
     # example power spectra:
     for c, alpha in enumerate(alphas):
-        path = sims_path / f'{model_cell}-contrastspectrum-{1000*alpha:03.0f}.npz'
-        plot_example(axe[0, c], axe[1, c], axe[2, c], s, path,
-                     cell, alpha, beatf1, beatf2, nfft, 100)
-    axe[1, 0].xscalebar(1, -0.1, 20, 'ms', ha='right')
-    axe[2, 0].legend(loc='center left', bbox_to_anchor=(0, -0.8),
-                     ncol=5, columnspacing=2)
+        path = sims_path / f'{model_cell}'
+        print(f'Example response for contrast {100*alpha:4.1f}%:')
+        plot_example(axe[0, c], axe[1, c], axe[2, c], axe[3, c], s, path,
+                     cell, alpha, beatf1, beatf2, eodf, nfft, trials)
+        print()
+    axe[2, 0].xscalebar(1, -0.1, 20, 'ms', ha='right')
+    axe[3, 3].legend(loc='center right', bbox_to_anchor=(1.05, -0.9),
+                     ncol=11, columnspacing=0.6, handletextpad=0)
     fig.common_yspines(axe[0, :])
     fig.common_yticks(axe[2, :])
+    fig.common_yticks(axe[3, :])
     fig.tag(axe[0, :], xoffs=-3, yoffs=1.6)
     
     # contrast dependence:

twobeats.py

@@ -0,0 +1,341 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pathlib import Path
+from scipy.stats import norm
+from scipy.optimize import curve_fit
+from spectral import rate
+from plotstyle import plot_style
+
+
+cell = '2021-08-03-ac-invivo-1'
+
+data_path = Path('data')
+
+sigma = 0.001
+
+
+def load_spikes(cell_path, f1=797, f2=631):
+    load = False
+    spikes = []
+    index = 0
+    with open(cell_path / 'threefish-spikes.dat') as sf:
+        for line in sf:
+            if load:
+                if ' before:' in line:
+                    t0 = 0.001*float(line.split(':')[1].strip().replace('ms', ''))
+                elif ' duration1 ' in line:
+                    t1 = 0.001*float(line.split(':')[1].strip().replace('ms', ''))
+                elif ' duration2 ' in line:
+                    t2 = 0.001*float(line.split(':')[1].strip().replace('ms', ''))
+                elif ' duration12 ' in line:
+                    t12 = 0.001*float(line.split(':')[1].strip().replace('ms', ''))
+                elif line.startswith('# index '):
+                    if len(spikes) > 0:
+                        spikes[-1] = np.array(spikes[-1])
+                    return spikes, eodf, df1, df2, t0, t1, t2, t12, index
+                elif line.startswith('# trial:'):
+                    if len(spikes) > 0:
+                        spikes[-1] = np.array(spikes[-1])
+                    spikes.append([])
+                elif len(line.strip()) > 0 and line[0] != '#':
+                    t = 0.001*float(line.strip())
+                    spikes[-1].append(t)
+            elif line.startswith('# index '):
+                index += 1
+            elif line.startswith('# EOD rate '):
+                eodf = float(line.split(':')[1].strip().replace('Hz', ''))
+            elif line.startswith('# Deltaf1 '):
+                df1 = float(line.split(':')[1].strip().replace('Hz', ''))
+            elif line.startswith('# Deltaf2 '):
+                df2 = float(line.split(':')[1].strip().replace('Hz', ''))
+                if abs(eodf + df1 - f1) < 1 and abs(eodf + df2 - f2) < 1:
+                    #print(f'EODf={eodf:6.1f}Hz, Df1={df1:6.1f}Hz, Df2={df2:6.1f}Hz, EODf1={eodf + df1:6.1f}Hz, EODf2={eodf + df2:6.1f}Hz')
+                    load = True
+    print(f'no spikes found for EODf1={f1:.1f}Hz and EODf2={f2:.1f}Hz')
+
+
+def load_am(cell_path, inx):
+    load = False
+    ams = []
+    index = 0
+    with open(cell_path / 'threefish-ams.dat') as sf:
+        for line in sf:
+            if load:
+                if line.startswith('# index '):
+                    if len(ams) > 0:
+                        ams[-1] = np.array(ams[-1])
+                    return ams
+                elif line.startswith('# EOD rate '):
+                    print(f'    EODf = {line.split(':')[1].strip()}')
+                elif line.startswith('# Deltaf1 '):
+                    print(f'    Df1  = {line.split(':')[1].strip()}')
+                elif line.startswith('# Deltaf2 '):
+                    print(f'    DF2  = {line.split(':')[1].strip()}')
+                elif line.startswith('# trial:'):
+                    if len(ams) > 0:
+                        ams[-1] = np.array(ams[-1])
+                    ams.append([])
+                elif len(line.strip()) > 0 and line[0] != '#':
+                    time, am = line.split()
+                    t = 0.001*float(time.strip())
+                    a = float(am.strip())
+                    ams[-1].append((t, a))
+            elif line.startswith('# index '):
+                index += 1
+                if inx == index:
+                    load = True
+    print(f'no AM found at index {inx}')
+
+
+def cosine(x, a, f, p, c):
+    return a*np.cos(2*np.pi*f*x + p) + c
+
+
+def two_cosine(x, a1, f1, p1, a2, f2, p2, c):
+    return a1*np.cos(2*np.pi*f1*x + p1) + a2*np.cos(2*np.pi*f2*x + p2) + c
+
+
+def am_phases(ams, eodf, df1, df2, t1, t2, t12):
+    twins = (t1, t2, t12)
+    dfs = ((df1,), (df2,), (df1, df2))
+    phases = np.zeros((len(ams), len(dfs) + 1))
+    for k in range(len(ams)):
+        t0 = 0
+        time = ams[k][:, 0]
+        am = ams[k][:, 1]
+        for i in range(len(twins)):
+            tw = twins[0]
+            t1 = t0 + tw
+            mask = (time >= t0) & (time <= t1)
+            tam = time[mask] - t0
+            aam = am[mask]
+            a = 0.5*(np.max(aam) - np.min(aam))
+            c = np.mean(aam)
+            tt = np.linspace(0, tw, 1000)
+            if len(dfs[i]) == 2:
+                popt = [a/2, dfs[i][0], 0, a/2, dfs[i][1], 0, c]
+                popt, _ = curve_fit(two_cosine, tam, aam, popt)
+                aa = two_cosine(tt, *popt)
+                phases[k, i] = popt[2] if popt[0] > 0 else popt[2] + np.pi
+                phases[k, i + 1] = popt[5] if popt[3] > 0 else popt[5] + np.pi
+            else:
+                popt = [a, dfs[i][0], 0, c]
+                popt, _ = curve_fit(cosine, tam, aam, popt)
+                aa = cosine(tt, *popt)
+                phases[k, i] = popt[2] if popt[0] > 0 else popt[2] + np.pi
+            t0 = t1
+    return phases
+
+
+def align_spikes(spikes, freqs, phases):
+    f1, f2 = freqs
+    if f1 is None and f2 is None:
+        return spikes
+    p1 = phases[0]
+    p2 = phases[1]
+    if f2 is None:
+        df = f1
+        p = p1
+    else:
+        df = f2
+        p = p2
+    for i in range(len(spikes)):
+        spikes[i] += p[i]/2/np.pi/df
+    return spikes
+
+
+def baseline_rate(spikes, t0, t1):
+    rates = []
+    for times in spikes:
+        c = np.sum((times > t0) & (times < t1))
+        rates.append(c/(t1 - t0))
+    return np.mean(rates)
+
+    
+def power_spectrum(spikes, tmax, dt=0.0005, nfft=512, p_ref=4000):
+    time = np.arange(0, tmax, dt)
+    if nfft > len(time):
+        print('nfft too large:', nfft, len(time))
+    freqs = np.fft.fftfreq(nfft, dt)
+    freqs = np.fft.fftshift(freqs)
+    segments = range(0, len(time) - nfft, nfft)
+    p_rr = np.zeros(len(freqs))
+    n = 0
+    for i, spiket in enumerate(spikes):
+        b, _ = np.histogram(spiket, time)
+        b = b / dt
+        for j, k in enumerate(segments):
+            fourier_r = np.fft.fft(b[k:k + nfft] - np.mean(b), n=nfft)
+            fourier_r = np.fft.fftshift(fourier_r)
+            p_rr += np.abs(fourier_r*np.conj(fourier_r))
+            n += 1
+    mask = freqs >= 0.0
+    freqs = freqs[mask]
+    scale = dt/nfft/n
+    p_rr = p_rr[mask]*scale
+    power = 10*np.log10(p_rr/p_ref)
+    return freqs, power
+
+
+def plot_symbols(ax, s, freqs):
+    f1, f2 = freqs
+    ax.show_spines('')
+    ax.add_artist(plt.Rectangle((-1, -0.5), 2, 1, color=s.colors['black']))
+    ax.harrow(1.6, 0, 1.3, **s.asLine)
+    ax.set_xlim(-6, 14)
+    ax.set_ylim(-1, 1)
+    if f1 is None and f2 is None:
+        ax.text(3.5, 0, '$r$', va='center')
+    else:
+        ax.harrow(-2.8, 0, 1.3, **s.asLine)
+        if f2 is None:
+            ax.text(-3.2, 0, '$s_1(t)$', ha='right', va='center')
+            ax.text(3.3, 0, '$r + r_1(t)$', va='center')
+        elif f1 is None:
+            ax.text(-3.2, 0, '$s_2(t)$', ha='right', va='center')
+            ax.text(3.3, 0, '$r + r_2(t)$', va='center')
+        else:
+            ax.text(-3.2, 0, '$s_1(t) + s_2(t)$', ha='right', va='center')
+            ax.text(3.3, 0, '$\\ne r + r_1(t) + r_2(t)$', va='center')
+
+
+def plot_stimulus(ax, s, tmax, eodf, freqs, c=0.1):
+    time = np.arange(0, tmax, 0.0001)
+    eod = np.cos(2*np.pi*eodf*time)
+    am = np.ones(len(time))
+    ams = {}
+    f1, f2 = freqs
+    label = '$f_{EOD}$'
+    dp = np.pi
+    if f1 is not None:
+        eod += c*np.cos(2*np.pi*(eodf + f1)*time + dp)
+        am += c*np.cos(2*np.pi*f1*time + dp)
+        ams = s.lsF02
+        label += r' \& $f_1$'
+    if f2 is not None:
+        eod += c*np.cos(2*np.pi*(eodf + f2)*time + dp)
+        am += c*np.cos(2*np.pi*f2*time + dp)
+        ams = s.lsF01
+        label += r' \& $f_2$'
+    if f1 is not None and f2 is not None:
+        ams = s.lsF01_2
+    ax.show_spines('')
+    ax.plot(1000*time, eod, **s.lsEOD)
+    if len(ams) > 0:
+        ax.plot(1000*time, am + 0.1, **ams)
+    ax.set_xlim(0, 1000*tmax)
+    ax.set_ylim(-1.02 - 2*c, 1.1 + 2*c)
+    ax.text(0.5, 1.2, label, ha='center', transform=ax.transAxes)
+
+
+def plot_raster(ax, s, spikes, tmin, tmax):
+    spikes_ms = [1000*(s[(s > tmin) & (s < tmax)] - tmin) for s in spikes]
+    ax.show_spines('')
+    ax.eventplot(spikes_ms, linelengths=0.8, **s.lsRaster)
+    ax.set_xlim(0, 1000*(tmax - tmin))
+
+    
+def plot_rate(ax, s, spikes, tmin, tmax, sigma=0.002):
+    time = np.arange(0, tmin + tmax, sigma/4)
+    r, rsd = rate(time, spikes, sigma)
+    mask = (time >= tmin) & (time <= tmax)
+    time = time[mask] - tmin
+    r = r[mask]
+    ax.show_spines('l')
+    ax.plot(1000*time, r, clip_on=False, **s.lsRate)
+    ax.set_xlim(0, 1000*(tmax - tmin))
+    ax.set_ylim(0, 550)
+    ax.set_ylabel('Rate', 'Hz')
+    ax.set_yticks_delta(200)
+
+    
+def plot_psd(ax, s, freqs, power, fmax, dt=0.0005, nfft=512):
+    # plot:
+    mask = freqs <= fmax
+    freqs = freqs[mask]
+    power = power[mask]
+    ax.show_spines('lb')
+    ax.plot(freqs, power, **s.lsPower)
+    ax.set_xlim(0, fmax)
+    ax.set_ylim(-20, 0)
+    ax.set_xlabel('Frequency', 'Hz')
+    ax.set_ylabel('Power', 'dB')
+    ax.set_yticks_delta(10)
+
+
+def mark_freq(ax, freqs, power, f, label, style, xoffs=10, yoffs=0, toffs=0, angle=0):
+    i = np.argmin(np.abs(freqs - abs(f)))
+    p = power[i]
+    f = freqs[i]
+    ax.plot(f, p + 1 + yoffs, clip_on=False, **style)
+    if label:
+        yoffs += 3 + toffs
+        if angle > 0:
+            yoffs -= 1
+        ax.text(f - xoffs, p + yoffs, label, color=style['color'], rotation=angle)
+        
+    
+if __name__ == '__main__':
+
+    spikes, eodf, df1, df2, t0, t1, t2, t12, index = load_spikes(data_path / cell)
+    print(f'Loaded spike data for cell {cell} @ index {index}:')
+    print(f'    EODf = {eodf:.1f}Hz')
+    print(f'    Df1  = {df1:.1f}Hz')
+    print(f'    Df2  = {df2:.1f}Hz')
+    print(f'    {len(spikes)} trials')
+
+    print(f'Load AMs for cell {cell} @ index {index}:')
+    ams = load_am(data_path / cell, index)
+
+    phases = am_phases(ams, eodf, df1, df2, t1, t2, t12)
+    
+    s = plot_style()
+    fig, axs = plt.subplots(5, 4, cmsize=(s.plot_width, 0.6*s.plot_width),
+                            height_ratios=[1, 0, 2, 1.3, 3, 0.7, 5])
+    fig.subplots_adjust(leftm=7.5, rightm=4, topm=1.5, bottomm=4,
+                        wspace=0.4, hspace=0.4)
+    fmax = 250
+    tmin = 0.106
+    tmax = 0.206
+    twins = [[-t0, 0], [t1, t1 + t2], [0, t1], [t1 + t2, t1 + t2 + t12]]
+    stim_freqs = [[None, None], [df2, None], [None, df1], [df1, df2]]
+    stim_phases = [[None, None], [phases[:, 1], None], [None, phases[:, 0]], [phases[:, 2], phases[:, 3]]]
+    base_rate = baseline_rate(spikes, *twins[0])
+    print(f'Baseline firing rate: {base_rate:.1f}Hz')
+    powers = []
+    for i in range(axs.shape[1]):
+        tstart, tend = twins[i]
+        plot_symbols(axs[0, i], s, stim_freqs[i])
+        plot_stimulus(axs[1, i], s, tmax - tmin, eodf, stim_freqs[i])
+        sub_spikes = [times[(times >= tstart) & (times <= tend)] - tstart for times in spikes]
+        freqs, power = power_spectrum(sub_spikes, tend - tstart)
+        powers.append(power)
+        plot_psd(axs[4, i], s, freqs, power, fmax)
+        sub_spikes = align_spikes(sub_spikes, stim_freqs[i], stim_phases[i])
+        plot_raster(axs[2, i], s, sub_spikes, tmin, tmax)
+        plot_rate(axs[3, i], s, sub_spikes, tmin, tmax, sigma)
+    mark_freq(axs[4, 0], freqs, powers[0], base_rate, f'$r={base_rate:.0f}$\\,Hz', s.psF0, 30)
+    mark_freq(axs[4, 1], freqs, powers[1], df2, f'$\\Delta f_1=f_1 - f_{{EOD}}={abs(df2):.0f}$\\,Hz', s.psF02)
+    mark_freq(axs[4, 1], freqs, powers[1], 2*df2, f'$2\\Delta f_1={abs(2*df2):.0f}$\\,Hz', s.psF02)
+    mark_freq(axs[4, 2], freqs, powers[2], df1, '', s.psF0)
+    mark_freq(axs[4, 2], freqs, powers[2], df1, f'$\\Delta f_2=f_2 - f_{{EOD}}={abs(df1):.0f}$\\,Hz',
+              s.psF01, 130, 1.5)
+    mark_freq(axs[4, 3], freqs, powers[3], df2, '', s.psF02)
+    mark_freq(axs[4, 3], freqs, powers[3], 2*df2, '', s.psF02)
+    mark_freq(axs[4, 3], freqs, powers[3], df1, '', s.psF0)
+    mark_freq(axs[4, 3], freqs, powers[3], df1, '', s.psF01, 130, 1.5)
+    mark_freq(axs[4, 3], freqs, powers[3], abs(df1) + abs(df2) - 2,
+              f'$\\Delta f_2 + \\Delta f_1$\n$={abs(df1) + abs(df2):.0f}$\\,Hz', s.psF012, 0, toffs=-1, angle=45)
+    mark_freq(axs[4, 3], freqs, powers[3], abs(df1) - abs(df2),
+              f'$\\Delta f_2 - \\Delta f_1$\n$={abs(df1) - abs(df2):.0f}$\\,Hz', s.psF01_2, 65, toffs=-2, angle=45)
+    fig.common_yticks(axs[3, :])
+    fig.common_yticks(axs[4, :])
+    axs[3, 0].xscalebar(1, 0, 20, 'ms', ha='right')
+    #axs[3, 0].scalebars(-0.03, 0, 20, 500, 'ms', 'Hz')
+    #axs[4, 0].yscalebar(-0.03, 0.5, 10, 'dB', va='center')
+    #fig.tag(axs.T)
+    fig.tag(axs[0])
+    fig.savefig()
+    #plt.show()
+    print()