nonlinearbaseline2025/nonlinearbaseline.tex

915 lines
116 KiB
TeX

\documentclass[10pt,letterpaper]{article}
\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},
Jan Grewe\textsuperscript{1},
Jan Benda\textsuperscript{1, 4, $\dagger$}}
\date{\normalsize
\textsuperscript{1} Institute for Neurobiology, Eberhard Karls Universit\"at T\"ubingen, Germany\\
\textsuperscript{2} Bernstein Center for Computational Neuroscience Berlin, Berlin, Germany\\
\textsuperscript{3} Department of Physics, Humboldt University Berlin, Berlin, Germany\\
\textsuperscript{4} Bernstein Center for Computational Neuroscience Tübingen, Germany\\
\textsuperscript{$\dagger$} corresponding author: \url{jan.benda@uni-tuebingen.de}}
\newcommand{\runningtitle}{}
%%%%% overall style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% number the lines:
\newcommand{\setlinenumbers}{false}
% use double spacing:
\newcommand{\setdoublespacing}{false}
% show authors:
\newcommand{\showauthors}{true}
% show keywords:
\newcommand{\showkeywords}{true}
% show author contributions:
\newcommand{\showcontributions}{true}
% figures and captions at end:
\newcommand{\figsatend}{false}
% remove images from figures:
\newcommand{\nofigs}{false}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% line numbering, spacing, authors, keywords %%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{ifthen}
\usepackage[nointegrals]{wasysym}
\ifthenelse{\boolean{\setlinenumbers}}{%
\usepackage{lineno}\linenumbers}%
{\newenvironment{linenomath}{}{}}
\ifthenelse{\boolean{\setdoublespacing}}{%
\usepackage{setspace}\doublespacing}{}
\ifthenelse{\boolean{\showauthors}}{}{\author{}\date{}}
\usepackage[inline]{enumitem}
\newenvironment{keywords}%
{\paragraph{Keywords}\begin{itemize*}[label={},itemjoin={~$|$\ },afterlabel={}]}%
{\end{itemize*}\par}
\usepackage{comment}
\ifthenelse{\boolean{\showkeywords}}{}{\excludecomment{keywords}}
\newenvironment{contributions}{\section{Author contributions}}{\par}
\ifthenelse{\boolean{\showcontributions}}{}{\excludecomment{contributions}}
%%%%% page style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[left=25mm,right=25mm,top=25mm,bottom=25mm]{geometry}
\pagestyle{myheadings}
\markright{\runningtitle}
%%%%% language %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[english]{babel}
%%%%% fonts %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\usepackage{pslatex} % nice font for pdf file
%%%%% section style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[sf,bf,it,big,clearempty]{titlesec}
\usepackage{titling}
\renewcommand{\maketitlehooka}{\sffamily\bfseries}
\renewcommand{\maketitlehookb}{\rmfamily\mdseries}
\setcounter{secnumdepth}{-1}
%%%%% math %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsmath}
%\usepackage{lualatex-math}
\usepackage{iftex}
% LuaTeX-spezifischen Code nur ausführen, wenn LuaTeX verwendet wird
\ifluatex
\usepackage{lualatex-math}
\fi
%%%%% units %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[mediumspace,mediumqspace,Gray]{SIunits} % \ohm, \micro
%%%%% graphics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{xcolor}
%\pagecolor{white}
\usepackage{graphicx}
\usepackage{tabularx}
\newcommand{\trh}{\rule[-1ex]{0pt}{3.5ex}}
\newcommand{\trt}{\rule{0pt}{2.5ex}}
\newcolumntype{L}[1]{>{\raggedright\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{C}[1]{>{\centering\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{R}[1]{>{\raggedleft\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\graphicspath{{../susceptibility1}}
%%%%% figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% captions:
\ifthenelse{\boolean{\setdoublespacing}}{%
\usepackage[format=plain,singlelinecheck=off,labelfont=bf,font={small,sf,doublespacing}]{caption}}%
{\usepackage[format=plain,singlelinecheck=off,labelfont=bf,font={small,sf}]{caption}}
% put caption on separate float:
\newcommand{\breakfloat}{\end{figure}\begin{figure}[t]}
% references to panels of a figure within the caption:
\newcommand{\figitem}[2][]{\ifthenelse{\equal{#1}{}}{\textsf{\bfseries #2}}{\textsf{\bfseries #2}}$_{\sf #1}$}
% references to panels of a figure within the text:
\newcommand{\panel}[2][]{\ifthenelse{\equal{#1}{}}{\textsf{#2}}{\textsf{#2}$_{\sf #1}$}}
% references to figures:
\newcommand{\fref}[1]{\textup{\ref{#1}}}
\newcommand{\subfref}[2]{\textup{\ref{#1}}\,\panel{#2}}
% references to figures in normal text:
\newcommand{\fig}{fig.}
\newcommand{\Fig}{Fig.}
\newcommand{\figs}{figs.}
\newcommand{\Figs}{Figs.}
\newcommand{\figref}[1]{\fig~\fref{#1}}
\newcommand{\Figref}[1]{\Fig~\fref{#1}}
\newcommand{\figsref}[1]{\figs~\fref{#1}}
\newcommand{\Figsref}[1]{\Figs~\fref{#1}}
\newcommand{\subfigref}[2]{\fig~\subfref{#1}{#2}}
\newcommand{\Subfigref}[2]{\Fig~\subfref{#1}{#2}}
\newcommand{\subfigsref}[2]{\figs~\subfref{#1}{#2}}
\newcommand{\Subfigsref}[2]{\Figs~\subfref{#1}{#2}}
% references to figures within brackets:
\newcommand{\figb}{fig.}
\newcommand{\Figb}{Fig.}
\newcommand{\figsb}{figs.}
\newcommand{\Figsb}{Figs.}
\newcommand{\figrefb}[1]{\figb~\fref{#1}}
\newcommand{\Figrefb}[1]{\Figb~\fref{#1}}
\newcommand{\figsrefb}[1]{\figsb~\fref{#1}}
\newcommand{\Figsrefb}[1]{\Figsb~\fref{#1}}
\newcommand{\subfigrefb}[2]{\figb~\subfref{#1}{#2}}
\newcommand{\Subfigrefb}[2]{\Figb~\subfref{#1}{#2}}
\newcommand{\subfigsrefb}[2]{\figsb~\subfref{#1}{#2}}
\newcommand{\Subfigsrefb}[2]{\Figsb~\subfref{#1}{#2}}
% figure placement:
\ifthenelse{\boolean{\figsatend}}{%
\setcounter{topnumber}{10}%
\setcounter{totalnumber}{10}%
\renewcommand{\textfraction}{0.0}%
\renewcommand{\topfraction}{1.0}%
\renewcommand{\floatpagefraction}{0.0}%
\usepackage[nolists,markers,figuresonly]{endfloat}%
\ifthenelse{\boolean{\nofigs}}{\renewcommand{\efloatseparator}{\mbox{}}}{}%
\newcommand{\figurecaptions}{\pagestyle{empty}\clearpage\processdelayedfloats\renewcommand{\figureplace}{}}}%
{\setcounter{topnumber}{1}%
\setcounter{totalnumber}{2}%
\renewcommand{\textfraction}{0.2}%
\renewcommand{\topfraction}{0.8}%
\renewcommand{\floatpagefraction}{0.7}%
\newcommand{\figurecaptions}{}%
}
\newcommand*{\captionc}[1]{\captionof{figure}{#1}}
% no images in figures:
\ifthenelse{\boolean{\nofigs}}{%
\newcommand{\showfigure}[1]{}}%
{\newcommand{\showfigure}[1]{#1}}
% bibliography ----------------------------------
\usepackage[round,colon]{natbib}
\renewcommand{\bibsection}{\section{References}}
\setlength{\bibsep}{0pt}
\setlength{\bibhang}{1.5em}
\bibliographystyle{jneurosci}
\usepackage[breaklinks=true,colorlinks=true,citecolor=blue!30!black,urlcolor=blue!30!black,linkcolor=blue!30!black]{hyperref}
%%%%% tables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% references to tables:
\newcommand{\tref}[1]{\textup{\ref{#1}}}
% references to tables in normal text:
\newcommand{\tab}{Tab.}
\newcommand{\Tab}{Table}
\newcommand{\tabs}{Tabs.}
\newcommand{\Tabs}{Tables}
\newcommand{\tabref}[1]{\tab~\tref{#1}}
\newcommand{\Tabref}[1]{\Tab~\tref{#1}}
\newcommand{\tabsref}[1]{\tabs~\tref{#1}}
\newcommand{\Tabsref}[1]{\Tabs~\tref{#1}}
% references to tables within the bracketed text:
\newcommand{\tabb}{Tab.}
\newcommand{\tabsb}{Tab.}
\newcommand{\tabrefb}[1]{\tabb~\tref{#1}}
\newcommand{\tabsrefb}[1]{\tabsb~\tref{#1}}
%%%%% equation references %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\eqref}[1]{(\ref{#1})}
\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][]{}
\newcommand{\notejb}[1]{\note[JB]{#1}}
\newcommand{\notejg}[1]{\note[JG]{#1}}
\newcommand{\notesr}[1]{\note[SR]{#1}}
\newcommand{\noteab}[1]{\note[AB]{#1}}
\newcommand{\notebl}[1]{\note[BL]{#1}}
\newcommand{\notems}[1]{\note[MS]{#1}}
%%%%% new commands %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\fstim}{\ensuremath{f_{\rm{stim}}}}
\newcommand{\feod}{\ensuremath{f_{\rm{EOD}}}}
\newcommand{\feodsolid}{\ensuremath{f_{\rm{EOD}}}}
\newcommand{\feodhalf}{\ensuremath{f_{\rm{EOD}}/2}}
\newcommand{\fbase}{\ensuremath{f_{\rm{base}}}}
\newcommand{\fbasesolid}{\ensuremath{f_{\rm{base}}}}
\newcommand{\fstimintro}{\ensuremath{\rm{EOD}_{2}}}
\newcommand{\feodintro}{\ensuremath{\rm{EOD}_{1}}}
\newcommand{\ffstimintro}{\ensuremath{f_{2}}}
\newcommand{\ffeodintro}{\ensuremath{f_{1}}}
\newcommand{\carrierinput}{\ensuremath{y(t)}}
\newcommand{\baseval}{134}
\newcommand{\bone}{\ensuremath{\Delta f_{1}}}
\newcommand{\btwo}{\ensuremath{\Delta f_{2}}}
\newcommand{\fone}{$f_{1}$}
\newcommand{\ftwo}{$f_{2}$}
\newcommand{\ff}{$f_{1}$--$f_{2}$}
\newcommand{\auc}{\rm{AUC}}
\newcommand{\n}{\ensuremath{N}}
\newcommand{\bvary}{\ensuremath{\Delta f_{1}}}
\newcommand{\avary}{\ensuremath{A(\Delta f_{1})}}
\newcommand{\afvary}{\ensuremath{A(f_{1})}}
\newcommand{\avaryc}{\ensuremath{A(c_{1})}}
\newcommand{\cvary}{\ensuremath{c_{1}}}
\newcommand{\fvary}{\ensuremath{f_{1}}}
\newcommand{\bstable}{\ensuremath{\Delta f_{2}}}
\newcommand{\astable}{\ensuremath{A(\Delta f_{2})}}
\newcommand{\astablec}{\ensuremath{A(c_{2})}}
\newcommand{\cstable}{\ensuremath{c_{2}}}
\newcommand{\fstable}{\ensuremath{f_{2}}}
\newcommand{\aeod}{\ensuremath{A(f_{\rm{EOD}})}}
\newcommand{\fbasecorrsolid}{\ensuremath{f_{\rm{BaseCorrected}}}}
\newcommand{\fbasecorr}{\ensuremath{f_{\rm{BaseCorrected}}}}
\newcommand{\ffall}{$f_{\rm{EOD}}$\&$f_{1}$\&$f_{2}$}
\newcommand{\ffvary}{$f_{\rm{EOD}}$\&$f_{1}$}%sum
\newcommand{\ffstable}{$f_{\rm{EOD}}$\&$f_{2}$}%sum
\newcommand{\colstableone}{blue}%sum
\newcommand{\colstabletwo}{cyan}%sum
\newcommand{\colvaryone}{brown}%sum
\newcommand{\colvarytwo}{red}%sum
\newcommand{\colstableoneb}{Blue}%sum
\newcommand{\colstabletwob}{Cyan}%sum
\newcommand{\colvaryoneb}{Brown}%sum Brown
\newcommand{\colvarytwob}{Red}%sum
\newcommand{\resp}{$A(f)$}%sum %\Delta
\newcommand{\respb}{$A(\Delta f)$}%sum %\Delta
\newcommand{\rocf}{98}%sum %\Delta
\newcommand{\roci}{36}%sum %\Delta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Beat combinations
\newcommand{\boneabs}{\ensuremath{|\Delta f_{1}|}}%sum
\newcommand{\btwoabs}{\ensuremath{|\Delta f_{2}|}}%sum
%\newcommand{\bsum}{\ensuremath{|\boneabs{} + \btwoabs{}|}}%sum
\newcommand{\bsum}{\ensuremath{\boneabs{} + \btwoabs{}}}%sum
\newcommand{\bdiff}{\ensuremath{|\boneabs{} - \btwoabs{}|}}%diff of both beat frequencies
\newcommand{\signalnoise}{$s_\xi(t)$}
\newcommand{\bsumb}{$\bsum{}=\fbase{}$}
\newcommand{\btwob}{$\Delta f_{2}=\fbase{}$}
\newcommand{\boneb}{$\Delta f_{1}=\fbase{}$}
\newcommand{\bsumbtwo}{$\bsum{}=2 \fbase{}$}
\newcommand{\bsumbc}{$\bsum{}=\fbasecorr{}$}
\newcommand{\bsume}{$\bsum{}=\feod{}$}
\newcommand{\bsumehalf}{$\bsum{}=\feod{}/2$}
\newcommand{\bdiffb}{$\bdiff{}=\fbase{}$}%diff of both beat frequencies
\newcommand{\bdiffbc}{$\bdiff{}=\fbasecorr{}$}%diff of both beat frequencies
\newcommand{\bdiffe}{$\bdiff{}=\feod{}$}%diff of both
\newcommand{\bdiffehalf}{$\bdiff{}=\feod{}/2$}%diff of both beat frequencies
%%%%%%%%%%%%%%%%%%%%%%% Frequency combinations
\newcommand{\fsum}{\ensuremath{f_{1} + f_{2}}}%sum
\newcommand{\fdiff}{\ensuremath{|f_{1}-f_{2}|}}%diff of
\newcommand{\fsumb}{$\fsum=\fbase{}$}%su\right m
\newcommand{\ftwob}{$f_{2}=\fbase{}$}%sum
\newcommand{\foneb}{$f_{1}=\fbase{}$}%sum
\newcommand{\ftwobc}{$f_{2}=\fbasecorr{}$}%sum
\newcommand{\fonebc}{$f_{1}=\fbasecorr{}$}%sum
\newcommand{\fsumbtwo}{$\fsum{}=2 \fbase{}$}%sum
\newcommand{\fsumbc}{$\fsum{}=\fbasecorr{}$}%sum
\newcommand{\fsume}{$\fsum{}=\feod{}$}%sum
\newcommand{\fsumehalf}{$\fsum{}=\feod{}/2$}%sum
\newcommand{\fdiffb}{$\fdiff{}=\fbase{}$}%diff of both beat frequencies
\newcommand{\fdiffbc}{$\fdiff{}=\fbasecorr{}$}%diff of both beat frequencies
\newcommand{\fdiffe}{$\fdiff{}=\feod{}$}%diff of both
\newcommand{\fdiffehalf}{$\fdiff{}=\feod{}/2$}%diff of both
\newcommand{\fctwo}{\ensuremath{f_{\rm{Female}}}}%sum
\newcommand{\fcone}{\ensuremath{f_{\rm{Intruder}}}}%sum
\newcommand{\bctwo}{\ensuremath{\Delta \fctwo{}}}%sum
\newcommand{\bcone}{\ensuremath{\Delta \fcone{}}}%sum
\newcommand{\bcsum}{\ensuremath{|\bctwo{}| + |\bcone{}|}}%sum
\newcommand{\abcsumabbr}{\ensuremath{A(\Delta f_{\rm{Sum}})}}%sum
\newcommand{\abcsum}{\ensuremath{A(\bcsum{})}}%sum
\newcommand{\abcone}{\ensuremath{A(\bcone{})}}%sum
\newcommand{\bctwoshort}{\ensuremath{\Delta f_{F}}}%sum
\newcommand{\bconeshort}{\ensuremath{\Delta f_{I}}}%sum
\newcommand{\bcsumshort}{\ensuremath{|\bctwoshort{} + \bconeshort{}|}}%sum
\newcommand{\abcsumshort}{\ensuremath{A(\bcsumshort{})}}%sum
\newcommand{\abconeshort}{\ensuremath{A(|\bconeshort{}|)}}%sum
\newcommand{\abcsumb}{\ensuremath{A(\bcsum{})=\fbasesolid{}}}%sum
\newcommand{\bcdiff}{\ensuremath{||\bctwo{}| - |\bcone{}||}}%diff of both beat frequencies
\newcommand{\bcsumb}{\ensuremath{\bcsum{} =\fbasesolid{}}}%su\right m
\newcommand{\bcsumbn}{\ensuremath{\bcsum{} \neq \fbasesolid{}}}%su\right m
\newcommand{\bctwob}{\ensuremath{\bctwo{} =\fbasesolid{}}}%sum
\newcommand{\bconeb}{\ensuremath{\bcone{} =\fbasesolid{}}}%sum
\newcommand{\bcsumbtwo}{\ensuremath{\bcsum{}=2 \fbasesolid{}}}%sum
\newcommand{\bcsumbc}{\ensuremath{\bcsum{}=\fbasecorr{}}}%sum
\newcommand{\bcsume}{\ensuremath{\bcsum{}=f_{\rm{EOD}}}}%sum
\newcommand{\bcsumehalf}{\ensuremath{\bcsum{}=f_{\rm{EOD}}/2}}%sum
\newcommand{\bcdiffb}{\ensuremath{\bcdiff{}=\fbasesolid{}}}%diff of both beat frequencies
\newcommand{\bcdiffbc}{\ensuremath{\bcdiff{}=\fbasecorr{}}}%diff of both beat frequencies
\newcommand{\bcdiffe}{\ensuremath{\bcdif{}f=f_{\rm{EOD}}}}%diff of both
\newcommand{\bcdiffehalf}{\ensuremath{\bcdiff{}=f_{\rm{EOD}}/2}}%diff of both
\newcommand{\burstcorr}{\ensuremath{{Corrected}}}
\newcommand{\cvbasecorr}{CV\ensuremath{_{BaseCorrected}}}
\newcommand{\cv}{CV\ensuremath{_{Base}}}%\cvbasecorr{}
\newcommand{\nli}{PNL\ensuremath{(\fbase{})}}%Fr$_{Burst}$
\newcommand{\nlicorr}{PNL\ensuremath{(\fbasecorr{}})}%Fr$_{Burst}$
\newcommand{\suscept}{$|\chi_{2}|$}
\newcommand{\susceptf}{$|\chi_{2}|(f_1, f_2)$}
\newcommand{\frcolor}{pink lines}
\newcommand{\rec}{\ensuremath{\rm{R}}}%{\ensuremath{con_{R}}}
\newcommand{\rif}{\ensuremath{\rm{RIF}}}
\newcommand{\ri}{\ensuremath{\rm{RI}}}%sum
\newcommand{\rf}{\ensuremath{\rm{RF}}}%sum
\newcommand{\withfemale}{\rm{ROC}\ensuremath{\rm{_{Female}}}}%sum \textit{
\newcommand{\wofemale}{\rm{ROC}\ensuremath{\rm{_{NoFemale}}}}%sum
\newcommand{\dwithfemale}{\rm{\ensuremath{\auc_{Female}}}}%sum CV\ensuremath{_{BaseCorrected}}
\newcommand{\dwofemale}{\rm{\ensuremath{\auc_{NoFemale}}}}%sum
\newcommand{\fp}{\ensuremath{\rm{FP}}}%sum
\newcommand{\cd}{\ensuremath{\rm{CD}}}%sum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\maketitle
%\paragraph{Short title:}
%\paragraph{Corresponding author:}Jan Grewe, E-mail: jan.grewe@uni-tuebingen.de
%\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.
\begin{keywords}
\item second-order susceptibility
\item electric fish
\item nonlinear coding
\end{keywords}
% Please keep the abstract below 300 words
\section{Abstract}
Neuronal processing is inherently nonlinear --- spiking thresholds or rectification at synapses is essential to neuronal computations. Nevertheless, linear response theory has been instrumental in understanding, for example, the impact of noise or neuronal synchrony on signal transmission, or the emergence of oscillatory activity, but is valid only at low stimulus amplitudes or large levels of intrinsic noise. At higher signal-to-noise ratios, however, nonlinear response components become relevant. Theoretical results for leaky integrate-and-fire neurons in the weakly nonlinear regime suggest strong responses at the sum of two input frequencies if these frequencies or their sum match the neuron's baseline firing rate.
We here analyze nonlinear responses in two types of primary electroreceptor afferents, the P-units of the active and the ampullary cells of the passive electrosensory system of the wave-type electric fish \textit{Apteronotus leptorhynchus}. In our combined experimental and modeling approach we identify the predicted nonlinear responses in low-noise P-units and much stronger in all 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.
% Such nonlinear responses boost responses to weak sinusoidal stimuli and are therefore of immediate relevance for wave-type electric fish that are exposed to superpositions of many frequencies in social contexts.
% Please keep the Author Summary between 150 and 200 words
% Use the first person. PLOS ONE authors please skip this step.
% Author Summary not valid for PLOS ONE submissions.
%\section{Author summary}
%Weakly electric fish use their self-generated electric field to detect a wide range of behaviorally relevant stimuli. Intriguingly, they show detection performances of stimuli that are (i) extremely weak and (ii) occur in the background of strong foreground signals, reminiscent of what is often described as the cocktail party problem. Such performances are achieved by boosting the signal detection through nonlinear mechanisms. We here analyze nonlinear encoding in two different populations of primary electrosensory afferents of the weakly electric fish. We derive the rules under which nonlinear effects can be observed in both electrosensory subsystems. In a combined experimental and modeling approach we generalize the approach of nonlinear susceptibility to systems that respond to amplitude modulations of a carrier signal.
\note{Ad Aertsen paper suggestions:}
\url{https://brainworks.biologie.uni-freiburg.de/1981/journal%20papers/Aertsen-BiolCyb-1981b.pdf} \url{https://brainworks.biologie.uni-freiburg.de/1983/journal%20papers/eggermont-quartrevbiophys-1983.pdf} \url{https://brainworks.biologie.uni-freiburg.de/1981/journal%20papers/The%20Phonocrome%20(1981).pdf}
\section{Introduction}
\begin{figure*}[t]
\includegraphics[width=\columnwidth]{lifsuscept}
\caption{\label{fig:lifresponse} First- and second-order response functions of the leaky integrate-and-fire model. \figitem{A} Magnitude of the first-order (linear) response function $|\chi_1(f)|$, also known as the ``gain'' function, quantifies the response amplitude relative to the stimulus amplitude, both measured at the same stimulus frequency. \figitem{B} Magnitude of the second-order (non-linear) response function $|\chi_2(f_1, f_2)|$ quantifies the response at the sum of two stimulus frequencies. For linear systems, the second-order response function is zero, because linear systems do not create new frequencies and thus there is no response at the sum of the two frequencies. The plots show the analytical solutions from \citet{Lindner2001} and \citet{Voronenko2017} with $\mu = 1.1$ and $D = 0.001$. Note that the leaky integrate-and-fire model is formulated without dimensions, frequencies are given in multiples of the inverse membrane time constant.}
\end{figure*}
We like to think about signal encoding in terms of linear relations with unique mapping of a given input value to a certain output of the system under consideration. Indeed, such linear methods, for example the transfer function or first-oder susceptibility shown in figure~\ref{fig:lifresponse}, have been widely and successfully applied to describe and predict neuronal responses and are an invaluable tool to characterize neural systems \citep{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 spider mechanoreceptors \citep{French2001}, mammalian visual systems \citep{Victor1977, Schanze1997}, locking in chinchilla auditory nerve fibers \citep{Temchin2005}, and bursting responses in paddlefish \citep{Neiman2011} have been demonstrated.
Noise in nonlinear systems, however, linearizes the system's response properties \citep{Yu1989, Chialvo1997}. \notejg{kann den Gegensatz nicht erkennen..., kann der noise Einschub vielleicht woanders hin oder ist das eine Weiterfuehrung des noise Gedankens? Dann kann das vice versa weg.} Vice versa, in the limit to small stimuli, nonlinear systems can be well described by linear response theory \citep{Roddey2000, Doiron2004, Rocha2007, Sharafi2013}. With increasing stimulus amplitude, the contribution of the second-order kernel of the Volterra series becomes more relevant. For these weakly nonlinear responses analytical expressions for the second-order susceptibility have been derived for leaky-integrate-and-fire (LIF) \citep{Voronenko2017} and theta model neurons \citep{Franzen2023}. In the suprathreshold regime, 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:lifresponse}{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:lifresponse}{B}). In experimental data, such structures in the second-order susceptibility have not been reported yet. \notejg{... one of two stimulus frequencies... includes the case that one equals but the other is something completely different}
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}.
\section{Results}
\begin{figure*}[t]
\includegraphics[width=\columnwidth]{motivation.pdf}
\caption{\label{fig:motivation} Nonlinearity in an electrophysiologically recorded P-unit of \lepto{} in a three-fish setting (cell identifier ``2021-08-03-ac"). Receiver with EOD frequency $\feod{} =664$\,Hz encounters fish with EOD frequencies $f_{1}=631$\,Hz and $f_{2}=797$\,Hz. Both 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 \fbase{} 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 \feod{}, $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.
}
\end{figure*}
Theoretical work on leaky integrate-and-fire and conductance-based models suggests a distinct structure of the second-order response function for neurons with low levels of intrinsic noise driven in the supra-threshold regime with low stimulus amplitudes (\figrefb{fig:lifresponse}, \citealp{Voronenko2017}). Here, we explored a large set of recordings of P-units and ampullary cells of the active and passive electrosensory systems of the brown ghost knifefish \textit{Apteronotus leptorhynchus} to search for such weakly nonlinear responses in real neurons. Additional simulations of LIF-based models of P-unit spiking put the experimental findings into context. We start with demonstrating the basic concepts using example P-units and models. \notejg{relativ viel Wiederholung des ob gesagten, kuerzen?}
\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 \feod{}) and spontaneously fires action potentials at the baseline rate \fbase{}. Accordingly, the power spectrum of the baseline activity has a peak at \fbase{} (\subfigrefb{fig:motivation}{A}). Superposition of the receiver's EOD with an EOD of another fish with frequency $f_1$ results in a beat, a periodic amplitude modulation of the receiver's EOD. The frequency of the beat is given by the difference frequency $\Delta f_1 = f_1 - \feod$ between the two fish. P-units encode this beat in their firing rate \citep{Bastian1981a, Barayeu2023} and consequently 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 easily identified by the clipping of the P-unit's firing rate at zero. Pairing the fish with another fish at a higher beat frequency $\Delta f_2 = f_2 - \feod > \Delta f_1$ results in a weaker response with a single peak in the response power spectrum, suggesting a linear response (\subfigrefb{fig:motivation}{C}). The weaker response to this beat can be explained by the beat tuning of the cell \citep{Walz2014}. Note, $\Delta f_2$ has been deliberately chosen to match the recorded P-unit's baseline firing rate.
When stimulating with both foreign signals simultaneously, additional peaks appear in the response power spectrum at the sum $\Delta f_1 + \Delta f_2$ and the difference frequency $\Delta f_2 - \Delta f_1$ (\subfigrefb{fig:motivation}{D}). Thus, the cellular response is not equal to the sum of the responses to the two beats presented separately. These 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\notejg{should not occur?}.
\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'', table~\ref{modelparams}) was stimulated with two sinewaves 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 repsonses 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}.
At very low stimulus contrasts (in the example cell less than approximately 1.5\,\% relative to the receiver's EOD amplitude) the spectrum has small peaks only at the beat frequencies (\subfigref{fig:regimes}{A}, green and purple). The amplitudes of these peaks initially increase linearly with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines), an indication of the linear response at lowest stimulus amplitudes.
This linear regime is followed by the weakly nonlinear regime (in the example cell between approximately 1.5\,\% and 4\,\% stimulus contrasts). In addition to the peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:regimes}{B}, orange and red). The amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines). Note, that we have chosen $\Delta f_2$ to match the baseline firing rate $f_{base}$ of the neuron.
At higher stimulus amplitudes, the linear response and the weakly-nonlinear response begin to deviate from their linear and quadratic dependency on amplitude (\subfigrefb{fig:regimes}{C \& 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.
\subsection{Nonlinear signal transmission in low-CV P-units}
Weakly nonlinear responses are expected in cells with sufficiently low intrinsic noise levels, i.e. low baseline CVs \citep{Voronenko2017}. P-units fire action potentials probabilistically phase-locked to the self-generated EOD \citep{Bastian1981a}. Skipping of EOD cycles leads to the characteristic multimodal ISI distribution with maxima at integer multiples of the EOD period (\subfigrefb{fig:punit}{A}). In this example, the baseline ISI distribution has a CV$_{\text{base}}$ of 0.2, which is at the lower end of the P-unit population \citep{Hladnik2023}. Spectral analysis of the baseline activity shows two major peaks: the first is located at the baseline firing rate \fbase, the second is located at the discharge frequency \feod{} of the electric organ and is flanked by two smaller peaks at $\feod \pm \fbase{}$ (\subfigref{fig:punit}{B}).
\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$ and an intermediate CV$_{\text{base}}$ of its interspike intervals, as indicated. \figitem{B} Power spectral density of the cell's baseline response with peaks at the cell's baseline firing rate $r$ and the fish's EOD frequency $f_{\text{EOD}}$. \figitem{C} Random amplitude modulation 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 AM 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{eq:susceptibility}, 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 weaker. \figitem{G} Second-order susceptibilities projected onto the diagonal (means of 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$), quantifies the height of this peak relative to the values in the vicinity \notejb{See equation XXX}. \figitem{H} ISI distributions (top) and second-order susceptibilities (bottom) of four more example P-units (``2021-06-18-ae'', ``2012-03-30-ah'', ``2018-08-24-ak'', ``2018-08-14-ac'') covering the range of baseline firing rates and CV$_{\text{base}}$s as indicated. The first two cells show an anti-diagonal, but not the full expected triangular structure. The second-order susceptibilities of the other two cells are mostly flat and consequently the SI($r$) values are close to one.}
\end{figure*}
Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:punit}{C}, top trace, red line), are commonly used to characterize stimulus-driven responses of sensory neurons using transfer functions (first-order susceptibility), spike-triggered averages, or stimulus-response coherences. Here, we additionally estimate the second-order susceptibility to quantify nonlinear encoding. P-unit spikes align more or less clearly with fluctuations in the RAM stimulus. A higher stimulus intensity, here a higher contrast of the RAM relative to the EOD amplitude (see methods), entrains the P-unit response more clearly (light and dark purple for low and high contrast stimuli, \subfigrefb{fig:punit}{C}). Linear encoding, quantified by the transfer function, \eqnref{linearencoding_methods}, is similar for the two RAM contrasts in this low-CV P-unit (\subfigrefb{fig:punit}{D}), as expected for a linear system. The first-order susceptibility is low for low frequencies, peaks in the range below 100\,Hz and then falls off again \citep{Benda2005}.
The second-order susceptibility, \eqnref{eq:susceptibility}, quantifies for each combination of two stimulus frequencies \fone{} and \ftwo{} the amplitude and phase of the stimulus-evoked response at the sum \fsum{} (and also the difference, \subfigrefb{fig:model_full}{A}). Large values of the second-order susceptibility indicate stimulus-evoked peaks in the response spectrum at the summed frequency that cannot be explained by linear response theory. Similar to the first-order susceptibility, the second-order susceptibility can be estimated directly from the response evoked by a RAM stimulus that stimulates the neuron with a whole range of frequencies simultaneously (\subfigsref{fig:punit}{E, F}). For LIF and theta neuron models driven in the supra-threshold regime, theory predicts nonlinear interactions between the two stimulus frequencies, when the two frequencies \fone{} and \ftwo{} or their sum \fsum{} exactly match the neuron's baseline firing rate \fbase{} \citep{Voronenko2017,Franzen2023}. Only then, additional stimulus-evoked peaks appear in the spectrum of the spiking response that would show up in the second-order susceptibility as a horizontal, a vertical, and an anti-diagonal line (\subfigrefb{fig:lifresponse}{B}, pink triangle in \subfigsref{fig:punit}{E, F}).
For the low-CV P-unit, we observe a band of slightly elevated second-order susceptibility for the low RAM contrast at \fsumb{} (yellowish anti-diagonal between pink edges, \subfigrefb{fig:punit}{E}). This structure vanishes for the stronger stimulus (\subfigref{fig:punit}{F}). Further, the overall level of the second-order susceptibility is reduced with increasing stimulus strength. To quantify the structural changes in the susceptibility matrices we projected the susceptibility values onto the diagonal by averaging over the anti-diagonals (\subfigrefb{fig:punit}{G}). At low RAM contrast this projected second-order susceptibility indeed has a small peak at \fbase{} (\subfigrefb{fig:punit}{G}, dot on top line). For the higher RAM contrast, however, this peak vanishes and the overall level of the second-order susceptibility is reduced (\subfigrefb{fig:punit}{G}). The reason behind this reduction is that a RAM with a higher contrast is not only a stimulus with an increased amplitude, but also increases the total noise in the system. Increased noise is known to linearize signal transmission \citep{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017} and thus the second-order susceptibility is expected to decrease.
In contrast, a high-CV P-unit (CV$_{\text{base}}=0.4$) does not exhibit pronounced nonlinearities even at low stimulus contrasts (\figrefb{fig:punithighcv}). Irrespective of the CV, neither of the two example P-units shows the complete expected structure of nonlinear interactions.
\subsection{Ampullary afferents exhibit strong nonlinear interactions}
\begin{figure*}[p]
\includegraphics[width=\columnwidth]{ampullaryexamplecell}
\caption{\label{fig:ampullary} Linear and nonlinear stimulus encoding in example ampullary afferents. \figitem{A} Interspike interval (ISI) distribution of the cell's baseline activity (cell identifier ``2012-04-26-ae''). The very low CV of the ISIs indicates almost perfect periodic spiking. \figitem{B} Power spectral density of baseline activity with peaks at the cell's baseline firing rate and its harmonics. Ampullary afferents do not respond to the fish's EOD frequency, $f_{\text{EOD}}$. \figitem{C} Band-limited white noise stimulus (top, red, with a cutoff frequency of 150\,Hz) added to the fish's self-generated electric field (no amplitude modulation) and spike raster of the evoked responses (bottom) for two stimulus contrasts as indicated (right). \figitem{D} Gain of the transfer function, \eqnref{linearencoding_methods}, of the responses to stimulation with 5\,\% (light green) and 10\,\% contrast (dark green) of 10\,s duration. \figitem{E, F} Absolute value of the second-order susceptibility, \eqnref{eq:susceptibility}, for both stimulus contrasts as indicated. Both show a strong anti-diagonal where the two stimulus frequencies add up to the afferent's baseline firing rate. \figitem{G} Projections of the second-order susceptibilities in \panel{E, F} onto the diagonal. \figitem{H} ISI distributions (top) and second-order susceptibilites (bottom) of four more example afferents (``2010-11-26-an'', ``2011-10-25-ac'', ``2011-02-18-ab'', and ``2014-01-16-aj'').}
\end{figure*}
Electric fish possess an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogenous electric stimuli. The population of ampullary afferents is much less heterogeneous, and known for the much lower CVs of their baseline ISIs (CV$_{\text{base}}=0.06$ to $0.22$, \citealp{Grewe2017}). Ampullary cells do not phase-lock to the EOD and the ISIs are unimodally distributed (\subfigrefb{fig:ampullary}{A}). As a consequence of the high regularity of their baseline spiking activity, the corresponding power spectrum shows distinct peaks at the baseline firing rate \fbase{} and its harmonics. Since the cells do not respond to the self-generated EOD, there is no peak at \feod{} (\subfigrefb{fig:ampullary}{B}). When driven by a low-contrast noise stimulus (note: this is no longer an AM stimulus, \subfigref{fig:ampullary}{C}), ampullary cells exhibit very pronounced bands in the second-order susceptibility, where \fsum{} is equal to \fbase{} or its harmonic (yellow diagonals in \subfigrefb{fig:ampullary}{E}), implying strong nonlinear response components at these frequency combinations (\subfigrefb{fig:ampullary}{G}, top). With higher stimulus contrasts these bands disappear (\subfigrefb{fig:ampullary}{F}), the projection onto the diagonal loses its distinct peak at \fsum{} and its overall level is reduced (\subfigrefb{fig:ampullary}{G}, bottom).
\subsection{Model-based estimation of the second-order susceptibility}
In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagonal of the second-order susceptibility, where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, which is in line with theoretical expectations \citep{Voronenko2017,Franzen2023}. However, a pronounced nonlinear response at frequencies \foneb{} or \ftwob{}, although predicted by theory, cannot be observed. In the following, we investigate how these discrepancies can be understood.
\begin{figure*}[p]
\includegraphics[width=\columnwidth]{noisesplit}
\caption{\label{fig:noisesplit} Estimation of second-order susceptibilities. \figitem{A} \suscept{} (right) estimated from $N=198$ 256\,ms long FFT segments of an electrophysiological recording of another P-unit (cell ``2017-07-18-ai'', $r=78$\,Hz, CV$_{\text{base}}=0.22$) driven with a RAM stimulus with contrast 5\,\% (left). \figitem[i]{B} \textit{Standard condition} of model simulations with intrinsic noise (bottom) and a RAM stimulus (top). \figitem[ii]{B} \suscept{} estimated from simulations of the cell's LIF model counterpart (cell ``2017-07-18-ai'', table~\ref{modelparams}) based on a similar number of $N=100$ FFT segments. As in the electrophysiological recording only a weak anti-diagonal is visible. \figitem[iii]{B} Same as \panel[ii]{B} but using $10^6$ FFT segments. Now, the expected triangular structure is revealed. \figitem[iv]{B} Convergence of the \suscept{} estimate as a function of FFT segments. \figitem{C} At a lower stimulus contrast of 1\,\% the estimate did not converge yet even for $10^6$ FFT segments. \figitem[i]{D} Same as in \panel[i]{B} but in the \textit{noise split} condition: there is no external RAM signal (red) driving the model. Instead, a large part (90\,\%) of the total intrinsic noise is treated as a signal and is presented as an equivalent amplitude modulation ($s_{\xi}(t)$, orange), while the intrinsic noise is reduced to 10\,\% of its original strength (bottom, see methods for details). \figitem[i]{D} 100 FFT segments are still not sufficient for estimating \suscept{}. \figitem[iii]{D} Simulating one million segments reveals the full expected trangular structure of the second-order susceptibility. \figitem[iv]{D} In the noise-split condition, the \suscept{} estimate converges already at about $10^{4}$ FFT segments.}
\end{figure*}
%\notejb{Since the model overestimated the sensitivity of the real P-unit, we adjusted the RAM contrast to 0.9\,\%, such that the resulting spike trains had the same CV as the electrophysiological recorded P-unit during the 2.5\,\% contrast stimulation (see table~\ref{modelparams} for model parameters).} \notejb{chi2 scale is higher than in real cell}
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, i.e. repetitions of the same RAM stimulus, are available. As a consequence, the cross-spectra, \eqnref{eq:crosshigh}, are insufficiently averaged and the full structure of the second-order susceptibility might be hidden in finite-data noise. This experimental limitation can be overcome by using a computational model for the P-unit, a stochastic leaky integrate-and-fire model with adaptation current, dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of a low-CV cell estimated from the same low number of FFT (fast fourier transform) segments as in the experiment ($\n{}=11$, compare \panel{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 (\subfigrefb{fig:noisesplit}\,\panel[iii]{B}). In this example, the estimate of the second-order susceptibility indeed improves. It gets less noisy and the diagonal at \fsum{} is emphasized. However, the expected vertical and horizontal lines at \foneb{} and \ftwob{} are still mainly missing.
Using a broadband stimulus increases the effective input-noise level. This may linearize signal transmission and suppress potential nonlinear responses \citep{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017}. Assuming that the intrinsic noise level in this P-unit is small enough, the full expected structure of the second-order susceptibility should appear in the limit of weak AMs. Again, this cannot be done experimentally, because the problem of insufficient averaging becomes even more severe for weak AMs (low contrast). In the model, however, we know the time course of the intrinsic noise and can use this knowledge to determine the susceptibilities by input-output correlations via the Furutsu-Novikov theorem \citep{Furutsu1963, Novikov1965}. This theorem, in its simplest form, states that the cross-spectrum $S_{x\eta}(\omega)$ of a Gaussian noise $\eta(t)$ driving a nonlinear system and the system's output $x(t)$ is proportional to the linear susceptibility according to $S_{x\eta}(\omega)=\chi(\omega)S_{\eta\eta}(\omega)$. Here $\chi(\omega)$ characterizes the linear response to an infinitely weak signal $s(t)$ in the presence of the background noise $\eta(t)$. Likewise, the nonlinear susceptibility can be determined in an analogous fashion from higher-order input-output cross-spectra (see methods, equations \eqref{eq:crosshigh} and \eqref{eq:susceptibility}) \citep{Egerland2020}. In line with an alternative derivation of the Furutsu-Novikov theorem \citep{Lindner2022}, we can split the total noise and consider a fraction of it as a stimulus. This allows us to calculate the susceptibility from the cross-spectrum between the output and this stimulus fraction of the noise. Adapting this approach to our P-unit model (see methods), we replace the intrinsic noise by an approximately equivalent RAM stimulus $s_{\xi}(t)$ and a weak remaining intrinsic noise $\sqrt{2D \, c_{\rm{noise}}}\;\xi(t)$ with $c_\text{noise} = 0.1$ (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \subfigrefb{fig:noisesplit}\,\panel[i]{C}). We tune the amplitude of the RAM stimulus $s_{\xi}(t)$ such that the output firing rate and variability (CV) are the same as in the baseline activity (i.e. full intrinsic noise $\sqrt{2D}\;\xi(t)$ in the voltage equation but no RAM) and compute the cross-spectra between the RAM part of the noise $s_{\xi}(t)$ and the output spike train. This procedure has two consequences: (i) by means of the cross-spectrum between the output and \signalnoise, which is a large fraction of the noise, the signal-to-noise ratio of the measured susceptibilities is drastically improved; (ii) the total noise in the system has been reduced (by what was before the external RAM stimulus $s(t)$), which makes the system more nonlinear. For both reasons we now see the expected nonlinear features in the second-order susceptibility for a sufficient number of segments (\subfigrefb{fig:noisesplit}\,\panel[iii]{C}), but not for a number of segments comparable to the experiment (\subfigrefb{fig:noisesplit}\,\panel[ii]{C}). In addition to the strong response for \fsumb{}, we now also observe pronounced nonlinear responses at \foneb{} and \ftwob{} (vertical and horizontal lines, \subfigrefb{fig:noisesplit}\,\panel[iii]{C}).
Note, that the increased number of segments goes along with a substantial reduction of second-order susceptibility values (\subfigrefb{fig:noisesplit}\,\panel[iii]{C}). Only for more than about $10^5$ segments does the estimate of the second-order susceptibility converge for most of the model cells.
\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'' (table~\ref{modelparams}) 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{eq:nli_equation2}) 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 (table~\ref{modelparams}). Model cells have been visually categorized based on the 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.}
\end{figure}
\subsection{Weakly nonlinear interactions in many model cells}
In the previous section we have shown one example cell for which we find in the corresponding model strong ridges in the second-order susceptibility where one \notejg{same here...} or the sum of two stimulus frequencies match the neuron's baseline firing rate (\figrefb{fig:noisesplit}\,\panel[iii]{C}). Using our 39 P-unit models, we now can explore how many P-unit model neurons show such a triangular structure.
By just looking at the second-order susceptibilities estimated using the noise-split method (first column of \figrefb{fig:modelsusceptcontrasts}) we can readily identify strong triangular patterns in 11 of the 39 model cells (28\,\%, see \figrefb{fig:modelsusceptcontrasts}\,\panel[i]{A}\&\panel[i]{B} for two examples). In another 5 cells (13\,\%) the triangle is much weaker and sits on top of a smooth bump of elevated second-order susceptibility (\figrefb{fig:modelsusceptcontrasts}\,\panel[i]{C} shows an example). The remaining 23 model cells (59\,\%) show no triangle (see \figrefb{fig:modelsusceptcontrasts}\,\panel[i]{D} for an example). This categorization is supported by the peakedness of the nonlinearity, \nli{} \eqnref{eq:nli_equation2}, which quantifies the height of the ridge where the stimulus frequencies add up to the neuron's baseline firing rate. Values above one indicate an elevated ridge. The absence of such a ridge results in values close to one. Indeed, the cells showing only a weak triangle (orange) arise out of values around one and the cells showing strong triangles (red) have consistently \nli{} values exceeding 1.9 (\figrefb{fig:modelsusceptcontrasts}\,\panel[i]{E}).
The \nli{} correlates with the cell's CV of its baseline interspike intervals ($r=-0.59$, $p<0.001$). The lower the cell's CV, the higher the \nli{} value and thus the stronger the triangular structure of its second-order susceptibility. The model cells with the most distinct triangular pattern in their second-order susceptibility are the ones with the lowest CVs, hinting at low intrinsic noise levels.
\subsection{Weakly nonlinear interactions vanish for higher stimulus contrasts}
As pointed out in the context of stimulation with sine waves (\figrefb{fig:regimes}), the weakly nonlinear regime can only be observed for sufficiently weak stimulus amplitudes. In the model cells we estimated second-order susceptibilities for RAM stimuli with a contrast of 1, 3, and 10\,\%. For the 10\,\% contrast, $10^5$ FFT segments were sufficient for estimating the susceptibility (\figrefb{fig:trialnr}\,\panel[iv]{E}). The estimate converged similarly as in the noise-split configuration (\figrefb{fig:trialnr}\,\panel[i]{E}). For 3\,\% contrast, one million segments were just sufficient (\figrefb{fig:trialnr}\,\panel[iii]{E}) and for 1\,\% contrast even $10^7$ segments were sometimes not enough for a good estimate (\figrefb{fig:trialnr}\,\panel[ii]{E}). This again demonstrates the detrimental effect of smaller signal-to-noise ratios on the estimate of the second-order susceptibility and the advantage of the noise-split method.
The estimates for 1\,\% contrast (\figrefb{fig:modelsusceptcontrasts}\,\panel[ii]{E}), however, were quite similar to the estimates from the noise-split method, corresponding to a stimulus contrast of 0\,\% ($r=0.92$, $p\ll 0.001$). Thus, RAM stimuli with 1\,\% contrast are sufficiently small to not destroy weakly nonlinear interactions by their linearizing effect. 48\,\% of the model cells have a \nli{} value greater than 1.2.
At a RAM contrast of 3\,\% the \nli{} values become smaller (\figrefb{fig:modelsusceptcontrasts}\,\panel[iii]{E}). Only 9 cells (23\,\%) have \nli{} values exceeding 1.2. Finally, at 10\,\% the \nli{} values of all cells drop below 1.2, except for two cells (5\,\%, \figrefb{fig:modelsusceptcontrasts}\,\panel[iv]{E}). The cell shown in \subfigrefb{fig:modelsusceptcontrasts}{A} is one of them. Interestingly, the baseline CVs of these two cells (0.16 and 0.23) are not the lowest one (0.11) of the model cells, but their sensitivities are quite low (small response modulations and thus small output noise). At 10\,\% contrast the \nli{} values are no longer correlated with the ones in the noise-split configuration ($r=0.23$, $p=0.15$). 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:lifresponse}. \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}
Estimating second-order susceptibilities reliably requires large numbers (millions) of FFT segments (\figrefb{fig:trialnr}). Electrophysiological measurements, however, suffer from limited recording durations and hence limited numbers of available FFT segments and estimating weakly nonlinear interactions from just ten segments appears futile. The question arises, to what extend such limited-data estimates are still informative?
The second-order susceptibility matrices that are based on only 10 sgements look flat and noisy, lacking the triangular structure \subfigref{fig:modelsusceptcontrasts}{B}. The anti-diagonal ridge, however, when the sum of the stimulus frequencies matches the neuron's baseline firing rate seems to be present whenever the converged estimate shows a clear triangular structure (compare \subfigref{fig:modelsusceptcontrasts}{B} and \subfigref{fig:modelsusceptcontrasts}{a}).
%A comparison of well converged estimates based on millions of segments with estimates based on just ten segments as a worst-case scenario (\subfigrefb{fig:modelsusceptlown}{A\&B}) seems hopeless on a first glance. These estimates using just ten segments look flat and noisy, no triangular structure is visible. However, the anti-diagonal ridge where the stimulus frequencies add up to the neuron's baseline firing rate seems to be present when the converged estimate shows a clear triangular structure.
The \nli{} characterizes the ridgeness of the second-oder susceptibility plane. Comparing it when based on 10 versus one or ten million segments for all 39 model cells (\subfigrefb{fig:modelsusceptlown}{C}) supports this impression. As we have seen in the context of \subfigref{fig:modelsusceptcontrasts}{E} for converged estimates, values greater than one indicate triangular structures in the second-order susceptibility. The \nli{} values based on just ten segments correlate quite well with the ones from the converged estimates for contrasts 1\,\% and 3\,\% ($r=0.9$, $p<0.001$). At a contrast of 10\,\% this correlation is weaker ($r=0.53$, $p<0.001$), because there are only two cells left with \nli{} values greater than one. Despite the good correlations, care has to be taken to set a threshold on the \nli{} values for deciding whether a triangular structure would emerge for a much higher number of segments. Because the estimates are noisier, there could be false positives for a too low threshold. Setting the threshold to 1.6 avoids false positives for the price of a few false negatives.
Trying to predict whether there is a triangular structure in the noise-split configuration is more difficult (\subfigrefb{fig:modelsusceptlown}{D}). The correlations are weaker and at a stimulus contrast of 10\,\% the correlation is no longer significant. One false positive arises at a contrast of 1\,\%, and false positives are absent for higher contrasts. False negatives increase with increasing contrast and at a contrast of 10\,\% all \nli{} values based on 10 segments are so low that just one triangular pattern can be predicted. This makes sense given the absent correlation between \nli{} values estimated at 10\,\% stimulus contrast and the noise-split configuration described above.
Overall, observing \nli{} values greater than at least 1.6, even for a number of FFT segments as low as ten, seems to be a reliable indication for a triangular structure in the second-order susceptibility at the corresponding stimulus contrast. Small stimulus contrasts of 1\,\% are less informative, because of the bad signal-to-noise ratio. Intermediate stimulus contrasts around 3\,\% seem to be optimal, because there, most cells still have a triangular structure in their susceptibility and the signal-to-noise ratio is better. At RAM stimulus contrasts of 10\,\% or higher the signal-to-noise ratio is even better, but only few cells remain with weak triangularly shaped susceptibilities that might be missed as a false positives. Note that increasing the number of segments used for estimating the susceptibilities to 100 or 1000 improves the situation only marginally (not shown).
\begin{figure*}[tp]
\includegraphics[width=\columnwidth]{dataoverview}
\caption{\label{fig:dataoverview} Nonlinear responses in P-units and
ampullary afferents. The second-order susceptibility is condensed
into the susceptibility index, SI($r$) \eqnref{eq:nli_equation},
that quantifies the relative amplitude of the ridge where the two
stimulus frequencies add up to the cell's baseline firing rate $r$
(see \subfigrefb{fig:punit}{G}). In both the models and the
experimental data, the SI($r$) was estimated based on 100 FFT
segments \notejb{dt and nfft}. The SI($r$) is plotted against the
cells' CV of its baseline interspike intervals (left column), the
response modulation (the standard deviation of firing rate evoked
by the band-limited white-noise stimulus) --- a measure of
effective stimulus strength (center column), and the CV of the
interspike intervals during stimulation with the white-noise
stimulus (right column). Pearson's correlation coefficient $R$ and
the number of data points $n$ are indicated; all correlations are
significant at a level below $p=0.0006$. Kernel-density estimates
of the distributions of the displayed quantities are plotted on
top and right. Data points are color coded by CV$_{\text{base}}$
or response modulation as indicated by the color bars. The
horizontal dashed line marks a threshold for SI($r$) values at 1.8
and the percentages to the right denote the fractions above and
below this threshold. \figitem{A} The SI($r$) of all 39 model
P-units (table~\ref{modelparams}) measured at contrasts of 1, 3,
and 10\,\% of RAM stimuli with a cutoff frequency of 300\,Hz. The
black square marks the cell from \subfigrefb{fig:noisesplit}{C},
the circles the four cells shown in
\subfigref{fig:modelsusceptcontrasts}{A--D}, and the triangle the
cell from \subfigref{fig:modelsusceptlown}{A--B}.
\figitem{B}
Electrophysiological data from 159 P-units. Each cell contributes
on average with 2 (min. 1, max. 6) RAM stimulus presentations to
the $n=329$ data points. The RAMs with cutoff frequency of 300\,Hz
were presented at contrasts ranging from 2.5\,\% to 20\,\%
(average 8\,\%). The two black triangles mark the responses of the
example P-unit from \subfigrefb{fig:punit}{E,F}, the circles the
other four examples from \subfigrefb{fig:punit}{H}, and the
triangle the unit from \subfigrefb{fig:noisesplit}{A}.
\figitem{C} Recordings from 30 ampullary afferents, each
contributing on average 3 (min. 1, max. 7) RAM stimulus
presentations to $n=89$ data points. Stimuli had a cutoff
frequency of 150\,Hz and their contrasts ranged from 2.5\,\% to
20\,\% (average 7\,\%). The two black triangles mark the
responses of the example ampullary afferent from
\subfigrefb{fig:ampullary}{E,F}, and the circles the other four
examples from \subfigrefb{fig:ampullary}{H}.}
\end{figure*}
\subsection{Low CVs and weak stimuli are associated with distinct nonlinearity in recorded electroreceptive neurons}
Now we are prepared to evaluate our pool of 221 P-units and 47 ampullary afferents recorded in 71 specimen.
\notejb{We need to state the number of trials}
\notejb{We need to know which contrasts}
% All the statements about nonlinear encoding in p-type and ampullary electroreceptor afferents based on single-cell examples shown above are supported by the analysis of our pool of 221 P-units and 47 ampullary afferents recorded in 71 specimen.
\notejb{\nli{} explanation is needed earlier}
For comparison across cells we summarize the structure of the second-order susceptibilities in a single number, the peakedness of the projected nonlinearity at \fbase{} (\nli{}) \eqnref{eq:nli_equation}, that characterizes the size of the expected peak of the projections of a \suscept{} matrix onto its diagonal at the baseline firing rate (e.g. \subfigref{fig:punit}{G}). \nli{} assumes high values when the peak at \fbase{} is pronounced relative to the median of projections onto the diagonal and is small when there is no distinct peak.
Three P-units stand out with \nli{} values exceeding two, but additional \notejb{XXX} cells have \nli{} values greater than 1.8. Based on our insights from the P-unit models these would have the expected triangular structure in their susceptibilities when estimated with a sufficiently high number of segments. However, we can only speculate how many of the cells with lower \nli{} values are false negatives. \notejb{because also many have been measured at too strong contrasts}.
The \nli{} values of the P-unit population correlate weakly with the CV of the baseline ISIs. Cells with lower baseline CVs tend to have more pronounced peaks in their projections than those that have high baseline CVs (\subfigrefb{fig:dataoverview}{A}). This negative correlation is more pronounced against the CV measured during stimulation (\subfigrefb{fig:dataoverview}{C}).
The effective stimulus strength also plays an important role. 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. P-units are heterogeneous in their sensitivity, their response modulations to the same stimulus contrast vary a lot \citep{Grewe2017}. Cells with weak responses to a stimulus, be it an insensitive cell or a weak stimulus, have higher \nli{} values and thus a more pronounced ridge in the second-order susceptibility at \fsumb{} in comparison to strongly responding cells that basically have flat second-order susceptibilities (\subfigrefb{fig:dataoverview}{E}). How pronounced nonlinear response components are in P-units thus depends on the baseline CV (a proxy for the internal noise level), and both the CV and response strength during stimulation (effective output noise).
%(Pearson's $r=-0.35$, $p<0.001$)221 P-units and 47 (Pearson's $r=-0.16$, $p<0.01$)
%In a P-unit population where each cell is represented not by several contrasts but by the lowest recorded contrast, \nli{} significantly correlates with the CV during baseline ($r=-0.17$, $p=0.01$), the response modulation ($r=-0.35$, $p<0.001$) and \fbase{} ($r=-0.32$, $p<0.001$).%, $\n{}=221$*, $\n{}=221$******, $\n{}=221$
The population of ampullary cells is generally more homogeneous, with lower baseline CVs than P-units. Accordingly, \nli{} values of ampullary cells are indeed much higher than in P-units by about a factor of ten. \notejb{XXX (percent)} ampullary cells with values greater than 1.8 would have a triangular structure in their second-order susceptibilities. Ampullary cells also show a negative correlation with baseline CV. Again, sensitive cells with strong response modulations are at the bottom of the distribution and have \nli{} values close to one (\subfigrefb{fig:dataoverview}{B, D}). The weaker the response modulation, because of less sensitive cells or weaker stimulus amplitudes, the stronger the nonlinear component of a cell's response (\subfigrefb{fig:dataoverview}{F}).
%(Pearson's $r=-0.35$, $p < 0.01$) (Pearson's $r=-0.59$, $p < 0.0001$)
\begin{figure*}[t]
\includegraphics[width=\columnwidth]{model_full.pdf}
\caption{\label{fig:model_full} Second-order susceptibility qualitatively predicts responses to sine-wave stimuli. \figitem[]{A} Absolute value of the second-order susceptibility, \eqnref{eq:susceptibility}, for both positive and negative frequencies. \susceptf{} was estimated from $N=10^6$ FFT segments of model simulations in the noise-split condition (cell ``2013-01-08-aa'', table~\ref{modelparams}). Dashed white lines mark zero frequencies. Nonlinear responses at \fsum{} are quantified in the upper right and lower left quadrants. Nonlinear responses at \fdiff{} are quantified in the upper left and lower right quadrants. The baseline firing rate of this cell was at $\fbase=120$\,Hz. The position of the orange/red letters corresponds to the beat frequencies used for the stimulation with pure sine waves in the subsequent panels and indicates the sum/difference of those beat frequencies. \figitem{B--E} Black line -- power spectral density of model simulations in response to stimulation with two pure sine waves, \fone{} and \ftwo, in addition to the receiving fish's own EOD (three-fish scenario). The contrast of beat beats is 2\,\%. Colored circles highlight the height of selected peaks in the power spectrum. Grey line -- power spectral density of model in the baseline condition. \figitem{B} The sum of the two beat frequencies match \fbase{}. \figitem{C} The difference of \fone{} and \ftwo{} match \fbase{}. \figitem{D} Only the first beat frequency matches \fbase{}. \figitem{E} None of the two beat frequencies matches \fbase{}.}
\end{figure*}
\subsection{Second-order susceptibility can explain nonlinear peaks in pure sinewave stimulation}
Using the RAM stimulation we found pronounced nonlinear responses in the limit to weak stimulus amplitudes. How do these findings relate to the situation of pure sinewave stimuli with finite amplitudes that approximate the interference of EODs of real animals? For the P-units, the relevant signals are the beat frequencies \fone{} and \ftwo{} that arise from the interference of either of the two foreign EODs with the receiving fish's own EOD (\figref{fig:motivation}). In the introductory example, the response power spectrum showed peaks from nonlinear interactions at the sum of the two beat frequencies (orange marker, \subfigrefb{fig:motivation}{D}) and at the difference between the two beat frequencies (red marker, \subfigrefb{fig:motivation}{D}). In this example, \ftwo{} was similar to \fbase{}, corresponding to the horizontal line of the second-order susceptibility estimated for a vanishing external RAM stimulus (\subfigrefb{fig:noisesplit}\,\panel[iii]{C}). In the three-fish example, there was a second nonlinearity at the difference between the two beat frequencies (red dot, \subfigrefb{fig:motivation}{D}), that is not covered by the so-far shown part of the second-order susceptibility (\subfigrefb{fig:noisesplit}\,\panel[iii]{C}), in which only the response at the sum of the two stimulus frequencies is addressed. % less prominent,
However, the second-order susceptibility \eqnref{eq:susceptibility} is a spectral measure that is based on Fourier transforms and thus is also defined for negative stimulus frequencies. The full \susceptf{} matrix is symmetric with respect to the origin. In the upper-right and lower-left quadrants of \susceptf{}, stimulus-evoked responses at \fsum{} are shown, whereas in the lower-right and upper-left quadrants nonlinear responses at the difference \fdiff{} are shown (\figref{fig:model_full}). The vertical and horizontal lines at \foneb{} and \ftwob{} are very pronounced in the upper-right quadrant of \subfigrefb{fig:model_full}{A} for nonlinear responses at \fsum{} and extend into the upper-left quadrant (representing \fdiff), where they fade out towards more negative $f_1$ frequencies. The peak in the response power-spectrum at \fdiff{} evoked by pure sine-wave stimulation (\subfigrefb{fig:motivation}{D}) is predicted by the horizontal line in the upper-left quadrant (\subfigrefb{fig:model_full}{A}, \citealp{Schlungbaum2023}).
Is it possible to predict nonlinear responses in a three-fish setting based on second-order susceptibility matrices estimated from RAM stimulation (\subfigrefb{fig:model_full}{A})? We test this by stimulating the same model with two weak beats (\subfigrefb{fig:model_full}{B--E}). If we choose a frequency combination in which the sum of the two beat frequencies is equal to the model's baseline firing rate \fbase{}, a peak at the sum of the two beat frequencies appears in the power spectrum of the response (\subfigrefb{fig:model_full}{B}), as expected from \suscept. If instead we choose two beat frequencies that differ from \fbase{}, a peak is present at the difference frequency (\subfigrefb{fig:model_full}{C}). If only one of the beat frequency matches \fbase{}, both, a peak at the sum and at the difference frequency are present in the P-unit response (\subfigrefb{fig:model_full}{D}). And if none of these conditions are met, neither a peak at the sum nor at the difference of the two beat frequencies appears (\subfigrefb{fig:model_full}{E}).
\notejb{This works only qualitatively. In fact, the beat matrix looks quite a bit different and also heavily depends on contrast. This will be addressed in a future paper.}
\section{Discussion}
Theoretical work \citep{Voronenko2017,Franzen2023} derived analytical expressions for weakly-nonlinear responses in LIF and theta model neurons driven by two sine waves with distinct frequencies. We here investigated such nonlinear responses in two types of electroreceptor afferents that differ in their intrinsic noise levels \citep{Grewe2017} using white-noise stimuli to estimate second-order susceptibilities. Following \citet{Voronenko2017} we expected to observe increased levels of second-order susceptibility where either of the stimulus frequencies alone or the sum of the stimulus frequencies matches the baseline firing rate ($f_1=\fbase{}$, $f_2=\fbase{}$ or \fsumb{}). We find traces of these nonlinear responses in most of the low-noise ampullary afferents and only those P-units with very low intrinsic noise levels. Complementary model simulations demonstrate, in the limit to vanishing stimulus amplitudes and extremely high number of repetitions, that the second order susceptibilities estimated from the electrophysiological data are indeed indicative of the theoretically expected weakly nonlinear responses. With this, we provide experimental evidence for nonlinear responses of a spike generator at low stimulus amplitudes.
% EOD locking:
% Phase-locking to the own field also leads to a representation of \feod{} in the P-unit firing rate (see \figref{fig:punit}\panel{B}) \citep{Sinz2020}.
% Weakly nonlinear responses versus saturation regime
\subsection{Intrinsic noise limits nonlinear responses}
The weakly nonlinear regime with its triangular pattern of elevated second-order susceptibility resides between the linear and a stochastic mode-locking regime. Too strong intrinsic noise linearizes the system and wipes out the triangular structure of the second-order susceptibility \citep{Voronenko2017}. Our electrophysiological recordings match this theoretical expectation. Only P-units with low coefficients of variation (CV $<$ 0.25) of the interspike-interval distribution in their baseline response show the expected nonlinearities (\figref{fig:punit}, \figref{fig:model_full}, \subfigref{fig:dataoverview}{A}). Such low-CV cells are rare among the 221 P-units used in this study. On the other hand, the majority of the ampullary cells have generally lower CVs (median of 0.12) and have an approximately ten-fold higher level of second-order susceptibilities where \fsumb{} (\figref{fig:ampullary}, \subfigrefb{fig:dataoverview}{B}).
The CV is a proxy for the intrinsic noise in the cells (\citealp{Vilela2009}, note however the effect of coherence resonance for excitable systems close to a bifurcation, \citealp{Pikovsky1997, Lindner2004}). In both cell types, we observe a negative correlation between the second-order susceptibility at \fsumb{} and the CV (\figrefb{fig:dataoverview}), indicating that it is the level of intrinsic noise that shapes nonlinear responses. These findings are in line with previous theoretical and experimental studies showing the linearizing effects of noise in nonlinear systems \citep{Roddey2000, Chialvo1997, Voronenko2017}. Increased intrinsic noise has been demonstrated to increase the CV and to reduce nonlinear phase-locking in vestibular afferents \citep{Schneider2011}. Reduced noise, on the other hand, has been associated with stronger nonlinearity in pyramidal cells of the ELL \citep{Chacron2006}. Only in cells with sufficiently low levels of intrinsic noise, weakly nonlinear responses can be observed.
\subsection{Linearization by white-noise stimulation}
Not only the intrinsic noise but also the stimulation with external white-noise linearizes the cells. This applies to both, the stimulation with AMs in P-units (\subfigrefb{fig:dataoverview}{E}) and direct stimulation in ampullary cells (\subfigrefb{fig:dataoverview}{F}). The stronger the effective stimulus, the less pronounced are the peaks in second-order susceptibility (see \subfigref{fig:punit}{E\&F} for a P-unit example and \subfigref{fig:ampullary}{E\&F} for an ampullary cell). This linearizing effect of noise stimuli limits the weakly nonlinear regime to small stimulus amplitudes. At higher stimulus amplitudes, however, other nonlinearities of the system eventually show up in the second-order susceptibility.
In order to characterize weakly nonlinear responses of the cells in the limit to vanishing stimulus amplitudes we utilized the Furutsu-Novikov theorem \citep{Novikov1965, Furutsu1963}. Following \citet{Lindner2022}, a substantial part of the intrinsic noise of a P-unit model \citep{Barayeu2023} is treated as signal. Performing this noise-split trick we can estimate the weakly nonlinear response without the linearizing effect of an additional external white noise stimulus. The model of a low-CV P-unit then shows the full nonlinear structure (\figref{fig:noisesplit}) known from analytical derivations and simulations of basic LIF and theta models driven with pairs of sine-wave stimuli \citep{Voronenko2017,Franzen2023}. Previous studies on second-order nonlinearities have not observed the weakly nonlinear regime, probably because of the linearizing effects of strong noise stimuli \citep{Victor1977, Schanze1997, Neiman2011}.
\subsection{Characterizing nonlinear coding from limited experimental data}
Estimating the Volterra series from limited experimental data is usually restricted to the first two or three orders, which might not be sufficient for a proper prediction of the neuronal response \citep{French2001}. A proper estimation of just the second-order susceptibility in the weakly nonlinear regime is challenging in electrophysiological experiments. Making assumptions about the nonlinearities in a system reduces the amount of data needed for parameter estimation. In particular, models combining linear filtering with static nonlinearities \citep{Chichilnisky2001} have been successful in capturing functionally relevant neuronal computations in visual \citep{Gollisch2009} as well as auditory systems \citep{Clemens2013}. Linear methods based on backward models for estimating the stimulus from neuronal responses, have been extensively used to quantify information transmission in neural systems \citep{Theunissen1996, Borst1999, Wessel1996, Machens2001}, because backward models do not need to generate action potentials that involve strong nonlinearities \citep{Rieke1999}.
\subsection{Nonlinear encoding in ampullary cells}
The afferents of the passive electrosensory system, the ampullary cells, exhibit strong second-order susceptibilities (\figref{fig:dataoverview}). Ampullary cells more or less directly translate external low-frequency electric fields into afferent spikes, much like in the standard LIF and theta models used by Lindner and colleagues \citep{Voronenko2017,Franzen2023}. Indeed, we observe in the ampullary cells similar second-order nonlinearities as in LIF models.
Ampullary stimuli originate from the muscle potentials induced by prey movement \citep{Kalmijn1974, Engelmann2010, Neiman2011fish}. For a single prey items such as \textit{Daphnia}, these potentials are often periodic but the simultaneous activity of a swarm of prey resembles Gaussian white noise \citep{Neiman2011fish}. Linear and nonlinear encoding in ampullary cells has been studied in great detail in the paddlefish \citep{Neiman2011fish}. The power spectrum of the baseline response shows two main peaks: One peak at the baseline firing frequency, a second one at the oscillation frequency of primary receptor cells in the epithelium, and interactions of both. Linear encoding in the paddlefish shows a gap at the epithelial oscillation frequency, instead, nonlinear responses are very pronounced there.
Ampullary 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 ampullaries increases with stimulus intensity while it disappears in our case (\subfigrefb{fig:dataoverview}{F}) indicating that paddlefish data have been recorded above the weakly-nonlinear regime.
The population of ampullary cells is very homogeneous with respect to the baseline rate (131$\pm$29\,Hz) and stimulus encoding properties \citep{Grewe2017}. This implies that, if the stimulus contains the appropriate frequency components that sum up to the baseline rate, there should be a nonlinear response at a frequency that is similar in the full population of ampullary cells (at the baseline firing frequency) that is outside the linear coding range. Postsynaptic cells integrating ampullary input might be able to extract such nonlinear responses. How such nonlinear effects might influence prey detection should be addressed in future studies.
\subsection{Nonlinear encoding in P-units}
Noise stimuli have the advantage that a range of frequencies can be measured with a single stimulus presentation and they have been successfully applied to characterize sensory coding in many systems \citep{French1973, Marmarelis1999, Borst1999, Chacron2005, Grewe2017}. In the real-world situation of wave-type electric fish, however, the natural stimuli encoded by P-units are periodic amplitude modulations of the self-generated electric field which arise from the superposition of the own and foreign EODs. Natural interactions usually occur between low numbers of close-by fish and thus the AMs are a mixture of a few distinct frequencies with substantial amplitudes \citep{Stamper2010,Fotowat2013, Henninger2020}. How informative are the second-order susceptibilities observed under noise stimulation for the encoding of distinct frequencies?
We here applied broad-band noise amplitude-modulation stimuli to estimate the second-order susceptibility.
%The total signal power in noise stimuli is uniformly distributed over a wide frequency band while the power is spectrally focused in pure sinewave stimuli.
Broadband noise stimuli introduce additional noise that linearizes the dynamics of the system. In contrast, pure sinewave stimulation is spectrally focussed and drives the system on the background of the intrinsic noise. This explains why we can observe nonlinear interactions between sinewave stimuli with distinct frequencies (\figrefb{fig:motivation}) although they are barely visible when stimulating with strong broad-band noise (\figrefb{fig:noisesplit}). As discussed above in the context of broad-band noise stimulation, we expect interactions with a cell's baseline frequency only for the very few P-units with low enough CVs of the baseline interspike intervals.
%We thus conclude that the presence of the anti-diagonal pattern in the \suscept{} matrix is sufficient to infer that the same nonlinear interactions happen here, in accordance with the the second-order susceptibilities found in LIF models without a carrier \citep{Voronenko2017, Schlungbaum2023}. This also validates the application of the white noise approach to characterize the full \suscept{} matrix instead of using all combinations of frequencies.
Extracting the AM from the stimulus already 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 in the limit of vanishing amplitudes and for stimulation with two distinct frequencies. Both settings have not been studied yet.
The encoding of secondary or social envelopes that arise from relative movement or the interaction of more than two animals \citep{Stamper2012Envelope} requires an additional nonlinearity in the system that was initially 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 electrosensory cocktail party. 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, \figref{fig:dataoverview}\panel{A}, \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 low-CV P-units 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.
\subsection{Conclusions}
We have demonstrated pronounced nonlinear responses in primary electrosensory afferents at weak stimulus amplitudes and sufficiently low intrinsic noise levels. The observed nonlinearities match the expectations from previous theoretical studies \citep{Voronenko2017,Franzen2023}. The resulting nonlinear components introduce spectral components not present in the original stimulus, but may provide an edge in the context of signal detection problems at stimulus amplitudes close to threshold \citep{Schlungbaum2023}. Electrosenory 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.
\section{Acknowledgements}
\notejb{Supported by DFG SPP XXX}
Tim Hladnik, Henriette Walz, Franziska Kuempfbeck, Fabian Sinz, Laura Seidler, Eva Vennemann, and Ibrahim Tunc recorded data.
\section{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 221 P-units and 47 ampullary electroreceptor afferents recorded in 71 weakly electric fish of the species \lepto{}. Electrophysiological recordings were performed on male and female weakly electric fish 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.
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{Experimental setup}
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{www.relacs.net}) running on a Linux computer was used for online spike and EOD detection, stimulus generation, and calibration. Recorded data was then stored on the hard drive for offline analysis.
\subsection{Identification of P-units and ampullary cells}
The neurons were classified into cell types during the recording by the experimenter. P-units were classified based on baseline firing rates of 50--450\,Hz, a clear phase-locking to the EOD, and by their responses to amplitude modulations of their own EOD \citep{Grewe2017, Hladnik2023}. Ampullary cells were classified based on firing rates of 80--200\,Hz, absent phase-locking to the EOD, and responses to low-frequency sinusoidal stimuli \citep{Grewe2017}. We here selected only those cells of which the neuron's baseline activity as well as the responses to frozen noise stimuli were recorded.
\subsection{Electric field recordings}
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.
\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, 300 or 400\,Hz. The stimulus intensity is given as the contrast, i.e. the standard deviation of the white noise stimulus in relation to the fish's EOD amplitude. The contrast varied between 1 and 20\,\%. Only cells with at least 10\,s of white noise stimulation were included for 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 EOD of the fish was multiplied with the desired random amplitude modulation profile (MXS-01M; npi electronics).
% and between 2.5 and 40\,\% for \eigen
\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}).
%sklearn \citep{scikitlearn2011},
\paragraph{Baseline analysis}\label{baselinemethods}
The baseline firing rate \fbase{} was calculated as the number of spikes divided by the duration of the baseline recording (on average 18\,s). The coefficient of variation (CV) of the interspike intervals (ISI) is their standard deviation relative to the average ISI: $\rm{CV} = \sqrt{\langle (ISI- \langle ISI \rangle) ^2 \rangle} / \langle ISI \rangle$. If the baseline was recorded several times in a recording, the average \fbase{} and CV were calculated.
\paragraph{White noise analysis} \label{response_modulation}
In the stimulus-driven case, the neuronal activity of the recorded cell is modulated around the average firing rate that is similar to \fbase{} and in that way encodes the time-course of the stimulus. Spiking activity
\begin{equation}
\label{eq:spikes}
x_k(t) = \sum_i\delta(t-t_{k,i})
\end{equation}
is recorded for each stimulus presentation $k$, as a train of times $t_{k,i}$ where action potentials occured. If only a single trial was recorded or is used for the analysis, we drop the trial index $k$.
The single-trial firing rate
\begin{equation}
r_k(t) = x_k(t) * K(t) = \int_{-\infty}^{+\infty} x_k(t') K(t-t') \, \mathrm{d}t'
\end{equation}
was estimated by convolving the spike train with a kernel. We used a Gaussian kernel
\begin{equation}
K(t) = {\scriptstyle \frac{1}{\sigma\sqrt{2\pi}}} e^{-\frac{t^2}{2\sigma^2}}
\end{equation}
with standard deviation $\sigma$ set to 2.5\,ms if not stated otherwise. Averaging over $n$ repeated stimulus presentations results in the trial averaged firing rate
\begin{equation}
\label{eq:rate}
r(t) = \left\langle r_k(t) \right\rangle _k = \frac{1}{n} \sum_{k=1}^n r_k(t)
\end{equation}
The average firing rate during stimulation, $r_s = \langle r(t) \rangle_t$, is given by the temporal average $\langle \cdot \rangle_t$ over the duration of the stimulus of the trial-averaged firing rate. To quantify how strongly a neuron was driven by the stimulus, we computed the response modulation as the standard deviation $\sigma_{s} = \sqrt{\langle (r(t)-r_s)^2\rangle_t}$ of the trial-averaged firing rate.
\paragraph{Spectral analysis}\label{susceptibility_methods}
The neuron is driven by the stimulus and thus its spiking response depends on the time course of the stimulus. To characterize the relation between stimulus $s(t)$ and response $x(t)$, we calculated the first- and second-order susceptibilities in the frequency domain.
Fast fourier transforms (FFT) $\tilde s_T(\omega)$ and $\tilde x_T(\omega)$ of $s(t)$ and $x(t)$, respectively, were computed according to $\tilde x_T(\omega) = \int_{0}^{T} \, x(t) e^{- i \omega t}\,dt$ for $T=0.5$\,s long segments with no overlap, resulting in a spectral resolution of 2\,Hz for the experimental data. For simplicity we use angular frequencies $\omega=2\pi f$ instead of frequencies $f$. In the experimental data, most stimuli had a duration of 10\,s and were chopped into 20 segments. Spectral measures were computed for single trials of neural responses, series of spike times, \eqnref{eq:spikes}. For spectral parameters used in simulations, see below.
\notejb{Update FFT parameter. Describe binary spiketrain normalized by dt and plain FFT. Power spectra scaled by dt/nfft, bispectra by dt**2/nfft.}
The power spectrum of the stimulus $s(t)$ was estimated as
\begin{equation}
\label{powereq}
S_{ss}(\omega) = \frac{\langle \tilde s_T(\omega) \tilde s_T^*(\omega)\rangle}{T}
\end{equation}
with $\tilde s^*$ being the complex conjugate of $\tilde s$ and $\langle \cdot \rangle$ denoting the average over the segments. The power spectrum of the spike trains $S_{xx}(\omega)$ was estimated accordingly. The cross-spectrum $S_{xs}(\omega)$ between stimulus and evoked spike trains was estimated according to
\begin{equation}
\label{cross}
S_{xs}(\omega) = \frac{\langle \tilde x_T(\omega) \tilde s_T^*(\omega)\rangle}{T}
\end{equation}
The first-order susceptibility (transfer function)
\begin{equation}
\label{linearencoding_methods}
\chi_{1}(\omega) = \frac{S_{xs}(\omega) }{S_{ss}(\omega) }
\end{equation}
was then computed from $S_{xs}(\omega)$ and $ S_{ss}(\omega)$.
The second-order cross-spectrum
\begin{equation}
\label{eq:crosshigh}
S_{xss} (\omega_{1},\omega_{2}) = \frac{\langle \tilde x_T(\omega_{1} + \omega_{2}) \tilde s_T^*(\omega_{1}) \tilde s_T^*(\omega_{2}) \rangle}{T}
\end{equation}
describes nonlinear interactions that generate responses at the sum and difference (at negative $\omega_2$) of two stimulus frequencies $\omega_1$ and $\omega_2$.
The second-order susceptibility
\begin{equation}
\label{eq:susceptibility}
\chi_{2}(\omega_{1}, \omega_{2}) = \frac{S_{xss} (\omega_{1},\omega_{2})}{2 S_{ss} (\omega_{1}) S_{ss} (\omega_{2})}
\end{equation}
normalizes the second-order cross-spectrum by the spectral power at the two stimulus frequencies.
% Applying the Fourier transform this can be rewritten resulting in:
% \begin{equation}
% \label{susceptibility}
% \chi_{2}(\omega_{1}, \omega_{2}) = \frac{TN \sum_{n=1}^N \int_{0}^{T} dt\,r_{n}(t) e^{-i(\omega_{1}+\omega_{2})t} \int_{0}^{T}dt'\,s_{n}(t')e^{i \omega_{1}t'} \int_{0}^{T} dt''\,s_{n}(t'')e^{i \omega_{2}t''}}{2 \sum_{n=1}^N \int_{0}^{T} dt\, s_{n}(t)e^{-i \omega_{1}t} \int_{0}^{T} dt'\,s_{n}(t')e^{i \omega_{1}t'} \sum_{n=1}^N \int_{0}^{T} dt\,s_{n}(t)e^{-i \omega_{2}t} \int_{0}^{T} dt'\,s_{n}(t')e^{i \omega_{2}t'}}
% \end{equation}
Throughout the manuscript we only show the absolute values of the complex-valued second-order susceptibility matrix and ignore the corresponding phases.
\paragraph{Nonlinearity index}\label{projected_method}
We expected to see elevated values in the second-order susceptibility at $\omega_1 + \omega_2 = \fbase$ \citep{Voronenko2017,Franzen2023}. To characterize this in a single number we computed the peakedness of the nonlinearity via
\begin{equation}
\label{eq:nli_equation}
\nli{} = \frac{\max D(\fbase{}-5\,\rm{Hz} \leq f \leq \fbase{}+5\,\rm{Hz})}{\mathrm{median}(D(f))}
\end{equation}
from the second-order susceptibilities estimated from experimental data. $D(f)$ is the projection of the absolute values of the second-order susceptibility matrix onto the diagonal, calculated by taking the mean of the anti-diagonal elements. The peak of $D(f)$ at the neuron's baseline firing rate $\fbase$ was found by finding the maximum of $D(f)$ in the range $\fbase \pm 5$\,Hz. This peak was then normalized by the median of $D(f)$ (\subfigrefb{fig:punit}{G}). Since in most experimentally measured cells the second-order susceptibilities was more or less flat, normalizing by the median is working well. If the same RAM was recorded several times in a cell, each trial resulted in a separate second-order susceptibility matrix. For the population statistics in \figref{fig:dataoverview} the mean of the resulting \nli{} values is used.
For the model simulations, the second-order susceptibilities where not flat but often showed a broader bump on which the peak at the neuron's baseline firing rate occured (see, for example, \figrefb{fig:modelsusceptcontrasts}, panels \panel[i]{B}, \panel[iii]{C}, and \panel[iv]{D}). Instead of normalizing by the median, we normalized the model simulations by the mean values of the projection $D(f)$ in small ranges close to the left and right of the neuron's baseline firing rate:
\begin{equation}
\label{eq:nli_equation2}
\nli{} = \frac{\max D(\fbase{}-5\,\rm{Hz} \leq f \leq \fbase{}+5\,\rm{Hz})}{\frac{1}{2}\left[\mathrm{mean}(D(\fbase{}-20\,\rm{Hz} \leq f \leq \fbase{}-10\,\rm{Hz})) + \mathrm{mean}(D(\fbase{}+10\,\rm{Hz} \leq f \leq \fbase{}+20\,\rm{Hz}))\right]}
\end{equation}
Using \eqnref{eq:nli_equation} instead of \eqnref{eq:nli_equation2} for the estimates of the second-order susceptibility based on $N=10$ segments does not change the results in \subfigref{fig:modelsusceptlown}{C \& D}. That is, experimental data based on low numbers of segments can be safely analyzed using the more robust \eqnref{eq:nli_equation}.
\begin{figure*}[t]
\includegraphics[width=\columnwidth]{flowchart.pdf}
\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 (table~\ref{modelparams}). \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 $\fbase$, $f_{EOD} - \fbase$, $f_{EOD}$, and $f_{EOD} + \fbase$. With RAM stimulation (\panel[ii]{E}), the peak at the baseline firing rate, $\fbase$, is washed out.}
\end{figure*}
\subsection{Leaky integrate-and-fire models for P-units}\label{lifmethods}
Modified leaky integrate-and-fire (LIF) models were constructed to reproduce the specific firing properties of P-units \citep{Chacron2001, Sinz2020}. The sole driving input into the P-unit model during baseline, i.e. when no external stimulus was given, is the fish's own EOD, modeled as a cosine wave
\begin{equation}
\label{eq:eod}
\carrierinput = y_{EOD}(t) = \cos(2\pi f_{EOD} t)
\end{equation}
with EOD frequency $f_{EOD}$ and an amplitude normalized to one.
To mimic the interaction with other fish, the EODs of a second or third fish with EOD frequencies $f_1$ and $f_2$, respectively, were added to the normalized EOD, \eqnref{eq:eod}, of the receiving fish according to their contrasts, $c_1$ and $c_2$ at the position of the receiving fish:
\begin{equation}
\label{eq:modelbeats}
\carrierinput = \cos(2\pi f_{EOD} t) + c_1 \cos(2\pi f_1 t) + c_2\cos(2\pi f_2 t)
\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
\begin{equation}
\label{eq:ram_equation}
y(t) = (1+ s(t)) \cos(2\pi f_{EOD} t)
\end{equation}
The contrast $c$ of the RAM is the standard deviation of the RAM relative to the amplitude of the receiving fish.
First, the input \carrierinput{} is thresholded by setting negative values to zero:
\begin{equation}
\label{eq:threshold2}
\lfloor \carrierinput \rfloor_0 = \left\{ \begin{array}{rcl} \carrierinput & ; & \carrierinput \ge 0 \\ 0 & ; & \carrierinput < 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
\begin{equation}
\label{eq:dendrite}
\tau_{d} \frac{d V_{d}}{d t} = -V_{d}+ \lfloor \carrierinput \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 dendritic voltage $V_d(t)$ is then fed into a stochastic leaky integrate-and-fire (LIF) model with adaptation,
\begin{equation}
\label{eq:LIF}
\tau_{m} \frac{d V_{m}}{d t} = - V_{m} + \mu + \alpha V_{d} - A + \sqrt{2D}\;\xi(t)
\end{equation}
where $\tau_{m}$ is the membrane time-constant, $\mu$ is a fixed bias current, $\alpha$ is a scaling factor for $V_{d}$, $A$ is an inhibiting adaptation current, and $\sqrt{2D}\;\xi(t)$ is a white noise with strength $D$. Note, that all state variables, membrane voltages $V_d$ and $V_m$ as well as the adaptation current, are dimensionless.
The adaptation current $A$ follows
\begin{equation}
\label{eq:adaptation}
\tau_{A} \frac{d A}{d t} = - A
\end{equation}
with adaptation time constant $\tau_A$ \citep{Benda2003,Benda2005,Benda2010}.
Whenever the membrane voltage $V_m(t)$ crosses the spiking threshold $\theta=1$, a spike was generated, $V_{m}(t)$ was reset to $0$, the adaptation current was incremented by $\Delta A$, and integration of $V_m(t)$ was paused for the duration of a refractory period $t_{ref}$ (\subfigrefb{flowchart}{D}):
\begin{equation}
\label{spikethresh}
V_m(t) \ge \theta \; : \left\{ \begin{array}{rcl} V_m & \mapsto & 0 \\ A & \mapsto & A + \Delta A/\tau_A \end{array} \right.
\end{equation}
The P-unit models were integrated by the Euler forward method with a time-step of $\Delta t = 0.05$\,ms. The intrinsic noise $\xi(t)$ (\eqnref{eq:LIF}, \subfigrefb{flowchart}{C}) was added by drawing a random number from a normal distribution $\mathcal{N}(0,\,1)$ with zero mean and standard deviation of one in each time step $i$. This number was multiplied with $\sqrt{2D}$ and divided by $\sqrt{\Delta t}$. For each trial of a simulation, the variables $A$, $V_{d}$ and $V_{m}$ were drawn from a distribution of initial values, estimated from a 100\,s long simulation of baseline activity after an initial 100\,s long integration that was discarded as a transient.
%\begin{equation}
% \label{eq:LIFintegration}
% V_{m_{i+1}} = V_{m_i} + \left(-V_{m_i} + \mu + \alpha V_{d_i} - A_i + \sqrt{\frac{2D}{\Delta t}}\mathcal{N}(0,\,1)_i\right) \frac{\Delta t}{\tau_m}
%\end{equation}
%\paragraph{Fitting the model to recorded P-units}
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 (table~\ref{modelparams}).
\notejb{add table with all 39 cells}
\begin{table*}[tp]
\caption{\label{modelparams} Model parameters of LIF models, fitted to 3 electrophysiologically recorded P-units \citep{Ott2020}.}
\begin{tabular}{lrrrrrrrr}
\hline
\bfseries cell & \bfseries $\mu$ & \bfseries $\beta$ & \bfseries $\tau_{m}$/ms & \bfseries $D$/$\mathbf{ms}$ & \bfseries $\tau_{A}$/ms & \bfseries $\Delta_A$ & \bfseries $\tau_{d}$/ms & \bfseries $t_{ref}$/ms \\\hline
2012-07-03-ak& $-1.32$& $10.6$& $1.38$& $0.001$& $96.05$& $0.01$& $1.18$& $0.12$ \\
2013-01-08-aa& $0.59$& $4.5$& $1.20$& $0.001$& $37.52$& $0.01$& $1.18$& $0.38$ \\
2018-05-08-ae& $-21.09$& $139.6$& $1.49$& $0.214$& $123.69$& $0.16$& $3.93$& $1.31$ \\
\hline
\end{tabular}
\end{table*}
In \figrefb{fig:noisesplit} and \figrefb{fig:model_full}, for each of the stimulus and response segments needed for the spectral analysis, \eqnsref{powereq}--\eqref{eq:susceptibility}, a simulation was run. The first second was discarded and the analysis was based on the last second of the data. The resulting spectra thus have a spectral resolution of 1\,Hz. To speed up model simulations for \figrefb{fig:modelsusceptcontrasts} and \figrefb{fig:trialnr} we simulated up to $10^6$ trials of 3\,s duration. In each trial the first 0.5\,s were skipped, yielding 10 FFT segments of $T=0.25$\,s (4\,Hz resolution).
\subsection{Noise split}
\label{intrinsicsplit_methods}
Based on the Furutsu-Novikov theorem \citep{Furutsu1963,Novikov1965,Lindner2022,Egerland2020}, we split the total noise, $\sqrt{2D}\;\xi(t)$, of a LIF model, \eqnref{eq:LIF}, into two parts. The first part is the intrinsic noise term, $\sqrt{2D \, c_{\rm{noise}}}\;\xi(t)$, whose strength is reduced by a factor $c_{\rm{noise}}=0.1$ (\subfigrefb{flowchart}\,\panel[iii]{C}). The second part replaces the now missing intrinsic noise by a driving input signal $s_{\xi}(t)$, a RAM stimulus with frequencies up to 300\,Hz (\subfigrefb{flowchart}\,\panel[iii]{A}). The LIF model with splitted noise then reads
\begin{eqnarray}
\label{eq:ram_split}
y(t) & = & (1+ s_\xi(t)) \cos(2\pi f_{EOD} t) \\
\label{eq:Noise_split_intrinsic_dendrite}
\tau_{d} \frac{d V_{d}}{d t} & = & -V_{d}+ \lfloor y(t) \rfloor_{0} \\
\label{eq:Noise_split_intrinsic}
\tau_{m} \frac{d V_{m}}{d t} & = & - V_{m} + \mu + \alpha V_{d} - A + \sqrt{2D \, c_{\rm{noise}}}\;\xi(t)
\end{eqnarray}
Both, the reduced intrinsic noise and the RAM stimulus, need to replace the original intrinsic noise. Because the RAM stimulus is band-limited and undergoes some transformations before it is added to the reduced intrinsic noise, it is not \textit{a priori} clear, what the amplitude of the RAM should be. We bisected the amplitude of $s_\xi(t)$, until the CV of the resulting interspike intervals matched the one of the original model's baseline activity. The second-order cross-spectra, \eqnref{eq:crosshigh}, were computed between the RAM stimulus $s_{\xi}(t)$ and the spike train $x(t)$ it evoked. In this way, the effective signal-to-noise ratio can be increased while maintaining the total noise in the system.
%\bibliographystyle{apalike}%alpha}%}%alpha}%apalike}
\bibliography{journalsabbrv,references}
% \bibliographystyle{apalike} %or any other style you like
%\bibliography{references}
%\bibliography{journalsabbrv,references}
\end{document}