benda/susceptibility1
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

improved methods

Browse Source
This commit is contained in:
Jan Benda 2024-09-24 18:30:05 +02:00
parent 4aca9b2cf2
commit 38ec09d91c
2 changed files with 94 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
references.bib
View File

@ -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.},

164
susceptibility1.tex
View File

@ -416,11 +416,11 @@ Neuronal processing is inherently nonlinear --- spiking thresholds or rectificat
\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 first-oder susceptibility or transfer function 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; 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.
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.
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.
Noise in the 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.
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}.
@ -496,7 +496,7 @@ In the example recordings shown above (\figsrefb{fig:cells_suscept} and \fref{fi
\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}
@ -568,14 +568,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}
@ -615,169 +616,157 @@ 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. 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. 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). 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. 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}
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.
The power spectrum of the stimulus $s(t)$ was calculated as
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) \cdot 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 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}
\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
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{med}(D(f))}
\nli{} = \frac{ \max D(\fbase{}-5\,\rm{Hz} \leq f \leq \fbase{}+5\,\rm{Hz})}{\mathrm{median}(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)$.
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 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}.
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.
\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.
In the model, the input \carrierinput{} was then first thresholded to model the synapse between the primary receptor cells and the afferent.
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)
\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
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{eifnl}
% f(V_m)= \Delta_V \text{e}^{\frac{V_m-1}{\Delta_V}}
% \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}
% \citep{Fourcaud-Trocme2003}, where $\Delta_V$ was varied from 0.001 to 0.1.
%, \figrefb{eif}
\begin{figure*}[t]
\includegraphics[width=\columnwidth]{flowchart.pdf}
@ -785,29 +774,22 @@ Whenever the membrane voltage $V_m(t)$ crossed the spiking threshold $\theta=1$
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}$:
\paragraph{Model fit 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.
\paragraph{Stimuli for the model}
The model neurons were driven with similar stimuli as the real neurons in the electrophysiological experiments. 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
\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}
\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}
according to their contrasts, $c_1$ and $c_2$ at the position of the receiving fish.
\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
Random amplitude modulations (RAMs) were generated by first generating the band-limited white noise stimulus, $s(t)$. For this, random amplitudes and phases were drawn from Gaussian distributions for each frequency component in the range from 0 to 300\,Hz in the Fourier domain. By means of the inverse Fourier transform, the time course of the RAM stimulus was generated. 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.
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.
\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.