diff --git a/README.md b/README.md index 6170a24..41f4ff0 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,12 @@ -Tabele with all cell information considering RAM analysis is saved in data/ as data_overview_mod.csv +Tabele with all cell information considering RAM overview analysis is saved in data/ as data_overview_mod.csv ``` .py frame = pd.read_csv('data_overview_mod.csv') ``` see data_overview_mod.py for loading and filtering. + +The susceptibility model is calculated in threefish\threefish\RAM\calc_model with calc_RAM_model() +The data single cell suceptiblity matrices are calculated in threefish\threefish\RAM\calc_model with calc_RAM_model() and the overview for data_overview_mod are calculated with calc_RAM_overview() + +You cannot use the develop folder (e.g. to change the cells chosen) from the susceptiblity folder. But all data and the folder structure is on anaconda in /home/rudnaya/Desktop/Masterarbeit/work/code +There you can rerun all scripts (the code is in there) and the data is saved automatically, and can be then retrieved in the susceptibility or twobeat folder by runing the redo paramter if defining yourself as the develop folder +See README in /home/rudnaya/Desktop/Masterarbeit/work/code \ No newline at end of file diff --git a/Settup.pdf b/Settup.pdf deleted file mode 100644 index 34ae5b1..0000000 Binary files a/Settup.pdf and /dev/null differ diff --git a/__pycache__/model_and_data.cpython-39.pyc b/__pycache__/model_and_data.cpython-39.pyc deleted file mode 100644 index 64bdeeb..0000000 Binary files a/__pycache__/model_and_data.cpython-39.pyc and /dev/null differ diff --git a/__pycache__/model_full.cpython-39.pyc b/__pycache__/model_full.cpython-39.pyc deleted file mode 100644 index 7127c7f..0000000 Binary files a/__pycache__/model_full.cpython-39.pyc and /dev/null differ diff --git a/__pycache__/plotstyle.cpython-39.pyc b/__pycache__/plotstyle.cpython-39.pyc deleted file mode 100644 index 4cbb046..0000000 Binary files a/__pycache__/plotstyle.cpython-39.pyc and /dev/null differ diff --git a/ampullary.pdf b/ampullary.pdf index d1b15dd..f85b166 100644 Binary files a/ampullary.pdf and b/ampullary.pdf differ diff --git a/cells_suscept.pdf b/cells_suscept.pdf index ef30935..af7d897 100644 Binary files a/cells_suscept.pdf and b/cells_suscept.pdf differ diff --git a/cells_suscept_high_CV.pdf b/cells_suscept_high_CV.pdf index 175af63..1989a6c 100644 Binary files a/cells_suscept_high_CV.pdf and b/cells_suscept_high_CV.pdf differ diff --git a/data_overview_mod.pdf b/data_overview_mod.pdf index 86d58d4..9d27092 100644 Binary files a/data_overview_mod.pdf and b/data_overview_mod.pdf differ diff --git a/flowchart.pdf b/flowchart.pdf index 4ea60f0..2715452 100644 Binary files a/flowchart.pdf and b/flowchart.pdf differ diff --git a/model_and_data.pdf b/model_and_data.pdf index 9bf98da..db6c70c 100644 Binary files a/model_and_data.pdf and b/model_and_data.pdf differ diff --git a/model_full.pdf b/model_full.pdf index 2551094..16ca0ae 100644 Binary files a/model_full.pdf and b/model_full.pdf differ diff --git a/motivation.pdf b/motivation.pdf index 0b87762..ff2a1ee 100644 Binary files a/motivation.pdf and b/motivation.pdf differ diff --git a/nonlin_regime.pdf b/nonlin_regime.pdf index d65406c..4f0e45d 100644 Binary files a/nonlin_regime.pdf and b/nonlin_regime.pdf differ diff --git a/plot_chi2.pdf b/plot_chi2.pdf index 6f85cbf..0abfcdb 100644 Binary files a/plot_chi2.pdf and b/plot_chi2.pdf differ diff --git a/references.bib b/references.bib index 76cb8fc..428fcf1 100644 --- a/references.bib +++ b/references.bib @@ -657,6 +657,18 @@ pages = {811--824} } +@ARTICLE{Benda2003, + AUTHOR = {Jan Benda and Andreas V. M. Herz}, + TITLE = {A universal model for spike-frequency adaptation.}, + YEAR = {2003}, + JOURNAL = NeuralComput, + VOLUME = {15}, + NUMBER = {11}, + PAGES = {2523--2564}, + PDF = {Benda2003a.pdf}, + URL = {http://neco.mitpress.org/cgi/content/abstract/15/11/2523} +} + @ARTICLE{Benda2005, AUTHOR = {Jan Benda and Andr\'e Longtin and Leonard Maler}, TITLE = {Spike-frequency adaptation separates transient communication signals from background oscillations.}, @@ -7001,3 +7013,20 @@ microelectrode recordings from visual cortex and functional implications.}, doi={10.1109/10.16447}} +@Article{Billah1990, + title={Numerical method for colored-noise generation and its application to a bistable system.}, + author={Billah, K. and Shinozuka, M.}, + year={1990}, + journal={Phys. Rev. A}, + volume={42}, + pages={7492–-7495}, +} + +@Article{Skorjanc2023, + title={Stimulator compensation and generation of Gaussian noise stimuli with defined amplitude spectra for studying input–output relations of sensory systems.}, + author={Ale\v{s} \v{S}korjanc and Marko Kreft and Jan Benda}, + year={2023}, + journal={J Comp Physiol A}, + volume={209}, + pages={361–-372}, +} diff --git a/susceptibility1.dvi b/susceptibility1.dvi deleted file mode 100644 index 5d3ef6d..0000000 Binary files a/susceptibility1.dvi and /dev/null differ diff --git a/susceptibility1.tex b/susceptibility1.tex index 3311140..cb31ce1 100644 --- a/susceptibility1.tex +++ b/susceptibility1.tex @@ -5,8 +5,8 @@ \author{Alexandra Barayeu\textsuperscript{1}, Maria Schlungbaum\textsuperscript{2,3}, Benjamin Lindner\textsuperscript{2,3}, - Jan Benda\textsuperscript{1, 4} - Jan Grewe\textsuperscript{1, $\dagger$}} + Jan Grewe\textsuperscript{1}, + Jan Benda\textsuperscript{1, 4, $\dagger$}} \date{\normalsize \textsuperscript{1} Institute for Neurobiology, Eberhard Karls Universit\"at T\"ubingen, Germany\\ @@ -208,13 +208,8 @@ %%%%% equation references %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \renewcommand{\eqref}[1]{(\ref{#1})} -\newcommand{\eref}[1]{\ref{#1}} -\newcommand{\eqn}{Eq} -\newcommand{\Eqn}{Eq} -\newcommand{\eqns}{Eqs} -\newcommand{\Eqns}{Eqs} -\newcommand{\Eqnref}[1]{\Eqn~\eqref{#1}} -\newcommand{\Eqnsref}[1]{\Eqns~\eqref{#1}} +\newcommand{\eqnref}[1]{Eq.~\eqref{#1}} +\newcommand{\eqnsref}[1]{Eqs.~\eqref{#1}} %%%%% species names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newcommand{\Lepto}{\emph{Apteronotus leptorhynchus}} @@ -395,7 +390,7 @@ %\paragraph{Conflict of interest:}The authors declare no conflict of interest. %\paragraph{Author Contributions:} All authors designed the study and discussed the results. AB performed the data analyses and modeling, AB and JG drafted the paper, all authors discussed and revised the manuscript. -\paragraph{Keywords:} population coding $|$ conduction delay $|$ heterogeneity $|$ electric fish $|$ mutual information +\paragraph{Keywords:} $|$ $|$ heterogeneity $|$ electric fish $|$ mutual information % Please keep the abstract below 300 words @@ -413,58 +408,54 @@ Neuronal processing is inherently nonlinear --- spiking thresholds or rectificat \begin{figure*}[t] \includegraphics[width=\columnwidth]{plot_chi2.pdf} + \notejg{units on the figure?} \caption{\label{fig:lifresponse} First- (linear) and second-order response functions of the leaky integrate-and-fire model. \figitem{A} Magnitude of the first-order response function $|\chi_1(f_1)|$, also known as ``gain'' quantifies the response amplitude relative to the stimulus amplitude, both measured at the stimulus frequency. \figitem{B} Magnitude of the second-order 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 summed frequencies. The plots show the analytical solutions from \citep{Lindner2001} and \citep{Voronenko2017} with $\mu = 1.1$ and $D = 0.001$.} \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 nervous 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. -Nonlinear processes are key to neuronal information processing. Decision-making is a fundamentally 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 used to describe 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. -At the heart of nonlinear system identification is the Volterra series \citep{Rieke1999}. Second-order kernels have been used to predict firing rate responses of catfish retinal ganglion cells \citep{Marmarelis1972}. -In the frequency domain, second-order kernels are known as second-order response functions or susceptibilities. They quantify the amplitude of the response at the sum and difference of two stimulus frequencies. Adding also third-order kernels, spike trains of spider mechanoreceptors have been predicted from sensory stimuli \citep{French2001}. The nonlinear nature of Y cells in contrast to the more linear responses of X cells in cat retinal ganglion cells has been demonstrated using second-order kernels \citep{Victor1977}. Interactions between different frequencies in the response of neurons in visual cortices of cats and monkeys have been studied using bispectra, the crucial constituent of the second-order susceptibility \citep{Schanze1997}. Locking of chinchilla auditory nerve fibers to pure tone stimuli is captured by second-order kernels \citep{Temchin2005}. In paddlefish ampullary afferents, bursting in response to strong, natural sensory stimuli boosts nonlinear responses in the bicoherence, the bispectrum normalized by stimulus and response spectra \citep{Neiman2011}. - -Noise linearizes nonlinear systems \citep{Yu1989, Chialvo1997} and therefore noisy neural systems can be well described by linear response theory in the limit of small stimulus amplitudes \citep{Roddey2000, Doiron2004, Rocha2007, Sharafi2013}. When increasing stimulus amplitude, at first the contribution of the second-order kernel of the Volterra series becomes relevant in addition to the linear one. For these weakly nonlinear responses of leaky-integrate-and-fire (LIF) neurons, an analytical expression for the second-order susceptibility has been derived \citep{Voronenko2017} in addition to its linear response function \citep{Lindner2001}. In the superthreshold regime, where the LIF generates a baseline firing rate in the absence of an external stimulus, the linear response function has a peak at the baseline firing rate and its harmonics (\subfigrefb{fig:lifresponse}{A}) and the second-order susceptibility shows very distinct ridges of elevated nonlinear responses, exactly where two stimulus frequencies equal or add up to the neuron's baseline firing rate (\subfigrefb{fig:lifresponse}{B}). In experimental data, such structures in the second-order susceptibility have not been reported yet. - -Here we study weakly nonlinear responses in two electrosensory systems in the wave-type electric fish \textit{Apteronotus leptorhynchus}. These fish generate a quasi-sinusoidal dipolar electric field (electric organ discharge, EOD). In communication contexts \citep{Walz2014, Henninger2018} the EODs of close-by fish superimpose and lead to amplitude modulations (AMs), called beats (two-fish interaction), and modulations of beats, called envelopes (multiple-fish interaction) \citep{Yu2005, Fotowat2013}. Therefore, stimuli with multiple distinct frequencies are part of the everyday life of wave-type electric fish \citep{Benda2020} and interactions of these frequencies in the electrosensory periphery are to be expected. P-type electroreceptor afferents of the tuberous electrosensory system, the P-units, use nonlinearities to extract and encode these AMs in their time-dependent firing rates \citep{Bastian1981a, Walz2014, Middleton2006, Stamper2012Envelope, Savard2011, Barayeu2023}. P-units are heterogeneous in their baseline firing rates as well as in their intrinsic noise levels, as quantified by the coefficient of variation (CV) of the interspike intervals (ISI) \citep{Grewe2017, Hladnik2023}. Low-CV P-units have a less noisy firing pattern that is closer to pacemaker firing, whereas high-CV P-units show a more irregular firing pattern that is closer to a Poisson process. On the other hand, ampullary cells of the passive electrosensory system are homogeneous in their response properties and have very low CVs \citep{Grewe2017}. - -Field observations have shown that courting males were able to react to distant intruder males despite the strong EOD of the nearby female \citep{Henninger2018}. Weakly nonlinear responses are of direct relevance in this setting since they could boost responses to the faint signal of the distant intruder and thus improve detection \citep{Schlungbaum2023}. +Noise in nonlinear systems, however, linearizes the system's response properties \citep{Yu1989, Chialvo1997}. Alternatively, 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 of leaky-integrate-and-fire (LIF) model neurons, an analytical expression for the second-order susceptibility has been derived \citep{Lindner2001, Voronenko2017}. In the suprathreshold regime, where the LIF generates a baseline firing rate in the absence of an external stimulus, the linear response function has a peak at the baseline firing rate and its harmonics (\subfigrefb{fig:lifresponse}{A}) and the second-order susceptibility shows very distinct ridges of elevated nonlinear responses, exactly where two stimulus frequencies equal or add up to the neuron's baseline firing rate (\subfigrefb{fig:lifresponse}{B}). In experimental data, such structures in the second-order susceptibility have not been reported yet. +Here we search for such weakly nonlinear responses in electroreceptors of the two electrosensory systems of the wave-type electric fish \textit{Apteronotus leptorhynchus}, i.e. the tuberous (active) and the ampullary (passive) electrosensory system. The p-type electroreceptors of the active system (P-units) are driven by the fish's high-frequency, quasi-sinusoidal electric organ discharges (EOD) and encode disturbances of it \citep{Bastian1981a}. The electroreceptors of the passive system are tuned to lower-frequency exogeneous electric fields such as caused by muscle activity of prey \citep{Kalmijn1974}. As different animals have different EOD-frequencies, being exposed to stimuli of multiple distinct frequencies is part of the animal's everyday life \citep{Benda2020,Henninger2020} and weakly non-linear interactions may occur in the electrosensory periphery. In communication contexts \citep{Walz2014, Henninger2018} the EODs of interacting fish superimpose and lead to periodic amplitude modulations (AMs or beats) of the receiver's EOD. Non-linear mechanisms in P-units, enable encoding of AMs in their time-dependent firing rates \citep{Bastian1981a, Walz2014, Middleton2006, Barayeu2023}. When multiple animals interact, the EOD interferences induce second-order amplitude modulations referred to as envelopes \citep{Yu2005, Fotowat2013, Stamper2012Envelope} and saturation non-linearities allow also for the encoding of these in the electrosensory periphery \citep{Savard2011}. Field observations have shown that courting males were able to react to the extremely weak signals of distant intruding males despite the strong foreground EOD of the nearby female \citep{Henninger2018}. Weakly nonlinear interactions at particular combinations of signals can be of immediate relevance in such settings as they could boost detectability of the faint signals \citep{Schlungbaum2023}. +%The population of P-units is heterogeneous with respect to their baseline firing rates as well as in their intrinsic noise levels, as quantified by the coefficient of variation (CV) of the interspike intervals (ISI) \citep{Grewe2017, Hladnik2023}. Low-CV P-units have a less noisy firing pattern that is closer to pacemaker firing, whereas high-CV P-units show a more irregular firing pattern that is closer to a Poisson process. The population of ampullary cells of the passive electrosensory system, on the other hand, is homogeneous in their response properties and CVs are low \citep{Grewe2017}. %\notejb{in Voronenko they talk about second-order response functions and not of susceptibilities}\notesr{Ja das ist wahr, aber so wie ich das verstehe sind das synonyme, oder etwa nicht?} - \section{Results} \begin{figure*}[t] \includegraphics[width=\columnwidth]{motivation.pdf} - \caption{\label{fig:motivation} Nonlinearity in an electrophysiologically recorded P-unit of \lepto{} in a three-fish setting (cell identifier ``2021-08-03-ac"). Receiver with EOD frequency $\feod{} =664$\,Hz encounters fish with EOD frequencies $f_{1}=631$\,Hz and $f_{2}=797$\,Hz. Both encountered fish lead to a beat contrast of 10\,\%. Top: Scheme of a nonlinear system. Second row: Interference of the receiver EOD with the EODs of other fish. Third row: Spike trains of the P-unit. Fourth row: Firing rate, retrieved as the convolution of the spike trains with a Gaussian kernel ($\sigma = 1$\,ms). Bottom row: Power spectrum of the firing rate. \figitem{A} Baseline condition: Only the receiver is present. The baseline firing rate \fbase{} dominates the power spectrum of the firing rate. \figitem{B} The receiver and the fish with EOD frequency $f_{1}=631$\,Hz are present. \figitem{C} The receiver and the fish with EOD frequency $f_{2}=797$\,Hz are present. \figitem{D} All three fish with the EOD frequencies \feod{}, $f_{1}$ and $f_{2}$ are present. Nonlinear peaks occur at the sum and difference of the two beat frequencies in the power spectrum of the firing rate. + \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 super-threshold regime with low stimulus amplitudes (\figrefb{fig:lifresponse}, \citealp{Voronenko2017}). Here, we re-analyze a large set of recordings of P-units and ampullary cells of the active and passive electrosensory systems of the brown ghost knifefish \textit{Apteronotus leptorhynchus} together with simulations of LIF-based models of P-unit spiking to search for such weakly nonlinear responses in real neurons. We start with a few example P-units and models to demonstrate the basic concepts. - +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. Our work is supported by simulations of LIF-based models of P-unit spiking. We start with demonstrating the basic concepts using example P-units and models. -\subsection{Nonlinear responses in P-units stimulated with two beat frequencies} -Without any external stimulation, a P-unit fires action potentials at a spontaneous baseline rate \fbase{} to the fish's own EOD of frequency \feod{}. Accordingly, a peak at \fbase{} is present in the power spectrum of this baseline activity (\subfigrefb{fig:motivation}{A}). Superposition of the receiver's EOD with an EOD of another fish with frequency $f_1$ results in a beat, a periodic amplitude modulation of the receiver's EOD. The frequency of the beat is given by the difference frequency $\Delta f_1 = f_1 - \feod$ between the two fish. P-units encode this beat in their firing rate \citep{Bastian1981a, Barayeu2023} and consequently, a peak at this beat frequency appears the the power spectrum of the response (\subfigrefb{fig:motivation}{B}). A second peak at the first harmonic indicates a nonlinear response that here is easily identified by the clipping of the P-unit's firing rate at zero. Pairing the fish with another fish with a higher beat frequency $\Delta f_2 = f_2 - \feod$ results in a weaker response with a single peak in the response power spectrum (\subfigrefb{fig:motivation}{C}). Note that $\Delta f_2$ has been chosen to match the P-unit's baseline firing rate. +\subsection{Nonlinear responses in P-units stimulated with two beat frequencies}\notejg{stimulated with two beats?} +Without external stimulation, a P-unit is driven by the fish's own EOD alone (with 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 \notejg{(only valid for $f_1 < f_{EOD}/2)$}. 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 with 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 that $\Delta f_2$ has been deliberately chosen to match the recorded P-unit's baseline firing rate. -When stimulating the fish with both frequencies, additional peaks appear in the response power spectrum at the sum \bsum{} and the difference frequency \bdiff{} (\subfigrefb{fig:motivation}{D}). Thus, the response is not equal to the sum of the responses to the two beats presented separately. These peaks at the sum and the difference of the two stimulus frequencies are a hallmark of nonlinear interactions that by definition are absent in linear systems. +When stimulating with both foreign signals simultaneously, additional peaks appear in the response power spectrum at the sum \bsum{} and the difference frequency \bdiff{} (\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. \subsection{Linear and weakly nonlinear regimes} \begin{figure*}[tp] \includegraphics[width=\columnwidth]{nonlin_regime.pdf} - \caption{\label{fig:nonlin_regime} Linear and nonlinear responses of a model P-unit in a three-fish setting in dependence on increasing stimulus amplitudes. The model P-unit was stimulated with two sinewaves of equal amplitude (contrast) at difference frequencies $\bone=30$\,Hz and $\btwo=130$\,Hz relative the receiver's EOD frequency. \btwo{} was set to match the baseline firing rate \fbase{} of the P-unit. The model used has the cell identifier 2013-01-08-aa (table~\ref{modelparams}). \figitem{A--D} Top: the stimulus, an amplitude modulation of the receiver's EOD resulting from the stimulation with the two sine waves. The contrasts of both beats increase from \panel{A} to \panel{D} as indicated. Middle: Spike raster of the model P-unit response. Bottom: power spectrum of the firing rate estimated from the spike raster with a Gaussian kernel ($\sigma=1$\,ms). \figitem{A} At very low stimulus contrasts the response is linear. The only peaks in the response spectrum are at the two stimulating beat frequencies (green and purple marker). \figitem{B} At higher stimulus contrasts, nonlinear responses appear at the sum and the difference of the stimulus frequencies (orange and red marker). \figitem{C} At even stronger stimulation additional peaks appear in the power spectrum. \figitem{D} At a contrast of 10\,\% the response at the sum o fthe stimulus frequencies almost disappears. \figitem{E} Amplitude of the linear (\bone{}, \btwo{}) and nonlinear (\bdiff{}, \bsum{}) responses of the model P-units as a function of beat contrast (contrasts increase equally for both beats).} + \caption{\label{fig:nonlin_regime} Linear and nonlinear responses of a model P-unit in a three-fish setting in dependence on increasing stimulus amplitudes. The model P-unit was stimulated with two sinewaves of equal amplitude (contrast) at difference frequencies $\bone=30$\,Hz and $\btwo=130$\,Hz relative the receiver's EOD frequency. \btwo{} was set to match the baseline firing rate \fbase{} of the P-unit. The model used has the cell identifier 2013-01-08-aa (table~\ref{modelparams}). \figitem{A--D} Top: the stimulus, an amplitude modulation of the receiver's EOD resulting from the stimulation with the two sine waves. The contrasts of both beats increase from \panel{A} to \panel{D} as indicated. Middle: Spike raster of the model P-unit response. Bottom: power spectrum of the firing rate estimated from the spike raster with a Gaussian kernel ($\sigma=1$\,ms). \figitem{A} At very low stimulus contrasts the response is linear. The only peaks in the response spectrum are at the two stimulating beat frequencies (green and purple marker). \figitem{B} At higher stimulus contrasts, nonlinear responses appear at the sum and the difference of the stimulus frequencies (orange and red marker). \figitem{C} At even stronger stimulation additional peaks appear in the power spectrum. \figitem{D} At a contrast of 10\,\% the response at the sum of the stimulus frequencies almost disappears. \figitem{E} Amplitude of the linear (\bone{}, \btwo{}) and nonlinear (\bdiff{}, \bsum{}) responses of the model P-units as a function of beat contrast (contrasts increase equally for both beats).} \end{figure*} -The beat stimuli in the example shown in \figref{fig:motivation} had just one amplitude that was not small. Whether this falls into the weakly nonlinear regime is not clear. In order to illustrate how the responses to two beat frequencies develop over a range of amplitudes we use an stochastic leaky-integrate-and-fire (LIF) model \citep{Barayeu2023}, that faithfully reproduces baseline firing properties and responses to step stimuli of a specific electrophysiologically measured P-unit. +The stimuli used in \figref{fig:motivation} had the same not-small amplitude. Whether this stimulus conditions falls into the weakly nonlinear regime as in \citet{Voronenko2017} is not clear. In order to illustrate how the responses to two beat frequencies develop over a range of amplitudes we use a stochastic leaky-integrate-and-fire (LIF) based P-unit model fitted to a specific electrophysiologically measured cell \citep{Barayeu2023}. + +At very low stimulus contrasts (less than approximately 0.5\,\% relative to the receiver's EOD amplitude) the spectrum has small peaks only at the beat frequencies (note, $|\Delta f_2 = f_{base}|$, \subfigref{fig:nonlin_regime}{A}, green and purple). The amplitudes of these peaks initially increase linearly with stimulus amplitude (\subfigref{fig:nonlin_regime}{E}), an indication of the linear response at lowest stimulus amplitudes. -At very low stimulus contrasts smaller than about 0.5\,\% only peaks at the stimulating beat frequencies appear in the response spectrum (\subfigref{fig:nonlin_regime}{A}). The amplitude of these peaks initially increase linearly with stimulus amplitude (\subfigref{fig:nonlin_regime}{E}). This is the linear response at low stimulus amplitudes. +The linear regime is followed by the weakly non-linear regime. In addition to peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:nonlin_regime}{B}, orange and red). The amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:nonlin_regime}{E}). -The linear regime is followed by the weakly non-linear regime. In addition to peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:nonlin_regime}{B}). THe amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:nonlin_regime}{E}). +At higher stimulus amplitudes (\subfigref{fig:nonlin_regime}{C \& D}) additional peaks appear in the response spectrum. The linear response and the weakly-nonlinear response start to deviate from their linear and quadratic dependence on amplitude (\subfigref{fig:nonlin_regime}{E}). The responses may even decrease for intermediate stimulus contrasts (\subfigref{fig:nonlin_regime}{D}). At high stimulus contrasts, additional nonlinearities, in particular clipping of the firing rate, in the system shape the responses. -At higher stimulus amplitudes (\subfigref{fig:nonlin_regime}{C \& D}) additional peaks appear in the response. The linear response and the weakly-nonlinear response start to deviate from their linear and waudratic dependence on amplitude (\subfigref{fig:nonlin_regime}{E}). The responses may even decrease for intermediate stimulus contrasts (\subfigref{fig:nonlin_regime}{D}). -For this example, we have chosen two specific stimulus (beat) frequencies. In particular one matching the P-unit's baseline firing rate. In the following, however, we are interested on how the non-linear responses depend on different combinations of stimulus frequencies in the weakly-nonlinear regime. +For this example, we have chosen two specific stimulus (beat) frequencies. One of these matching the P-unit's baseline firing rate. In the following, however, we are interested in how the non-linear responses depend on different combinations of stimulus frequencies in the weakly nonlinear regime. \subsection{Nonlinear signal transmission in low-CV P-units} @@ -473,43 +464,43 @@ Weakly nonlinear responses are expected in cells with sufficiently low intrinsic \begin{figure*}[t] \includegraphics[width=\columnwidth]{cells_suscept.pdf} -\caption{\label{fig:cells_suscept} Linear and nonlinear stimulus encoding in a low-CV P-unit (cell identifier ``2010-06-21-ai"). \figitem{A} Interspike interval (ISI) distribution of the cell's baseline activity, i.e. the cell is driven only by the unperturbed own electric field. The low CV of the ISIs indicates quite regular firing. \figitem{B} Power spectral density of the baseline response with peaks at the cell's baseline firing rate \fbase{} and the fish's EOD frequency \feod{}. \figitem{C} Random amplitude modulation stimulus (top, with cutoff frequency of 300\,Hz) and evoked responses (spike raster, bottom) of the same P-unit. The stimulus contrast (right) measures the strength of the AM. \figitem{D} Gain of the transfer function (first-order susceptibility), \Eqnref{linearencoding_methods}, computed from the responses to 10\,\% (light purple) and 20\,\% contrast (dark purple) RAM stimulation of 10\,s duration. \figitem{E, F} Absolute value of the second-order susceptibility, \Eqnref{eq:susceptibility}, for both the low and high stimulus contrast. Pink triangles mark vertical, horizontal, and diagonal lines where \fone, \ftwo{} or \fsum{} are equal to \fbase{}. \figitem{G} Second-order susceptibilities projected onto the diagonal (means of all anti-diagonals of the matrices shown in \panel{E, F}). Dots mark \fbase{}, horizontal dashed lines mark medians of the projected susceptibilities.} +\caption{\label{fig:cells_suscept} Linear and nonlinear stimulus encoding in a low-CV P-unit (cell identifier ``2010-06-21-ai"). \figitem{A} Interspike interval (ISI) distribution of the cell's baseline activity, i.e. the cell is driven only by the unperturbed own electric field. The low CV of the ISIs indicates quite regular firing. \figitem{B} Power spectral density of the baseline response with peaks at the cell's baseline firing rate \fbase{} and the fish's EOD frequency \feod{}. \figitem{C} Random amplitude modulation stimulus (top, with cutoff frequency of 300\,Hz) and evoked responses (spike raster, bottom) of the same P-unit. The stimulus contrast (right) measures the strength of the AM. \figitem{D} Gain of the transfer function (first-order susceptibility), \eqnref{linearencoding_methods}, computed from the responses to 10\,\% (light purple) and 20\,\% contrast (dark purple) RAM stimulation of 10\,s duration. \figitem{E, F} Absolute value of the second-order susceptibility, \eqnref{eq:susceptibility}, for both the low and high stimulus contrast. Pink triangles mark vertical, horizontal, and diagonal lines where \fone, \ftwo{} or \fsum{} are equal to \fbase{}. \figitem{G} Second-order susceptibilities projected onto the diagonal (means of all anti-diagonals of the matrices shown in \panel{E, F}). Dots mark \fbase{}, horizontal dashed lines mark medians of the projected susceptibilities.} \end{figure*} -Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:cells_suscept}{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:cells_suscept}{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:cells_suscept}{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% \notejb{Cite Moe paper?}. +Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:cells_suscept}{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:cells_suscept}{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:cells_suscept}{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 \notejb{Cite Moe paper and Benda2005?}. -The second-order susceptibility, \Eqnref{eq:susceptibility}, quantifies the amplitude and phase of the stimulus-evoked response at the sum \fsum{} for each combination of two stimulus frequencies \fone{} and \ftwo{}. 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 simulates the neuron with a whole range of frequencies simultaneously (\subfigsref{fig:cells_suscept}{E, F}). For a LIF driven in the super-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}. 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:cells_suscept}{E, F}). +The second-order susceptibility, \eqnref{eq:susceptibility}, quantifies the amplitude and phase of the stimulus-evoked response at the sum \fsum{} for each combination of two stimulus frequencies \fone{} and \ftwo{}. 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:cells_suscept}{E, F}). For a LIF-model 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}. 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:cells_suscept}{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:cells_suscept}{E}). This structure vanishes for the stronger stimulus (\subfigref{fig:cells_suscept}{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 susceptibility values onto the diagonal by taking the mean over all anti-diagonals (\subfigrefb{fig:cells_suscept}{G}). For the low RAM contrast this projected second-order susceptibility indeed has a small peak at \fbase{} (\subfigrefb{fig:cells_suscept}{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:cells_suscept}{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. +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:cells_suscept}{E}). This structure vanishes for the stronger stimulus (\subfigref{fig:cells_suscept}{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:cells_suscept}{G}). At low RAM contrast this projected second-order susceptibility indeed has a small peak at \fbase{} (\subfigrefb{fig:cells_suscept}{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:cells_suscept}{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:cells_suscept_high_CV}). +In contrast, a high-CV P-unit (CV$_{\text{base}}=0.4$) does not exhibit pronounced nonlinearities even at low stimulus contrasts (\figrefb{fig:cells_suscept_high_CV}). 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*}[t] \includegraphics[width=\columnwidth]{ampullary.pdf} - \caption{\label{fig:ampullary} Linear and nonlinear stimulus encoding in an ampullary afferent (cell identifier ``2012-04-26-ae"). \figitem{A} Interspike interval (ISI) distribution of the cell's baseline activity. 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. \figitem{C} Bad-limited white noise stimulus (top, with a cutoff frequency of 150\,Hz) added to the fish's self-generated electric field 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 2\,\% (light green) and 20\,\% 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. Pink triangles indicate the baseline firing rate. \figitem{G} Projections of the second-order susceptibilities in \panel{E, F} onto the diagonal. } + \caption{\label{fig:ampullary} Linear and nonlinear stimulus encoding in an ampullary afferent (cell identifier ``2012-04-26-ae"). \figitem{A} Interspike interval (ISI) distribution of the cell's baseline activity. 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. \figitem{C} Bad-limited white noise stimulus (top, with a cutoff frequency of 150\,Hz) added to the fish's self-generated electric field 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 2\,\% (light green) and 20\,\% 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. Pink triangles indicate the baseline firing rate. \figitem{G} Projections of the second-order susceptibilities in \panel{E, F} onto the diagonal. } \end{figure*} -Irrespective of the CV, neither of the two example P-units shows the complete expected structure of nonlinear interactions. Electric fish possess an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogenous electric stimuli. The population of ampullary afferents is much less heterogeneous, and known for the much lower CVs of their baseline ISIs (CV$_{\text{base}}=0.06$--$0.22$) \citep{Grewe2017}. Ampullary cells do not phase-lock to the EOD and the ISIs are unimodally distributed (\subfigrefb{fig:ampullary}{A}). As a consequence of the high regularity of their baseline spiking activity, the corresponding power spectrum shows distinct peaks at the baseline firing rate \fbase{} and its harmonics. Since the cells do not fire phase-locked to the EOD, there is no peak at \feod{} (\subfigrefb{fig:ampullary}{B}). When driven by a low-contrast noise stimulus (note: this is not an AM but a stimulus that is added to the self-generated EOD, \subfigref{fig:ampullary}{C}), ampullary cells exhibit very pronounced bands in the second-order susceptibility, where \fsum{} is equal to \fbase{} or its harmonic (yellow diagonals in \subfigrefb{fig:ampullary}{E}), implying strong nonlinear response components at these frequency combinations (\subfigrefb{fig:ampullary}{G}, top). With higher stimulus contrasts these bands disappear (\subfigrefb{fig:ampullary}{F}), the projection onto the diagonal loses its distinct peak at \fsum{} and its overall level is reduced (\subfigrefb{fig:ampullary}{G}, bottom). +Electric fish possess an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogenous electric stimuli. The population of ampullary afferents is much less heterogeneous, and known for the much lower CVs of their baseline ISIs (CV$_{\text{base}}=0.06$--$0.22$) \citep{Grewe2017}. Ampullary cells do not phase-lock to the EOD and the ISIs are unimodally distributed (\subfigrefb{fig:ampullary}{A}). As a consequence of the high regularity of their baseline spiking activity, the corresponding power spectrum shows distinct peaks at the baseline firing rate \fbase{} and its harmonics. Since the cells do not 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 nonlinear structure} -In the example recordings shown above (\figsrefb{fig:cells_suscept} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagonal of the second-order susceptibility, where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, which is in line with theoretical expectations \citep{Voronenko2017}. However, a pronounced nonlinear response at frequencies \foneb{} or \ftwob{}, although predicted by theory, cannot be observed. Here we investigate how these discrepancies can be understood. +In the example recordings shown above (\figsrefb{fig:cells_suscept} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagonal of the second-order susceptibility, where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, which is in line with theoretical expectations \citep{Voronenko2017}. However, a pronounced nonlinear response at frequencies \foneb{} or \ftwob{}, although predicted by theory, cannot be observed. In the following, we investigate how these discrepancies can be understood. \begin{figure*}[t] \includegraphics[width=\columnwidth]{model_and_data.pdf} - \caption{\label{model_and_data} Estimation of second-order susceptibilities in the limit of weak stimuli. \figitem{A} \suscept{} estimated from $N=11$ 0.5\,s long trials of an electrophysiological recording of another low-CV P-unit (cell 2012-07-03-ak, $\fbase=120$\,Hz, CV=0.20) driven with a weak RAM stimulus with contrast 2.5\,\%. Pink edges mark the baseline firing rate where enhanced nonlinear responses are expected. \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 2012-07-03-ak, table~\ref{modelparams}) based on the same number of trials as in the electrophysiological recording. \figitem[iii]{B} Same as \panel[ii]{B} but using $10^6$ stimulus repetitions. \figitem[i-iii]{C} Same as in \panel[i-iii]{B} but in the \textit{noise split} condition: there is no external RAM signal 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 (\signalnoise, center), while the intrinsic noise is reduced to 10\,\% of its original strength (see methods for details). In addition to one million trials, this reveals the full expected structure of the second-order susceptibility.} + \caption{\label{model_and_data} Estimation of second-order susceptibilities in the limit of weak stimuli. \figitem{A} \suscept{} estimated from $N=11$ 0.5\,s long segments of an electrophysiological recording of another low-CV P-unit (cell 2012-07-03-ak, $\fbase=120$\,Hz, CV=0.20) driven with a weak RAM stimulus with contrast 2.5\,\%. Pink edges mark the baseline firing rate where enhanced nonlinear responses are expected. \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 2012-07-03-ak, table~\ref{modelparams}) based on the same number of trials as in the electrophysiological recording. \figitem[iii]{B} Same as \panel[ii]{B} but using $10^6$ stimulus repetitions. \figitem[i-iii]{C} Same as in \panel[i-iii]{B} but in the \textit{noise split} condition: there is no external RAM signal 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 (\signalnoise, center), while the intrinsic noise is reduced to 10\,\% of its original strength (see methods for details). In addition to one million trials, this reveals the full expected structure of the second-order susceptibility.} \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 reason could be simply too little data for a good estimate of the second-order susceptibility. Electrophysiological recordings are limited in time, and therefore responses to only a limited number of trials, i.e. repetitions of the same RAM stimulus, are available. As a consequence, the cross-spectra, \Eqnref{eq:crosshigh}, are insufficiently averaged and the full structure of the second-order susceptibility might be hidden in finite-data noise. This experimental limitation can be overcome by using a computational model for the P-unit, a stochastic leaky integrate-and-fire model with adaptation current and dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of a low-CV cell estimated from the same low number of trials as in the experiment ($\n{}=11$, compare \panel{A} and \panel[ii]{B} in \figrefb{model_and_data}). +One reason could be simply too little data for a good estimate of the second-order susceptibility. Electrophysiological recordings are limited in time, and therefore responses to only a limited number of trials, i.e. repetitions of the same RAM stimulus, are available. As a consequence, the cross-spectra, \eqnref{eq:crosshigh}, are insufficiently averaged and the full structure of the second-order susceptibility might be hidden in finite-data noise. This experimental limitation can be overcome by using a computational model for the P-unit, a stochastic leaky integrate-and-fire model with adaptation current, dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of a low-CV cell estimated from the same low number of trials as in the experiment ($\n{}=11$, compare \panel{A} and \panel[ii]{B} in \figrefb{model_and_data}). -In simulations of the model, we can increase the number of trials beyond what would be experimentally possible, here to one million (\subfigrefb{model_and_data}\,\panel[iii]{B}). The estimate of the second-order susceptibility indeed improves. It gets less noisy and the diagonal at \fsum{} is emphasized. However, the expected vertical and horizontal lines at \foneb{} and \ftwob{} are still mainly missing. +In model simulations we can increase the number of trials beyond what would be experimentally possible, here to one million (\subfigrefb{model_and_data}\,\panel[iii]{B}). The estimate of the second-order susceptibility indeed improves. It gets less noisy and the diagonal at \fsum{} is emphasized. However, the expected vertical and horizontal lines at \foneb{} and \ftwob{} are still mainly missing. -Using a broadband stimulus increases the effective input-noise level and this may linearize signal transmission and suppress potential nonlinear responses \citep{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017}. Assuming that the intrinsic noise level in this P-unit is small enough, the full expected structure of the second-order susceptibility should appear in the limit of weak AMs. Again, this cannot be done experimentally, because the problem of insufficient averaging becomes even more severe for weak AMs (low contrast). In the model, however, we know the time course of the intrinsic noise and can use this knowledge to determine the susceptibilities by input-output correlations via the Furutsu-Novikov theorem \citep{Furutsu1963, Novikov1965}. This theorem, in its simplest form, states that the cross-spectrum $S_{x\eta}(\omega)$ of a Gaussian noise $\eta(t)$ driving a nonlinear system and the system's output $x(t)$ is proportional to the linear susceptibility according to $S_{x\eta}(\omega)=\chi(\omega)S_{\eta\eta}(\omega)$. Here $\chi(\omega)$ characterizes the linear response to an infinitely weak signal $s(t)$ in the presence of the background noise $\eta(t)$. Likewise, the nonlinear susceptibility can be determined in an analogous fashion from higher-order input-output cross-spectra (see methods, equations \eqref{eq:crosshigh} and \eqref{eq:susceptibility}) \citep{Egerland2020}. In line with an alternative derivation of the Furutsu-Novikov theorem \citep{Lindner2022}, we can split the total noise and consider a fraction of it as a stimulus. This allows us to calculate the susceptibility from the cross-spectrum between the output and this stimulus fraction of the noise. Adapting this approach to our P-unit model (see methods), we replace the intrinsic noise by an approximately equivalent RAM stimulus $s_\xi(t)$ and a weak remaining intrinsic noise $\sqrt{2D \, c_{\rm{noise}}} \cdot \xi(t)$ with $c_\text{noise} = 0.1$ (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \subfigrefb{model_and_data}\,\panel[i]{C}). We tune the amplitude of the RAM stimulus $s_\xi(t)$ such that the output firing rate and variability (CV) are the same as in the baseline activity (i.e. full intrinsic noise $\sqrt{2D}\xi(t)$ in the voltage equation but no RAM) and compute the cross-spectra between the RAM part of the noise $s_\xi(t)$ and the output spike train. This procedure has two consequences: (i) by means of the cross-spectrum between the output and \signalnoise, which is a large fraction of the noise, the signal-to-noise ratio of the measured susceptibilities is drastically improved; (ii) the total noise in the system has been reduced (by what was before the external RAM stimulus $s(t)$), which makes the system more nonlinear. For both reasons we now see the expected nonlinear features in the second-order susceptibility for a sufficient number of trials (\subfigrefb{model_and_data}\,\panel[iii]{C}), but not for a number of trials comparable to the experiment (\subfigrefb{model_and_data}\,\panel[ii]{C}). In addition to the strong response for \fsumb{}, we now also observe pronounced nonlinear responses at \foneb{} and \ftwob{} (vertical and horizontal lines, \subfigrefb{model_and_data}\,\panel[iii]{C}). +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{model_and_data}\,\panel[i]{C}). We tune the amplitude of the RAM stimulus $s_{\xi}(t)$ such that the output firing rate and variability (CV) are the same as in the baseline activity (i.e. full intrinsic noise $\sqrt{2D}\;\xi(t)$ in the voltage equation but no RAM) and compute the cross-spectra between the RAM part of the noise $s_{\xi}(t)$ and the output spike train. This procedure has two consequences: (i) by means of the cross-spectrum between the output and \signalnoise, which is a large fraction of the noise, the signal-to-noise ratio of the measured susceptibilities is drastically improved; (ii) the total noise in the system has been reduced (by what was before the external RAM stimulus $s(t)$), which makes the system more nonlinear. For both reasons we now see the expected nonlinear features in the second-order susceptibility for a sufficient number of trials (\subfigrefb{model_and_data}\,\panel[iii]{C}), but not for a number of trials comparable to the experiment (\subfigrefb{model_and_data}\,\panel[ii]{C}). In addition to the strong response for \fsumb{}, we now also observe pronounced nonlinear responses at \foneb{} and \ftwob{} (vertical and horizontal lines, \subfigrefb{model_and_data}\,\panel[iii]{C}). Note, that the increased number of trials goes along with a substantial reduction of second-order susceptibility values (\subfigrefb{model_and_data}\,\panel[iii]{C}), that saturate in its peak values for $N>10^6$ (\figrefb{fig:trialnr}). This demonstrates the limited reliability of an estimate of the second-order susceptibility that is based on 11 trials only. However, we would like to point out that already the limited number of trials used in the experiments reveals key features of the nonlinear response. @@ -517,30 +508,30 @@ With high levels of intrinsic noise, we would not expect the nonlinear response \subsection{Second-order susceptibility can explain nonlinear peaks in pure sinewave stimulation} -We estimated the second-order susceptibility of P-unit responses using RAM stimuli. In particular, we found pronounced nonlinear responses in the limit of weak stimulus amplitudes. How do these findings relate to the situation of two 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 \bone{} and \btwo{} 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, $\Delta f_{2}$ was similar to \fbase{}, corresponding to the horizontal line of the second-order susceptibility estimated for a vanishing external RAM stimulus (\subfigrefb{model_and_data}\,\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{model_and_data}\,\panel[iii]{C}), in which only the response at the sum of the two stimulus frequencies is addressed. % less prominent, +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 \bone{} and \btwo{} 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, $\Delta f_{2}$ was similar to \fbase{}, corresponding to the horizontal line of the second-order susceptibility estimated for a vanishing external RAM stimulus (\subfigrefb{model_and_data}\,\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{model_and_data}\,\panel[iii]{C}), in which only the response at the sum of the two stimulus frequencies is addressed. % less prominent, \begin{figure*}[t] \includegraphics[width=\columnwidth]{model_full.pdf} - \caption{\label{fig:model_full} Using second-order susceptibility to predict responses to sine-wave stimuli. \figitem[]{A} Absolute value of the second-order susceptibility, \Eqnref{eq:susceptibility}, for both positive and negative frequencies. \susceptf{} was estimated from $N=10^6$ trials of model simulations in the noise-split condition (cell 2013-01-08-aa, see table~\ref{modelparams} for model parameters). White lines indicate zero frequencies. Nonlinear responses at \fsum{} are quantified in the upper right and lower left quadrants. Nonlinear responses at \fdiff{} are quantified in the upper left and lower right quadrants. The baseline firing rate of this cell was at $\fbase=120$\,Hz. The position of the orange/red letters corresponds to the beat frequencies used for the stimulation with pure sine waves in the subsequent panels and indicates the sum/difference of those beat frequencies. \figitem{B--E} Black line -- power spectral density of model simulations in response to stimulation with two pure sine waves, \fone{} and \ftwo, in addition to the receiving fish's own EOD (three-fish scenario). The contrast of beat beats is 0.02. Colored circles highlight the height of selected peaks in the power spectrum. Grey line -- power spectral density of model in the baseline condition. \figitem{B} The sum of the two beat frequencies match \fbase{}. \figitem{C} The difference of \fone{} and \ftwo{} match \fbase{}. \figitem{D} Only the first beat frequency matches \fbase{}. \figitem{E} None of the two beat frequencies matches \fbase{}.} + \caption{\label{fig:model_full} Using second-order susceptibility to predict responses to sine-wave stimuli. \figitem[]{A} Absolute value of the second-order susceptibility, \eqnref{eq:susceptibility}, for both positive and negative frequencies. \susceptf{} was estimated from $N=10^6$ trials of model simulations in the noise-split condition (cell 2013-01-08-aa, see table~\ref{modelparams} for model parameters). White lines indicate zero frequencies. Nonlinear responses at \fsum{} are quantified in the upper right and lower left quadrants. Nonlinear responses at \fdiff{} are quantified in the upper left and lower right quadrants. The baseline firing rate of this cell was at $\fbase=120$\,Hz. The position of the orange/red letters corresponds to the beat frequencies used for the stimulation with pure sine waves in the subsequent panels and indicates the sum/difference of those beat frequencies. \figitem{B--E} Black line -- power spectral density of model simulations in response to stimulation with two pure sine waves, \fone{} and \ftwo, in addition to the receiving fish's own EOD (three-fish scenario). The contrast of beat beats is 0.02. Colored circles highlight the height of selected peaks in the power spectrum. Grey line -- power spectral density of model in the baseline condition. \figitem{B} The sum of the two beat frequencies match \fbase{}. \figitem{C} The difference of \fone{} and \ftwo{} match \fbase{}. \figitem{D} Only the first beat frequency matches \fbase{}. \figitem{E} None of the two beat frequencies matches \fbase{}.} \end{figure*} -However, the second-order susceptibility \Eqnref{eq:susceptibility} is a spectral measure that is based on Fourier transforms and thus is also defined for negative stimulus frequencies. The full \susceptf{} matrix is symmetric with respect to the origin. In the upper-right and lower-left quadrants of \susceptf{}, stimulus-evoked responses at \fsum{} are shown, whereas in the lower-right and upper-left quadrants nonlinear responses at the difference \fdiff{} are shown (\figref{fig:model_full}). The vertical and horizontal lines at \foneb{} and \ftwob{} are very pronounced in the upper-right quadrant of \subfigrefb{fig:model_full}{A} for nonlinear responses at \fsum{} and extend into the upper-left quadrant (representing \fdiff) where they fade out towards more negative $f_1$ frequencies. The peak in the response power-spectrum at \fdiff{} evoked by pure sine-wave stimulation (\subfigrefb{fig:motivation}{D}) is predicted by the horizontal line in the upper-left quadrant (\subfigrefb{fig:model_full}{A}, \citealp{Schlungbaum2023}). -Is it possible based on the second-order susceptibility estimated by means of RAM stimuli (\subfigrefb{fig:model_full}{A}) to predict nonlinear responses in a three-fish setting? We can test this by stimulating the same model with two beats with weak amplitudes (\subfigrefb{fig:model_full}{B--E}). If we choose a frequency combination where the sum of the two beat frequencies is equal to the model's baseline firing rate \fbase{}, a peak at the sum of the two beat frequencies appears in the power spectrum of the response (\subfigrefb{fig:model_full}{B}), as expected from \suscept. If instead we choose two beat frequencies that differ by \fbase{}, a peak is present at the difference frequency (\subfigrefb{fig:model_full}{C}). If only one beat frequency is equal to \fbase{}, both a peak at the sum and at the difference frequency are present in the P-unit response (\subfigrefb{fig:model_full}{D}). And if none of these conditions are met, neither a peak at the sum nor at the difference of the two beat frequencies appears (\subfigrefb{fig:model_full}{E}). +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}). \begin{figure*}[tp] \includegraphics[width=\columnwidth]{data_overview_mod.pdf} - \caption{\label{fig:data_overview} Nonlinear responses in P-units and ampullary cells. The second-order susceptibility is condensed into the peakedness of the nonlinearity, \nli{} \Eqnref{eq:nli_equation}, that relates the amplitude of the projected susceptibility at a cell's baseline firing rate to its median (see \subfigrefb{fig:cells_suscept}{G}). Each of the recorded neurons contributes at maximum with two stimulus contrasts. Black squares and circles highlight recordings conducted in a single cell. Squares in \panel{A, C, E} correspond to the cell in \figrefb{fig:cells_suscept} and circles to the cell in \figrefb{fig:cells_suscept_high_CV}. Squares in \panel{B, D, F} correspond to the cell in \figrefb{fig:ampullary}. \figitem{A, B} There is a negative correlation between the CV during baseline and \nli. \figitem{C, D} There is a negative correlation between the CV during stimulation and \nli. \figitem{E, F} \nli{} is plotted against the response modulation, (see methods), an indicator of the subjective stimulus strength for a cell. There is a negative correlation between response modulation and \nli. Restricting the analysis to the weakest stimulus that was presented to each unique neuron, does not change the results. The number of unique neurons is 221 for P-units and 45 for ampullary cells. + \caption{\label{fig:data_overview} Nonlinear responses in P-units and ampullary cells. The second-order susceptibility is condensed into the peakedness of the nonlinearity, \nli{} \eqnref{eq:nli_equation}, that relates the amplitude of the projected susceptibility at a cell's baseline firing rate to its median (see \subfigrefb{fig:cells_suscept}{G}). Each of the recorded neurons contributes at maximum with two stimulus contrasts. Black squares and circles highlight recordings conducted in a single cell. Squares in \panel{A, C, E} correspond to the cell in \figrefb{fig:cells_suscept} and circles to the cell in \figrefb{fig:cells_suscept_high_CV}. Squares in \panel{B, D, F} correspond to the cell in \figrefb{fig:ampullary}. \figitem{A, B} There is a negative correlation between the CV during baseline and \nli. \figitem{C, D} There is a negative correlation between the CV during stimulation and \nli. \figitem{E, F} \nli{} is plotted against the response modulation, (see methods), an indicator of the subjective stimulus strength for a cell. There is a negative correlation between response modulation and \nli. Restricting the analysis to the weakest stimulus that was presented to each unique neuron, does not change the results. The number of unique neurons is 221 for P-units and 45 for ampullary cells. % The two example P-units shown before are highlighted with dark markers in \subfigrefb{fig:data_overview}{A, C, E} (squares -- \figrefb{fig:cells_suscept}, circles -- \figrefb{fig:cells_suscept_high_CV}). % Several of the recorded neurons contribute with two samples to the population analysis as their responses have been recorded to two different contrasts of the same RAM stimulus. Higher stimulus contrasts lead to a stronger drive and thus stronger response modulations (see color code bar in \subfigref{fig:data_overview}{A}, see methods). % The example cell shown above (\figref{fig:ampullary}) was recorded at two different stimulus intensities and the \nli{} values are highlighted with black squares. } \end{figure*} -%\Eqnref{response_modulation} \subsection{Low CVs and weak stimuli are associated with strong nonlinearity} -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 specimens. For comparison across cells we summarize the structure of the second-order susceptibilities in a single number, the peakedness of the nonlinearity \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:cells_suscept}{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. The \nli{} values of the P-unit population depend weakly on the CV of the baseline ISI distribution. Cells with lower baseline CVs tend to have more pronounced peaks in their projections than those that have high baseline CVs (\subfigrefb{fig:data_overview}{A}). This negative correlation is more pronounced against the CV measured during stimulation (\subfigrefb{fig:data_overview}{C}). +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. 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:cells_suscept}{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. The \nli{} values of the P-unit population depend weakly on the CV of the baseline ISI distribution. Cells with lower baseline CVs tend to have more pronounced peaks in their projections than those that have high baseline CVs (\subfigrefb{fig:data_overview}{A}). This negative correlation is more pronounced against the CV measured during stimulation (\subfigrefb{fig:data_overview}{C}). -The effective stimulus strength also plays an important role. We quantify the effect a stimulus has 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:data_overview}{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). +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:data_overview}{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$ @@ -571,14 +562,15 @@ The population of ampullary cells is generally more homogeneous, with lower base %Even though the second-order susceptibilities here were estimated from data and models with a modulated (EOD) carrier (\figrefb{fig:model_full}) they are in good accordance with the second-order susceptibilities found in LIF models without a carrier \citep{Voronenko2017, Schlungbaum2023}. -%\,\panel[iii]{C} +\notejb{Insert short summary here} + \subsection{Theory applies to systems with and without carrier} -Theoretical work \citep{Voronenko2017} explained analytically the occurrence of nonlinear products when a LIF model neuron is stimulated with pure sine waves. To investigate whether the same mechanisms occur in electroreceptor afferents which are driven by AMs of a carrier and not by pure sine-waves, we followed the previous approach and quantified the second-order susceptibility from responses to white-noise stimuli \citep{Voronenko2017, Egerland2020, Neiman2011fish, Nikias1993}. We expected to see elevated second-order susceptibility where either of the foreign signals matches the baseline firing rate ($f_1=\fbase{}$ or $f_2=\fbase{}$) or when the sum equals the baseline firing rate of the neuron (\fsumb{}) creating a triangular pattern of elevated \suscept{} e.g.\,\subfigref{model_and_data}\,\panel[iii]{C}. Indeed, we find traces of the same nonlinearities in the neuronal responses of p-type electroreceptor afferents. The nonlinear pattern observed in the experimental data, however, matches the expectations only partially and only in a subset of neurons (\figsref{fig:cells_suscept} and \ref{fig:ampullary}). Nevertheless, the theory holds also for systems that are driven by AMs of a carrier and is thus more widely applicable. +Theoretical work \citep{Voronenko2017} explained analytically the occurrence of nonlinear products when a LIF model neuron is stimulated with pure sine waves. To investigate whether the same mechanisms occur in electroreceptor afferents which are driven by AMs of a carrier and not by pure sine-waves, we followed the previous approach and quantified the second-order susceptibility from responses to white-noise stimuli \citep{Voronenko2017, Egerland2020, Neiman2011fish, Nikias1993}. We expected to see elevated second-order susceptibility where either of the foreign signals matches the baseline firing rate ($f_1=\fbase{}$ or $f_2=\fbase{}$) or when the sum equals the baseline firing rate of the neuron (\fsumb{}) creating a triangular pattern of elevated \suscept{} (e.g.\,\subfigrefb{model_and_data}\,\panel[iii]{C}). Indeed, we find traces of the same nonlinearities in the neuronal responses of p-type electroreceptor afferents. The nonlinear pattern observed in the experimental data, however, matches the expectations only partially and only in a subset of neurons (\figsref{fig:cells_suscept} and \ref{fig:ampullary}). Nevertheless, the theory holds also for systems that are driven by AMs of a carrier and is thus more widely applicable. \subsection{Intrinsic noise limits nonlinear responses} -Only those P-units that exhibit low coefficients of variation (CV) of the interspike-interval distribution (\subfigref{fig:cells_suscept}{A}) in their unperturbed baseline response show the expected nonlinearities (\subfigref{fig:data_overview}{A}). Such low-CV cells are rare among the 221 P-units that we used in this study. The afferents of the passive electrosensory system, the ampullary cells, however, have generally lower CVs and show a much clearer nonlinearity pattern than the low-CV P-unit exemplified here (compare \figsref{fig:cells_suscept} and \ref{fig:ampullary}). The single ampullary cell featured in \figref{fig:ampullary} is representative of the majority of ampullary cells analyzed here. All ampullary cells have CVs below 0.4 with a median around 0.12 and the observed \nli{}s are 10-fold higher than in P-units. +Only those P-units that exhibit low coefficients of variation (CV) of the interspike-interval distribution (\subfigref{fig:cells_suscept}{A}) in their unperturbed baseline response show the expected nonlinearities (\subfigref{fig:data_overview}{A}). Such low-CV cells are rare among the 221 P-units that we used in this study. The afferents of the passive electrosensory system, the ampullary cells, however, have generally lower CVs and show a much clearer nonlinearity pattern than the low-CV P-unit exemplified here (compare \figsref{fig:cells_suscept} and \ref{fig:ampullary}). The single ampullary cell featured in \figref{fig:ampullary} is representative of the majority of ampullary cells analyzed here. All ampullary cells have CVs below 0.4 with a median around 0.12 and the observed \nli{}s are about 10-fold higher than in P-units. -The CV serves as a proxy for the intrinsic noise in the cells. In both cell types, we observe a negative correlation between \nli{} and the CV, indicating that it is the level of intrinsic noise that plays a role here. These findings are in line with previous studies that propose that noise linearizes the system \citep{Roddey2000, Chialvo1997, Voronenko2017}. More intrinsic noise has been demonstrated to increase the CV and 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}. Further support for the notion of noise limiting the nonlinearity comes from our P-unit LIF model that faithfully reproduces P-unit activity \citep{Barayeu2023}. We can use this model and apply a noise-split \citep{Lindner2022} based on the Furutsu-Novikov theorem \citep{Novikov1965, Furutsu1963}, to increase the signal-to-noise ratio in the cell while keeping the overall response variability constant (see methods). Treating 90\,\% of the total noise as a signal and simulating large numbers of trials uncovers the full nonlinear structure (\figref{model_and_data}) seen in LIF neurons and the analytical derivations when driven with sine-wave stimuli \citep{Voronenko2017}. +The CV is a proxy for the intrinsic noise in the cells \notejb{Cite Lindner IF by rate and CV}. In both cell types, we observe a negative correlation between \nli{} and the CV, indicating that it is the level of intrinsic noise that plays a role here. These findings are in line with previous studies that propose that noise linearizes the system \citep{Roddey2000, Chialvo1997, Voronenko2017}. More intrinsic noise has been demonstrated to increase the CV and 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}. Further support for the notion of noise limiting the nonlinearity comes from our P-unit LIF model that faithfully reproduces P-unit activity \citep{Barayeu2023}. We can use this model and apply a noise-split \citep{Lindner2022} based on the Furutsu-Novikov theorem \citep{Novikov1965, Furutsu1963}, to increase the signal-to-noise ratio in the cell while keeping the overall response variability constant (see methods). Treating 90\,\% of the total noise as a signal and simulating large numbers of trials uncovers the full nonlinear structure (\figref{model_and_data}) seen in LIF neurons and the analytical derivations when driven with sine-wave stimuli \citep{Voronenko2017}. % \subsection{Noise stimulation approximates the real three-fish interaction} @@ -618,274 +610,217 @@ auditory nerve fibers and such nonlinear effects might also be expected in the a \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{}. The original electrophysiological recordings were performed on male and female weakly electric fish of the species \lepto{} that were obtained from a commercial supplier for tropical fish (Aquarium Glaser GmbH, Rodgau, -Germany). The fish were 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. +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 with 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). +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 reapplication 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 (\figrefb{fig:setup}, blue triangle). 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. +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 and a clear phase-locking to the EOD and, 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. +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} - The electric field of the fish was recorded in two ways: 1. we measured the so-called \textit{global EOD} with two vertical carbon rods ($11\,\centi\meter$ long, 8\,mm diameter) in a head-tail configuration (\figrefb{fig:setup}, green bars). The electrodes were placed isopotential to the stimulus. This signal was differentially amplified with a 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). 2. The so-called \textit{local EOD} was measured with 1\,cm-spaced silver wires 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, \figrefb{fig:setup}, red markers). This local measurement recorded the combination of the fish's own field and the applied stimulus and thus serves as a proxy of the transdermal potential that drives the electroreceptors. +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} -The stimulus was isolated from the 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 (\figrefb{fig:setup}, gray bars). The stimulus was calibrated with respect to the local EOD. -\begin{figure*}[t] - \includegraphics[width=\columnwidth]{Settup} - \caption{\label{fig:setup} Electrophysiolocical recording setup. The fish, depicted as a black scheme and surrounded by isopotential lines, was positioned in the center of the tank. Blue triangle -- electrophysiological recordings were conducted in the posterior anterior lateral line nerve (pALLN). Gray horizontal bars -- electrodes for the stimulation. Green vertical bars -- electrodes to measure the \textit{global EOD} placed isopotential to the stimulus, i.e. recording fish's unperturbed EOD. Red dots -- electrodes to measure the \textit{local EOD} picking up the combination of fish's EOD and the stimulus. The local EOD was measured with a distance of 1 \,cm between the electrodes. All measured signals were amplified, filtered, and stored for offline analysis.} -\end{figure*} +\subsection{Stimulation}\label{rammethods} +ELectric stimuli were isolated from the 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. \notejb{attenuator?} -\subsection{White noise stimulation}\label{rammethods} -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 cell recordings 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). +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 performed with Python 3 using the packages matplotlib \citep{Hunter2007}, numpy \citep{Walt2011}, scipy \citep{scipy2020}, pandas \citep{Mckinney2010}, nixio \citep{Stoewer2014}, and thunderfish (\url{https://github.com/bendalab/thunderfish}). +\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 thunderfish (\url{https://github.com/bendalab/thunderfish}). %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) was calculated as the standard deviation of the interspike intervals (ISI) divided by 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. +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. -The time-dependent response of the neuron was estimated from the spiking activity -\begin{equation}\label{eq:spikes} -x_k(t) = \sum_i\delta(t-t_{k,i}) +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} -recorded for each stimulus presentation, $k$, by kernel convolution with a Gaussian kernel +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} -K(t) = \scriptstyle \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{t^2}{2\sigma^2}} + r_k(t) = x_k(t) * K(t) = \int_{-\infty}^{+\infty} x_k(t') K(t-t') \, \mathrm{d}t' \end{equation} -with $\sigma$ the standard deviation of the Gaussian which was set to 2.5\,ms if not stated otherwise. For each trial $k$ the $x_k(t)$ is convolved with the kernel $K(t)$ - +was estimated by convolving the spike train with a kernel. We used a Gaussian kernel \begin{equation} - r_k(t) = x_k(t) * K(t) = \int_{-\infty}^{+\infty} x_k(t') K(t-t') \, \mathrm{d}t' \;, +K(t) = {\scriptstyle \frac{1}{\sigma\sqrt{2\pi}}} e^{-\frac{t^2}{2\sigma^2}} \end{equation} -where $*$ denotes the convolution. $r(t)$ is then calculated as the across-trial average +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} - r(t) = \left\langle r_k(t) \right\rangle _k. + \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} -To quantify how strongly the neuron is driven by the stimulus we quantified the response modulation as the standard deviation $\sigma_{M} = \sqrt{\langle (r(t)-\langle r(t) \rangle_t)^2\rangle_t}$, where $\langle \cdot \rangle_t$ indicates averaging over time. +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 the spiking response $x(t)$, \Eqnref{eq:spikes}, depends on the stimulus $s(t)$. To investigate the relation between stimulus and response we calculated the first- and second-order susceptibility of the neuron to the stimulus in the frequency domain. The Fourier transforms of $s(t)$ and $x(t)$ are denoted as $\tilde s(\omega)$ and $\tilde x(\omega)$ and were calculated according to $\tilde x(\omega) = \int_{0}^{T} \, x(t) \cdot e^{- i \omega t}\,dt$, with $T$ being the signal duration. Stimuli had a duration of 10\,s and spectra of stimulus and response were calculated in separate segments of 0.5\,s with no overlap resulting in a spectral resolution of 2\,Hz. +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. + +Fourier transforms $\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 simplicity we use angular frequencies $\omega=2\pi f$ instead of frequencies $f$. 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}. -The power spectrum of the stimulus $s(t)$ was calculated as +The power spectrum of the stimulus $s(t)$ was estimated as \begin{equation} \label{powereq} - \begin{split} - S_{ss}(\omega) = \frac{\langle \tilde s(\omega) \tilde s^* (\omega)\rangle}{T} - \end{split} + S_{ss}(\omega) = \frac{\langle \tilde s_T(\omega) \tilde s_T^*(\omega)\rangle}{T} \end{equation} -with $\tilde s^* $ being the complex conjugate and $\langle ... \rangle$ denoting averaging over the segments. The power spectrum of the spike trains $S_{xx}(\omega)$ was calculated accordingly. The cross-spectrum $S_{xs}(\omega)$ between stimulus and evoked spike trains was calculated according to +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} - \begin{split} - S_{xs}(\omega) = \frac{\langle \tilde x(\omega) \tilde s^* (\omega)\rangle}{T} - \end{split} + S_{xs}(\omega) = \frac{\langle \tilde x_T(\omega) \tilde s_T^*(\omega)\rangle}{T} \end{equation} -From $S_{xs}(\omega)$ and $ S_{ss}(\omega)$ we calculated the linear susceptibility (transfer function) as + +The first-order susceptibility (transfer function) \begin{equation} \label{linearencoding_methods} - \begin{split} - \chi_{1}(\omega) = \frac{S_{xs}(\omega) }{S_{ss}(\omega) } - \end{split} + \chi_{1}(\omega) = \frac{S_{xs}(\omega) }{S_{ss}(\omega) } \end{equation} -The second-order cross-spectrum that depends on the two frequencies $\omega_1$ and $\omega_2$ was calculated according to +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 (\omega_{1}+\omega_{2}) \tilde s^* (\omega_{1})\tilde s^* (\omega_{2}) \rangle}{T} + 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} -The second-order susceptibility was calculated by dividing the higher-order cross-spectrum by the spectral power at the respective frequencies. +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} - %\begin{split} \chi_{2}(\omega_{1}, \omega_{2}) = \frac{S_{xss} (\omega_{1},\omega_{2})}{2S_{ss} (\omega_{1}) S_{ss} (\omega_{2})} - %\end{split} \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} -% \begin{split} % \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{split} % \end{equation} -The absolute value of a second-order susceptibility matrix is visualized in \figrefb{fig:model_full}. There the upper right and the lower left quadrants characterize the nonlinearity in the response $x(t)$ at the sum frequency of the two input frequencies. The lower right and upper left quadrants characterize the nonlinearity in the response $x(t)$ at the difference of the input frequencies. +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 expect to see nonlinear susceptibility when $\omega_1 + \omega_2 = \fbase{}$. To characterize this we calculated the peakedness of the nonlinearity (PNL) as - \begin{equation} - \label{eq:nli_equation} - \nli{} = \frac{ \max D(\fbase{}-5\,\rm{Hz} \leq f \leq \fbase{}+5\,\rm{Hz})}{\mathrm{med}(D(f))} - \end{equation} -For this index, the second-order susceptibility matrix was projected onto the diagonal $D(f)$, by taking the mean of the anti-diagonals. The peakedness at the frequency $\fbase{}$ in $D(f)$ was quantified by finding the maximum of $D(f)$ in the range $\fbase{} \pm 5$\,Hz (\subfigrefb{fig:cells_suscept}{G}) and dividing it by the median of $D(f)$. +We expected to see elevated values in the second-order susceptibility at $\omega_1 + \omega_2 = \fbase$ \citep{Voronenko2017}. To characterize this in a single number we computed the peakedness of the nonlinearity (PNL) defined as +\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} +where $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:cells_suscept}{G}). + +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:data_overview} the mean of the resulting \nli{} values is used. -If the same frozen noise was recorded several times in a cell, each noise repetition resulted in a separate second-order susceptibility matrix. The mean of the corresponding \nli{} values was used for the population statistics in \figref{fig:data_overview}. + +\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}\label{lifmethods} -Leaky integrate-and-fire (LIF) models with a carrier 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 +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 the EOD frequency $f_{EOD}$ and an amplitude normalized to one. +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$. -In the model, the input \carrierinput{} was then first thresholded to model the synapse between the primary receptor cells and the afferent. +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} +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. + +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} -$\lfloor \carrierinput \rfloor_{0}$ denotes the threshold operation that sets negative values to zero (\subfigrefb{flowchart}{A}). - -The resulting receptor signal was then low-pass filtered to approximate passive signal conduction in the afferent's dendrite (\subfigrefb{flowchart}{B}) +(\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} -with $\tau_{d}$ as the dendritic time constant. Dendritic low-pass filtering was necessary to reproduce the loose coupling of P-unit spikes to the EOD while maintaining high sensitivity at small amplitude modulations. Because the input was dimensionless, the dendritic voltage was dimensionless, too. The combination of threshold and low-pass filtering extracts the amplitude modulation of the input \carrierinput. +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)$ was the input to a leaky integrate-and-fire (LIF) model +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) + \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$. All state variables except $\tau_m$ are dimensionless. +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$ followed +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$. +with adaptation time constant $\tau_A$ \citep{Benda2003,Benda2005,Benda2010}. -Whenever the membrane voltage $V_m(t)$ crossed 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}). +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 static nonlinearity $f(V_m)$ was equal to zero for the LIF. In the case of an exponential integrate-and-fire model (EIF), this function was set to -% \begin{equation} -% \label{eifnl} -% f(V_m)= \Delta_V \text{e}^{\frac{V_m-1}{\Delta_V}} -% \end{equation} -% \citep{Fourcaud-Trocme2003}, where $\Delta_V$ was varied from 0.001 to 0.1. -%, \figrefb{eif} - -\begin{figure*}[t] - \includegraphics[width=\columnwidth]{flowchart.pdf} - \caption{\label{flowchart} - Components of the P-unit model. The main steps of the model are illustrated in the left column. The three other columns show the corresponding signals in three different settings: (i) the baseline situation, no external stimulus, only the fish's self-generated EOD (i.e. the carrier) is present. (ii) RAM stimulation, the carrier is amplitude modulated with a weak (2\,\% contrast) band-limited white-noise stimulus. (iii) Noise split condition in which 90\,\% of the internal noise is used as a driving RAM stimulus scaled with the correction factor $\rho$ (see text). Note that the mean firing rate and the CV of the ISI distribution is the same in this and the baseline condition. As an example, simulations of the model for cell ``2012-07-03-ak'' are shown (see table~\ref{modelparams} for model parameters). \figitem{A} Thresholding: a simple linear threshold was applied to the EOD carrier, \Eqnref{eq:eod}. \figitem{B} Subsequent dendritic low-pass filtering attenuates the carrier and carves out the AM signal. \figitem{C} Gaussian white-noise is added to the signal in \panel{B}. Note the reduced internal noise amplitude in the noise split (iii) condition. \figitem{D} Spiking output of the LIF model in response to the sum of B and C. \figitem{E} Power spectra of the LIF neuron's spiking activity. Under the baseline condition (\panel[i]{E}) there are several peaks, from left to right, at the baseline firing rate $\fbase{}$, $f_{EOD} - \fbase{}$, $f_{EOD}$, and $f_{EOD} + \fbase{}$. In the stimulus-driven regime (\panel[ii]{E}), there is only a peak at \feod, while under the noise split condition (\panel[iii]{E}) the same peaks as in the baseline condition are present.} -\end{figure*} - -\subsection{Numerical implementation} -The model's ODEs 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}$: -\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} - -\subsection{Model parameters}\label{paramtext} -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-like increases and decreases in EOD amplitude (onset-state and steady-state responses, effective adaptation time constant). For each simulation, the start parameters $A$, $V_{d}$ and $V_{m}$ were drawn from a random starting value distribution, estimated from a 100\,s baseline simulation after an initial 100\,s of simulation that was discarded as a transient. - -\subsection{Stimuli for the model} -The model neurons were driven with similar stimuli as the real neurons in the original experiments. To mimic the interaction with one or two foreign animals the receiving fish's EOD (\Eqnref{eq:eod}) was normalized to an amplitude of one and the respective EODs of a second or third fish were added.%\Eqnref{ eq.\,\ref{eq:eod} - -The random amplitude modulation (RAM) input to the model was created by drawing random amplitude and phases from Gaussian distributions for each frequency component in the range 0--300 Hz. An inverse Fourier transform was applied to get the final amplitude RAM time course. The input to the model was then -\begin{equation} - \label{eq:ram_equation} - y(t) = (1+ s(t)) \cdot \cos(2\pi f_{EOD} t) -\end{equation} -From each simulation 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. -% \subsection{Second-order susceptibility analysis of the model} -% %\subsubsection{Model second-order nonlinearity} - -% The second-order susceptibility in the model was calculated with \Eqnref{eq:susceptibility}, resulting in matrices as in \figrefb{model_and_data} and \figrefb{fig:model_full}. For this, the model neuron was presented the input $x(t)$ for 2\,s, with the first second being dismissed as the transient. The second-order susceptibility calculation was performed on the last second, resulting in a frequency resolution of 1\,Hz. - -\subsection{Model noise split into a noise and a stimulus component}\label{intrinsicsplit_methods}% the Furutsu-Novikov Theorem with the same correlation function -According to previous works \citep{Lindner2022} the total noise of a LIF model ($\xi$) can be split up into several independent noise processes. Here we split the internal noise into two parts: (i) One part is treated as a driving input signal $s_\xi(t)$, a RAM stimulus where frequencies above 300\,Hz are discarded (\Eqnref{eq:ram_split}), and used to calculate the cross-spectra in \Eqnref{eq:crosshigh} and (ii) the remaining noise $\sqrt{2D \, c_{\rm{noise}}} \cdot \xi(t)$ that is treated as pure noise (\Eqnref{eq:Noise_split_intrinsic}). In this way, the effective signal-to-noise ratio can be increased while maintaining the total noise in the system. -%\sqrt{\rho \, 2D \,c_{\rm{signal}}} \cdot \xi(t) - -%(1-c_{\rm{signal}})\cdot\xi$c_{\rm{noise}} = 1-c_{\rm{signal}}$ -%c_{\rm{signal}} \cdot \xi -\begin{equation} - \label{eq:ram_split} - y(t) = (1+ s_\xi(t)) \cdot \cos(2\pi f_{EOD} t) -\end{equation} - -\begin{equation} - \label{eq:Noise_split_intrinsic_dendrite} - \tau_{d} \frac{d V_{d}}{d t} = -V_{d}+ \lfloor y(t) \rfloor_{0} -\end{equation} - - -\begin{equation} - \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}}} \cdot \xi(t) -\end{equation} -% das stimmt so, das c kommt unter die Wurzel! - +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 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{Noise_split_intrinsic} -% V_{m_{i+1}} = V_{m_i} + \left(-V_{m_i} + \mu + \alpha V_{d_i} - A_i + \sqrt{\frac{2D c_{\rm{noise}}}{\Delta t}}\mathcal{N}(0,\,1)_i\right) \frac{\Delta t}{\tau_m} +% \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}). - -In the here used model a small portion of the original noise was assigned to the noise component ($c_{\rm{noise}} = 0.1$, \subfigrefb{flowchart}\,\panel[iii]{C}) and a big portion used as the signal component ($c_{\rm{signal}} = 0.9$, \subfigrefb{flowchart}\,\panel[iii]{A}). For the noise split to be valid \citep{Lindner2022} both components must add up to the initial 100\,\% of the total noise and the baseline properties as the firing rate and the CV of the model are maintained. This is easily achieved in a model without a carrier if the condition $c_{\rm{signal}}+c_{\rm{noise}}=1$ is satisfied. The situation here is more complicated. After the original noise was split into a signal component with $c_{\rm{signal}}$, the frequencies above 300\,Hz were discarded and the signal strength was reduced after the dendritic low pass filtering. To compensate for these transformations the initial signal component was multiplied with the factor $\rho$, keeping the baseline CV (only carrier) and the CV during the noise split comparable, and resulting in $s_\xi(t)$. $\rho$ was found by bisecting the plane of possible factors and minimizing the error between the CV during baseline and stimulation. - -%that was found by minimizing the error between the -%Furutsu-Novikov Theorem \citep{Novikov1965, Furutsu1963}\Eqnref{eq:ram_split}, (red in \subfigrefb{flowchart}\,\panel[iii]{A}) bisecting the space of possible $\rho$ scaling factors -%$\rho$ a scaling factor that compensates (see below) for the signal transformations the amplitude modulation stimulus undergoes in the model, i.e. the threshold and the dendritic lowpass. -%In our case the model has a carrier (the fish's self-generated EOD) and we thus want to drive the model with an amplitude modulation stimulus - - -% See section \ref{lifmethods} for model and parameter description. -\begin{table*}[hp!] - \caption{\label{modelparams} Model parameters of LIF models, fitted to 2 electrophysiologically recorded P-units \citep{Ott2020}.} - \begin{center} +\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 $\beta$ & \bfseries $\tau_{m}$/ms & \bfseries $\mu$ & \bfseries $D$/$\mathbf{ms}$ & \bfseries $\tau_{A}$/ms & \bfseries $\Delta_A$ & \bfseries $\tau_{d}$/ms & \bfseries $t_{ref}$/ms \\\hline -2012-07-03-ak& $10.6$& $1.38$& $-1.32$& $0.001$& $96.05$& $0.01$& $1.18$& $0.12$ \\ -2013-01-08-aa& $4.5$& $1.20$& $0.59$& $0.001$& $37.52$& $0.01$& $1.18$& $0.38$ \\ -2018-05-08-ae& $139.6$& $1.49$& $-21.09$& $0.214$& $123.69$& $0.16$& $3.93$& $1.31$ \\ - \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{center} -\end{table*}% 2013-01-08-aa % 2012-07-03-ak - -%\notejb{to methods: ``Note that the signal component \signalnoise{} is added as an amplitude modulation and is thus limited with respect to its spectral content by the Nyquist frequency of the carrier, half the EOD frequency. It thus has a reduced high-frequency content as compared to the intrinsic noise. Adding these discarded high-frequency components to the intrinsic noise does not affect the results here (not shown).''} - +\end{table*} + +\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. -% Either type in your references using -% \begin{thebibliography}{} -% \bibitem{} -% Text -% \end{thebibliography} -% -% or -% -% Compile your BiBTeX database using our plos2015.bst -% style file and paste the contents of your .bbl file -% here. See http://journals.plos.org/plosone/s/latex for -% step-by-step instructions. -% %\bibliographystyle{apalike}%alpha}%}%alpha}%apalike} \bibliography{journalsabbrv,references} % \bibliographystyle{apalike} %or any other style you like %\bibliography{references} %\bibliography{journalsabbrv,references} + + +\setcounter{figure}{0} +\renewcommand{\thefigure}{S\arabic{figure}} + \newpage \section{Supporting information} @@ -896,7 +831,7 @@ CVs in P-units can range up to 1.5 \citep{Grewe2017, Hladnik2023}. We show the s \label{S1:highcvpunit} \begin{figure*}[!ht] \includegraphics[width=\columnwidth]{cells_suscept_high_CV.pdf} - \caption{\label{fig:cells_suscept_high_CV} Response of experimentally measured noisy P-units (cell identifier ``2018-08-24-af") with a relatively high CV of 0.34 to RAM stimuli with two different contrasts. \figitem{A} Interspike intervals (ISI) distribution during baseline. \figitem{B} Baseline power spectrum. \figitem{C} Top: EOD carrier (gray) with RAM (red). Center: Spike trains in response to the 5\,\% RAM contrast. Bottom: Spike trains in response to the 10\,\% RAM contrast. \figitem{D} First-order susceptibility (\Eqnref{linearencoding_methods}). \figitem{E} Absolute value $|\chi_2(f_1, f_2)|$ of the second-order susceptibility, \Eqnref{eq:susceptibility}, for the 5\,\% RAM contrast. Pink lines -- edges of the structure when \fone, \ftwo{} or \fsum{} are equal to \fbase{}. \figitem{F} $|\chi_2(f_1, f_2)|$ for the 10\,\% RAM contrast. \figitem{G} Projected diagonals, calculated as the mean of the anti-diagonals of the matrices in \panel{E--F}. Gray dots: \fbase{}. Dashed lines: Medians of the projected diagonals.} + \caption{\label{fig:cells_suscept_high_CV} Response of a noisy P-units (cell ``2018-08-24-af") with a relatively high baseline CV of 0.34 to RAM stimuli with two different contrasts. \figitem{A} Interspike intervals (ISI) distribution during baseline. \figitem{B} Baseline power spectrum. \figitem{C} Top: EOD carrier (gray) with RAM (red). Center: Spike trains in response to the 5\,\% RAM contrast. Bottom: Spike trains in response to the 10\,\% RAM contrast. \figitem{D} First-order susceptibility (\eqnref{linearencoding_methods}). \figitem{E} Absolute value $|\chi_2(f_1, f_2)|$ of the second-order susceptibility, \eqnref{eq:susceptibility}, for the 5\,\% RAM contrast. Pink lines -- edges of the structure when \fone, \ftwo{} or \fsum{} are equal to \fbase{}. \figitem{F} $|\chi_2(f_1, f_2)|$ for the 10\,\% RAM contrast. \figitem{G} Projected diagonals, calculated as the mean of the anti-diagonals of the matrices in \panel{E--F}. Gray dots: \fbase{}. Dashed lines: Medians of the projected diagonals.} \end{figure*} @@ -907,29 +842,3 @@ CVs in P-units can range up to 1.5 \citep{Grewe2017, Hladnik2023}. We show the s \end{figure*} \end{document} - - - - -%\begin{itemize} -%\item \notejb{ \citep{French1973} Derivation of the Fourier transformed kernels measured with white noise.} -%\item \notejb{ \citep{French1976} Technical issues and tests of Fourier transformed kernels measured with white noise.} -%\item \notejb{ \citep{Victor1977} Cat retinal ganglion cells, gratings with sum of 6 or 8 sinusoids. X - versus Y cells. Peak at f1 == f2 in Y cells. X-cells rather linear. Discussion of mechanism, where a nonlinearity comes in along the pathway} -%\item \notejb{ \citep{Marmarelis1972} Temporal 2nd order kernels, how well do kernels predict responses, catfish retinal ganglion cells} -%\item \notejb{ \citep{Marmarelis1973} Temporal 2nd order kernels, how well do kernels predict responses} -%\item \notejb{ \citep{Victor1988} Cat retinal ganglion cells, the sum of sinusoids, very technical, one measurement similar to \citep{Victor1977}.} -%\item \notejb{\citep{Nikias1993} Third order spectra or bispectra. Very technical overview to higher order spectra} -%\item \notejb{ \citep{Mitsis2007} Spider mechanoreceptor. Linear filters, multivariate nonlinearity, and threshold. The second-order kernel is needed for this. Gaussian noise stimuli.} -%\item \notejb{ \citep{French2001} Time kernels up to 3rd order for predicting spider mechanoreceptor responses (spikes!)} -%\item \notejb{ \citep{French1999} Review on time domain nonlinear systems identification} -%\item \notejb{ \citep{Temchin2005, RecioSpinosa2005} 2nd order Wiener kernel for predicting chinchilla auditory nerve fiber firing rate responses. Strong 2nd order blob at characteristic frequency} -%\item \notejb{ \citep{Schanze1997} lots of bispectra, visual cortex MUA recordings} - -%\item \notejb{ \citep{Theunissen1996} Linear backward stimulus reconstruction in the context of information theory/signal-to-noise ratios} -%\item \notejb{ \citep{Wessel1996} Same as Theunissen1996 but for P-units} -%\item \notejb{ \citep{Neiman2011} cross bispectrum, bicoherence, mutual information, saturating nonlinearities, `` ampullary electroreceptors of paddlefish are perfectly suited to linearly encode weak low-frequency stimuli.''} - -%\item \notejb{ \citep{Chichilnisky2001} Linear Nonlinear Poisson model} -%\item \notejb{ \citep{Gollisch2009} Linear Nonlinear models in retina} -%\item \notejb{ \citep{Clemens2013} Grasshoppper model for female preferences} -%\end{itemize} diff --git a/trialnr.pdf b/trialnr.pdf index a7d105a..e728aa0 100644 Binary files a/trialnr.pdf and b/trialnr.pdf differ