From 1e1f66b1ea40843ac4a3f5deae1f21ea2eccc0d2 Mon Sep 17 00:00:00 2001 From: Jan Benda Date: Mon, 30 Jun 2025 10:59:22 +0200 Subject: [PATCH] finished spectral analysis in methods --- nonlinearbaseline.tex | 135 +++++++++++++++++++++--------------------- 1 file changed, 68 insertions(+), 67 deletions(-) diff --git a/nonlinearbaseline.tex b/nonlinearbaseline.tex index 15fb2da..d067462 100644 --- a/nonlinearbaseline.tex +++ b/nonlinearbaseline.tex @@ -399,6 +399,8 @@ \item nonlinear coding \end{keywords} +\notejb{SI ist ein SNR wenn mit r stimuliert wird! Discuss this when introduced.} + % Please keep the abstract below 300 words \section{Abstract} Neuronal processing is inherently nonlinear --- spiking thresholds or rectification at synapses is essential to neuronal computations. Nevertheless, linear response theory has been instrumental in understanding, for example, the impact of noise or neuronal synchrony on signal transmission, or the emergence of oscillatory activity, but is valid only at low stimulus amplitudes or large levels of intrinsic noise. At higher signal-to-noise ratios, however, nonlinear response components become relevant. Theoretical results for leaky integrate-and-fire neurons in the weakly nonlinear regime suggest strong responses at the sum of two input frequencies if these frequencies or their sum match the neuron's baseline firing rate. @@ -471,12 +473,12 @@ P-units fire action potentials probabilistically phase-locked to the self-genera \begin{figure*}[p] \includegraphics[width=\columnwidth]{punitexamplecell} - \caption{\label{fig:punit} Linear and nonlinear stimulus encoding in example P-units. \figitem{A} Interspike interval (ISI) distribution of a cell's baseline activity, i.e. the cell is driven only by the unperturbed own electric field (cell identifier ``2020-10-27-ag''). This cell has a rather high baseline firing rate $r=405$\,Hz and an intermediate CV$_{\text{base}}=0.49$ of its interspike intervals. \figitem{B} Power spectral density of the cell's baseline response with marked peaks at the cell's baseline firing rate $r$ and the fish's EOD frequency $f_{\text{EOD}}$. \figitem{C} Random amplitude modulation (RAM) stimulus (top, red, with cutoff frequency of 300\,Hz) and evoked responses (spike raster, bottom) of the same P-unit for two different stimulus contrasts (right). The stimulus contrast quantifies the standard deviation of the RAM relative to the fish's EOD amplitude. \figitem{D} Gain of the transfer function (first-order susceptibility), \eqnref{linearencoding_methods}, computed from the responses to 10\,\% (light blue) and 20\,\% contrast (dark blue) RAM stimulation of 5\,s duration. \figitem{E} Absolute value of the second-order susceptibility, \eqnref{eq:susceptibility}, for both the low and high stimulus contrast. At the lower stimulus contrast an anti-diagonal where the sum of the two stimulus frequencies equals the neuron's baseline frequency clearly sticks out of the noise floor. \figitem{F} At the higher contrast, the anti-diagonal is much weaker. \figitem{G} Second-order susceptibilities projected onto the diagonal (averages over all anti-diagonals of the matrices shown in \panel{E, F}). The anti-diagonals from \panel{E} and \panel{F} show up as a peak close to the cell's baseline firing rate $r$. The susceptibility index, SI($r$), quantifies the height of this peak relative to the values in the vicinity \notejb{See equation XXX}. \figitem{H} ISI distributions (top) and second-order susceptibilities (bottom) of two more example P-units (``2021-06-18-ae'', ``2017-07-18-ai'') showing an anti-diagonal, but not the full expected triangular structure. \figitem{I} Most P-units, however, have a flat second-order susceptibility and consequently their SI($r$) values are close to one (cell identifiers ``2018-08-24-ak'', ``2018-08-14-ac'').} + \caption{\label{fig:punit} Linear and nonlinear stimulus encoding in example P-units. \figitem{A} Interspike interval (ISI) distribution of a cell's baseline activity, i.e. the cell is driven only by the unperturbed own electric field (cell identifier ``2020-10-27-ag''). This cell has a rather high baseline firing rate $r=405$\,Hz and an intermediate CV$_{\text{base}}=0.49$ of its interspike intervals. \figitem{B} Power spectral density of the cell's baseline response with marked peaks at the cell's baseline firing rate $r$ and the fish's EOD frequency $f_{\text{EOD}}$. \figitem{C} Random amplitude modulation (RAM) stimulus (top, red, with cutoff frequency of 300\,Hz) and evoked responses (spike raster, bottom) of the same P-unit for two different stimulus contrasts (right). The stimulus contrast quantifies the standard deviation of the RAM relative to the fish's EOD amplitude. \figitem{D} Gain of the transfer function (first-order susceptibility), \eqnref{linearencoding_methods}, computed from the responses to 10\,\% (light blue) and 20\,\% contrast (dark blue) RAM stimulation of 5\,s duration. \figitem{E} Absolute value of the second-order susceptibility, \eqnref{chi2}, for both the low and high stimulus contrast. At the lower stimulus contrast an anti-diagonal where the sum of the two stimulus frequencies equals the neuron's baseline frequency clearly sticks out of the noise floor. \figitem{F} At the higher contrast, the anti-diagonal is much weaker. \figitem{G} Second-order susceptibilities projected onto the diagonal (averages over all anti-diagonals of the matrices shown in \panel{E, F}). The anti-diagonals from \panel{E} and \panel{F} show up as a peak close to the cell's baseline firing rate $r$. The susceptibility index, SI($r$), quantifies the height of this peak relative to the values in the vicinity \notejb{See equation XXX}. \figitem{H} ISI distributions (top) and second-order susceptibilities (bottom) of two more example P-units (``2021-06-18-ae'', ``2017-07-18-ai'') showing an anti-diagonal, but not the full expected triangular structure. \figitem{I} Most P-units, however, have a flat second-order susceptibility and consequently their SI($r$) values are close to one (cell identifiers ``2018-08-24-ak'', ``2018-08-14-ac'').} \end{figure*} Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:punit}{C}, top trace, red line), have been commonly used to characterize stimulus-driven responses of sensory neurons using transfer functions (first-order susceptibility), spike-triggered averages, or stimulus-response coherences. Here, we additionally estimate from existing recordings the second-order susceptibility to quantify nonlinear encoding. P-unit spikes align more or less clearly with fluctuations in the RAM stimulus. A higher stimulus intensity, here a higher contrast of the RAM relative to the EOD amplitude (see methods), entrains the P-unit response more clearly (light and dark blue for low and high contrast stimuli, respectively, \subfigrefb{fig:punit}{C}). Linear encoding, quantified by the first-order susceptibility or transfer function, \eqnref{linearencoding_methods}, is similar for the two RAM contrasts in this low-CV P-unit (\subfigrefb{fig:punit}{D}), as expected for a linear system. The first-order susceptibility is low for low frequencies, peaks in the range below 100\,Hz and then falls off again \citep{Benda2005}. -The second-order susceptibility, \eqnref{eq:susceptibility}, quantifies for each combination of two stimulus frequencies \fone{} and \ftwo{} the amplitude and phase of the stimulus-evoked response at the sum \fsum{} (and also the difference, \subfigrefb{fig:model_full}{A}). Large values of the second-order susceptibility indicate stimulus-evoked peaks in the response spectrum at the summed frequency that cannot be explained by linear response theory. Similar to the first-order susceptibility, the second-order susceptibility can be estimated directly from the response evoked by a RAM stimulus that stimulates the neuron with a whole range of frequencies simultaneously (\subfigsref{fig:punit}{E, F}). For LIF and theta neuron models driven in the supra-threshold regime, theory predicts nonlinear interactions between the two stimulus frequencies, when the two frequencies \fone{} and \ftwo{} or their sum \fsum{} exactly match the neuron's baseline firing rate $r$ \citep{Voronenko2017,Franzen2023}. Only then, additional stimulus-evoked peaks appear in the spectrum of the spiking response that would show up in the second-order susceptibility as a horizontal, a vertical, and an anti-diagonal line (\subfigrefb{fig:lifresponse}{B}). +The second-order susceptibility, \eqnref{chi2}, quantifies for each combination of two stimulus frequencies \fone{} and \ftwo{} the amplitude and phase of the stimulus-evoked response at the sum \fsum{} (and also the difference, \subfigrefb{fig:model_full}{A}). Large values of the second-order susceptibility indicate stimulus-evoked peaks in the response spectrum at the summed frequency that cannot be explained by linear response theory. Similar to the first-order susceptibility, the second-order susceptibility can be estimated directly from the response evoked by a RAM stimulus that stimulates the neuron with a whole range of frequencies simultaneously (\subfigsref{fig:punit}{E, F}). For LIF and theta neuron models driven in the supra-threshold regime, theory predicts nonlinear interactions between the two stimulus frequencies, when the two frequencies \fone{} and \ftwo{} or their sum \fsum{} exactly match the neuron's baseline firing rate $r$ \citep{Voronenko2017,Franzen2023}. Only then, additional stimulus-evoked peaks appear in the spectrum of the spiking response that would show up in the second-order susceptibility as a horizontal, a vertical, and an anti-diagonal line (\subfigrefb{fig:lifresponse}{B}). For the example P-unit, we observe a ridge of elevated second-order susceptibility for the low RAM contrast at \fsumb{} (yellowish anti-diagonal, \subfigrefb{fig:punit}{E}). This structure is less prominent for the stronger stimulus (\subfigref{fig:punit}{F}). Further, the overall level of the second-order susceptibility is reduced with increasing stimulus strength. To quantify the structural changes in the susceptibility matrices we projected the susceptibility values onto the diagonal (white dashed line) by averaging over the anti-diagonals (\subfigrefb{fig:punit}{G}). At low RAM contrast this projection indeed has a distinct peak close to the neuron's baseline firing rate (\subfigrefb{fig:punit}{G}, dot on top line). For the higher RAM contrast this peak is much smaller and the overall level of the second-order susceptibility is reduced (\subfigrefb{fig:punit}{G}). The reason behind this reduction is that a RAM with a higher contrast is not only a stimulus with an increased amplitude, but also increases the total noise in the system. Increased noise is known to linearize signal transmission \citep{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017} and thus the second-order susceptibility is expected to decrease. @@ -487,7 +489,7 @@ Overall we observed in 17\,\% of the 159 P-units ridges where the stimulus frequ \begin{figure*}[p] \includegraphics[width=\columnwidth]{ampullaryexamplecell} - \caption{\label{fig:ampullary} Linear and nonlinear stimulus encoding in example ampullary afferents. \figitem{A} Interspike interval (ISI) distribution of the cell's baseline activity (cell identifier ``2012-05-15-ac''). The very low CV of the ISIs indicates almost perfect periodic spiking. \figitem{B} Power spectral density of baseline activity with peaks at the cell's baseline firing rate and its harmonics. Ampullary afferents do not respond to the fish's EOD frequency, $f_{\text{EOD}}$ --- a sharp peak at $f_{\text{EOD}}$ is missing. \figitem{C} Band-limited white noise stimulus (top, red, with a cutoff frequency of 150\,Hz) added to the fish's self-generated electric field (no amplitude modulation!) and spike raster of the evoked responses (bottom) for two stimulus contrasts as indicated (right). \figitem{D} Gain of the transfer function, \eqnref{linearencoding_methods}, of the responses to stimulation with 5\,\% (light green) and 10\,\% contrast (dark green) of 10\,s duration. \figitem{E, F} Absolute value of the second-order susceptibility, \eqnref{eq:susceptibility}, for both stimulus contrasts as indicated. Both show a clear anti-diagonal where the two stimulus frequencies add up to the afferent's baseline firing rate. \figitem{G} Projections of the second-order susceptibilities in \panel{E, F} onto the diagonal. \figitem{H} ISI distributions (top) and second-order susceptibilites (bottom) of three more example afferents with clear anti-diagonals (``2010-11-26-an'', ``2010-11-08-aa'', ``2011-02-18-ab''). \figitem{I} Some ampullary afferents do not show any structure in their second-order susceptibility (``2014-01-16-aj'').} + \caption{\label{fig:ampullary} Linear and nonlinear stimulus encoding in example ampullary afferents. \figitem{A} Interspike interval (ISI) distribution of the cell's baseline activity (cell identifier ``2012-05-15-ac''). The very low CV of the ISIs indicates almost perfect periodic spiking. \figitem{B} Power spectral density of baseline activity with peaks at the cell's baseline firing rate and its harmonics. Ampullary afferents do not respond to the fish's EOD frequency, $f_{\text{EOD}}$ --- a sharp peak at $f_{\text{EOD}}$ is missing. \figitem{C} Band-limited white noise stimulus (top, red, with a cutoff frequency of 150\,Hz) added to the fish's self-generated electric field (no amplitude modulation!) and spike raster of the evoked responses (bottom) for two stimulus contrasts as indicated (right). \figitem{D} Gain of the transfer function, \eqnref{linearencoding_methods}, of the responses to stimulation with 5\,\% (light green) and 10\,\% contrast (dark green) of 10\,s duration. \figitem{E, F} Absolute value of the second-order susceptibility, \eqnref{chi2}, for both stimulus contrasts as indicated. Both show a clear anti-diagonal where the two stimulus frequencies add up to the afferent's baseline firing rate. \figitem{G} Projections of the second-order susceptibilities in \panel{E, F} onto the diagonal. \figitem{H} ISI distributions (top) and second-order susceptibilites (bottom) of three more example afferents with clear anti-diagonals (``2010-11-26-an'', ``2010-11-08-aa'', ``2011-02-18-ab''). \figitem{I} Some ampullary afferents do not show any structure in their second-order susceptibility (``2014-01-16-aj'').} \end{figure*} Electric fish possess an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogenous electric stimuli. The population of ampullary afferents is much less heterogeneous, and known for the much lower CVs of their baseline ISIs ($0.06 < \text{CV}_{\text{base}} < 0.22$, \citealp{Grewe2017}). Ampullary cells do not phase-lock to the high-frequency 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 $r$ and its harmonics. Since the cells do not respond to the self-generated EOD, there is no sharp peak at \feod{} (\subfigrefb{fig:ampullary}{B}). When driven by a band-limited white noise stimulus (note: for ampullary afferents this is not an AM stimulus, \subfigref{fig:ampullary}{C}), ampullary afferents exhibit very pronounced ridges in the second-order susceptibility, where $f_1 + f_2$ is equal to $r$ or its harmonics (yellow anti-diagonals in \subfigrefb{fig:ampullary}{E--H}), implying strong nonlinear response components at these frequency combinations (\subfigrefb{fig:ampullary}{G}, top). With higher stimulus contrasts these bands get weaker (\subfigrefb{fig:ampullary}{F}), the projection onto the diagonal loses its distinct peak at $r$, and its overall level is reduced (\subfigrefb{fig:ampullary}{G}, bottom). Some ampullary afferents (27\,\% of 30 afferents), however, do not show any such structure in their second-order susceptibility (\subfigrefb{fig:ampullary}{I}). @@ -501,18 +503,18 @@ In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampull \caption{\label{fig:noisesplit} Estimation of second-order susceptibilities. \figitem{A} \suscept{} (right) estimated from $N=198$ 256\,ms long FFT segments of an electrophysiological recording of another P-unit (cell ``2017-07-18-ai'', $r=78$\,Hz, CV$_{\text{base}}=0.22$) driven with a RAM stimulus with contrast 5\,\% (left). \figitem[i]{B} \textit{Standard condition} of model simulations with intrinsic noise (bottom) and a RAM stimulus (top). \figitem[ii]{B} \suscept{} estimated from simulations of the cell's LIF model counterpart (cell ``2017-07-18-ai'', table~\ref{modelparams}) based on the same RAM contrast and number of $N=100$ FFT segments. As in the electrophysiological recording only a weak anti-diagonal is visible. \figitem[iii]{B} Same as \panel[ii]{B} but using $10^6$ FFT segments. Now, the expected triangular structure is revealed. \figitem[iv]{B} Convergence of the \suscept{} estimate as a function of FFT segments. \figitem{C} At a lower stimulus contrast of 1\,\% the estimate did not converge yet even for $10^6$ FFT segments. The triangular structure is not revealed yet. \figitem[i]{D} Same as in \panel[i]{B} but in the \textit{noise split} condition: there is no external RAM signal (red) driving the model. Instead, a large part (90\,\%) of the total intrinsic noise is treated as a signal and is presented as an equivalent amplitude modulation ($s_{\xi}(t)$, orange, 10.6\,\% contrast), while the intrinsic noise is reduced to 10\,\% of its original strength (bottom, see methods for details). \figitem[i]{D} 100 FFT segments are still not sufficient for estimating \suscept{}. \figitem[iii]{D} Simulating one million segments reveals the full expected trangular structure of the second-order susceptibility. \figitem[iv]{D} In the noise-split condition, the \suscept{} estimate converges already at about $10^{4}$ FFT segments.} \end{figure*} -One simple reason could be the lack of data, i.e. the estimation of the second-order susceptibility is not good enough. Electrophysiological recordings are limited in time, and therefore only a limited number of trials, here repeated presentations of the same frozen RAM stimulus, are available. In our data set we have 1 to 199 trials (median: 10) of RAM stimuli with a duration ranging from 2 to 5\,s (median: 8\,s), resulting in a total duration of 30 to 400\,s. Using a resolution of 0.5\,ms and FFT segments of 512 samples this yields 105 to 1520 available FFT segments for a specific RAM stimulus. As a consequence, the cross-spectra, \eqnref{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 the P-unit estimated from the same low number of FFT (fast fourier transform) segments as in the experiment ($N=100$, compare faint anti-diagonal in the bottom left corner of the second-order susceptibility in \panel[ii]{A} and \panel[ii]{B} in \figrefb{fig:noisesplit}). +One simple reason could be the lack of data, i.e. the estimation of the second-order susceptibility is not good enough. Electrophysiological recordings are limited in time, and therefore only a limited number of trials, here repeated presentations of the same frozen RAM stimulus, are available. In our data set we have 1 to 199 trials (median: 10) of RAM stimuli with a duration ranging from 2 to 5\,s (median: 8\,s), resulting in a total duration of 30 to 400\,s. Using a resolution of 0.5\,ms and FFT segments of 512 samples this yields 105 to 1520 available FFT segments for a specific RAM stimulus. As a consequence, the cross-spectra, \eqnref{crossxss}, are insufficiently averaged and the full structure of the second-order susceptibility might be hidden in finite-data noise. This experimental limitation can be overcome by using a computational model for the P-unit, a stochastic leaky integrate-and-fire model with adaptation current, dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of the P-unit estimated from the same low number of FFT (fast fourier transform) segments as in the experiment ($N=100$, compare faint anti-diagonal in the bottom left corner of the second-order susceptibility in \panel[ii]{A} and \panel[ii]{B} in \figrefb{fig:noisesplit}). In model simulations we can increase the number of FFT segments beyond what would be experimentally possible, here to one million (\figrefb{fig:noisesplit}\,\panel[iii]{B}). Then the estimate of the second-order susceptibility indeed improves. It gets less noisy, the diagonal at $f_ + f_2 = r$ is emphasized, and the vertical and horizontal ridges at $f_1 = r$ and $f_2 = r$ are revealed. Increasing the number of FFT segments also reduces the order of magnitude of the susceptibility estimate until close to one million the estimate levels out at a low values (\subfigrefb{fig:noisesplit}\,\panel[iv]{B}). At a lower stimulus contrast of 1\,\% (\subfigrefb{fig:noisesplit}{C}), however, one million FFT segments are still not sufficient for the estimate to converge (\figrefb{fig:noisesplit}\,\panel[iv]{C}). Still only a faint anti-diagonal is visible (\figrefb{fig:noisesplit}\,\panel[iii]{C}). -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. As we just have seen, 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$ \notejb{$c$ is already used for contrast} (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \figrefb{fig:noisesplit}\,\panel[i]{D}). We tuned the amplitude of the RAM stimulus $s_{\xi}(t)$ such that the output firing rate and variability (CV of interspike intervals) 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 and thus the estimate converges already at about ten thousand FFT segments (\figrefb{fig:noisesplit}\,\panel[iv]{D}); (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 FFT segments (\figrefb{fig:noisesplit}\,\panel[iii]{D}), but not for a number of segments comparable to the experiment (\figrefb{fig:noisesplit}\,\panel[ii]{D}). In addition to the strong response for \fsumb{}, we now also observe pronounced nonlinear responses at \foneb{} and \ftwob{} (vertical and horizontal lines, \figrefb{fig:noisesplit}\,\panel[iii]{D}). +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. As we just have seen, 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{crossxss} and \eqref{chi2}) \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$ \notejb{$c$ is already used for contrast} (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \figrefb{fig:noisesplit}\,\panel[i]{D}). We tuned the amplitude of the RAM stimulus $s_{\xi}(t)$ such that the output firing rate and variability (CV of interspike intervals) 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 and thus the estimate converges already at about ten thousand FFT segments (\figrefb{fig:noisesplit}\,\panel[iv]{D}); (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 FFT segments (\figrefb{fig:noisesplit}\,\panel[iii]{D}), but not for a number of segments comparable to the experiment (\figrefb{fig:noisesplit}\,\panel[ii]{D}). In addition to the strong response for \fsumb{}, we now also observe pronounced nonlinear responses at \foneb{} and \ftwob{} (vertical and horizontal lines, \figrefb{fig:noisesplit}\,\panel[iii]{D}). \begin{figure}[p] \includegraphics[width=\columnwidth]{modelsusceptcontrasts} - \caption{\label{fig:modelsusceptcontrasts}Dependence of second order susceptibility on stimulus contrast. \figitem{A} Second-order susceptibilities estimated for increasing stimulus contrasts of $c=0, 1, 3$ and $10$\,\% as indicated ($N=10^7$ FFT segments for $c=1$\,\%, $N=10^6$ segments for all other contrasts). $c=0$\,\% refers to the noise-split configuration (limit to vanishing external RAM signal, see \subfigrefb{fig:noisesplit}{D}). Shown are simulations of the P-unit model cell ``2017-07-18-ai'' (table~\ref{modelparams}) with a baseline firing rate of $82$\,Hz and CV$_{\text{base}}=0.23$. The cell shows a clear triangular pattern in its second-order susceptibility even up to a contrast of $10$\,\%. Note, that for $c=1$\,\% (\panel[ii]{D}), the estimate did not converge yet. \figitem{B} Cell ``2012-12-13-ao'' (baseline firing rate of $146$\,Hz, CV$=0.23$) also has strong interactions at its baseline firing rate that survive up to a stimulus contrast of $3$\,\%. \figitem{C} Model cell ``2012-12-20-ac'' (baseline firing rate of $212$\,Hz, CV$=0.26$) shows a weak triangular structure in the second-order susceptibility that vanishes for stimulus contrasts larger than $1$\,\%. \figitem{D} Cell ``2013-01-08-ab'' (baseline firing rate of $218$\,Hz, CV$=0.55$) does not show any triangular pattern in its second-order susceptibility. Nevertheless, interactions between low stimulus frequencies become substantial at higher contrasts. \figitem{E} The presence of an elevated second-order susceptibility where the stimulus frequency add up to the neuron's baseline frequency, can be identified by the susceptibility index (SI($r$), \eqnref{eq:nli_equation2}) greater than one (horizontal black line). The SI($r$) (density to the right) is plotted as a function of the model neuron's baseline CV for all $39$ model cells (table~\ref{modelparams}). Model cells have been visually categorized based on the strong (11 cells) or weak (5 cells) presence of a triangular pattern in their second-order susceptibility estimated in the noise-split configuration (legend). The cells from \panel{A--D} are marked by black circles. Pearson's correlation coefficients $r$, the corresponding significance level $p$ and regression line (dashed gray line) are indicated. The higher the stimulus contrast, the less cells show weakly nonlinear interactions as expressed by the triangular structure in the second-order susceptibility.} + \caption{\label{fig:modelsusceptcontrasts}Dependence of second order susceptibility on stimulus contrast. \figitem{A} Second-order susceptibilities estimated for increasing stimulus contrasts of $c=0, 1, 3$ and $10$\,\% as indicated ($N=10^7$ FFT segments for $c=1$\,\%, $N=10^6$ segments for all other contrasts). $c=0$\,\% refers to the noise-split configuration (limit to vanishing external RAM signal, see \subfigrefb{fig:noisesplit}{D}). Shown are simulations of the P-unit model cell ``2017-07-18-ai'' (table~\ref{modelparams}) with a baseline firing rate of $82$\,Hz and CV$_{\text{base}}=0.23$. The cell shows a clear triangular pattern in its second-order susceptibility even up to a contrast of $10$\,\%. Note, that for $c=1$\,\% (\panel[ii]{D}), the estimate did not converge yet. \figitem{B} Cell ``2012-12-13-ao'' (baseline firing rate of $146$\,Hz, CV$=0.23$) also has strong interactions at its baseline firing rate that survive up to a stimulus contrast of $3$\,\%. \figitem{C} Model cell ``2012-12-20-ac'' (baseline firing rate of $212$\,Hz, CV$=0.26$) shows a weak triangular structure in the second-order susceptibility that vanishes for stimulus contrasts larger than $1$\,\%. \figitem{D} Cell ``2013-01-08-ab'' (baseline firing rate of $218$\,Hz, CV$=0.55$) does not show any triangular pattern in its second-order susceptibility. Nevertheless, interactions between low stimulus frequencies become substantial at higher contrasts. \figitem{E} The presence of an elevated second-order susceptibility where the stimulus frequency add up to the neuron's baseline frequency, can be identified by the susceptibility index (SI($r$), \eqnref{siindex2}) greater than one (horizontal black line). The SI($r$) (density to the right) is plotted as a function of the model neuron's baseline CV for all $39$ model cells (table~\ref{modelparams}). Model cells have been visually categorized based on the strong (11 cells) or weak (5 cells) presence of a triangular pattern in their second-order susceptibility estimated in the noise-split configuration (legend). The cells from \panel{A--D} are marked by black circles. Pearson's correlation coefficients $r$, the corresponding significance level $p$ and regression line (dashed gray line) are indicated. The higher the stimulus contrast, the less cells show weakly nonlinear interactions as expressed by the triangular structure in the second-order susceptibility.} \end{figure} \subsection{Weakly nonlinear interactions in many model cells} @@ -520,7 +522,7 @@ In the previous section we have shown one example cell for which we find in the By just looking at the second-order susceptibilities estimated using the noise-split method (first column of \figrefb{fig:modelsusceptcontrasts}) we can readily identify strong triangular patterns in 11 of the 39 model cells (28\,\%, see \figrefb{fig:modelsusceptcontrasts}\,\panel[i]{A}\&\panel[i]{B} for two examples). In another 5 cells (13\,\%) the triangle is much weaker and sits on top of a smooth bump of elevated second-order susceptibility (\figrefb{fig:modelsusceptcontrasts}\,\panel[i]{C} shows an example). The remaining 23 model cells (59\,\%) show no triangle (see \figrefb{fig:modelsusceptcontrasts}\,\panel[i]{D} for an example). -This categorization is supported by the susceptibility index, SI($r$), \eqnref{eq:nli_equation2}, which quantifies the height of the ridge where the stimulus frequencies add up to the neuron's baseline firing rate relative to the background. Values above one indicate an elevated ridge. The absence of such a ridge results in values close to one. Indeed, the cells showing only a weak triangle (orange) arise out of values around one and the cells showing strong triangles (red) have consistently SI($r$) values exceeding 1.8 (\figrefb{fig:modelsusceptcontrasts}\,\panel[i]{E}). +This categorization is supported by the susceptibility index, SI($r$), \eqnref{siindex2}, which quantifies the height of the ridge where the stimulus frequencies add up to the neuron's baseline firing rate relative to the background. Values above one indicate an elevated ridge. The absence of such a ridge results in values close to one. Indeed, the cells showing only a weak triangle (orange) arise out of values around one and the cells showing strong triangles (red) have consistently SI($r$) values exceeding 1.8 (\figrefb{fig:modelsusceptcontrasts}\,\panel[i]{E}). The SI($r$) correlates with the cell's CV of its baseline interspike intervals ($r=-0.60$, $p<0.001$). The lower the cell's CV$_{\text{base}}$, the higher the SI($r$) value and thus the stronger the triangular structure of its second-order susceptibility. The model cells with the most distinct triangular pattern in their second-order susceptibility are the ones with the lowest CVs, hinting at low intrinsic noise levels. @@ -546,7 +548,7 @@ Overall, observing SI($r$) values greater than about 1.8, even for a number of F \includegraphics[width=\columnwidth]{dataoverview} \caption{\label{fig:dataoverview} Nonlinear responses in P-units and ampullary afferents. The second-order susceptibility is condensed - into the susceptibility index, SI($r$) \eqnref{eq:nli_equation}, + into the susceptibility index, SI($r$) \eqnref{siindex}, that quantifies the relative amplitude of the ridge where the two stimulus frequencies add up to the cell's baseline firing rate $r$ (see \subfigrefb{fig:punit}{G}). In both the models and the @@ -588,7 +590,7 @@ Overall, observing SI($r$) values greater than about 1.8, even for a number of F \end{figure*} \subsection{Low CVs and weak responses predict weakly nonlinear responses} -Now we are prepared to evaluate our pool of 39 P-unit model cells, 172 P-units, and 30 ampullary afferents recorded in 80 specimen of \textit{Apteronotus leptorhynchus}. For comparison across cells we condensed the structure of the second-order susceptibilities into SI($r$) values, \eqnref{eq:nli_equation}. In order to make the data comparable, both model and experimental SI($r$) estimates, \eqnref{eq:nli_equation}, are based on 100 FFT segments. +Now we are prepared to evaluate our pool of 39 P-unit model cells, 172 P-units, and 30 ampullary afferents recorded in 80 specimen of \textit{Apteronotus leptorhynchus}. For comparison across cells we condensed the structure of the second-order susceptibilities into SI($r$) values, \eqnref{siindex}. In order to make the data comparable, both model and experimental SI($r$) estimates, \eqnref{siindex}, are based on 100 FFT segments. In the P-unit models, each model cell contributed with three RAM stimulus presentations with contrasts of 1, 3, and 10\,\%, resulting in $n=117$ data points. 19 (16\,\%) had SI($r$) values larger than 1.8, indicating the expected ridges at the baseline firing rate in their second-order susceptiility. The lower the cell's baseline CV, i.e. the less intrinsic noise, the higher the SI($r$) (\figrefb{fig:dataoverview}\,\panel[i]{A}). @@ -604,13 +606,13 @@ The population of ampullary cells is generally more homogeneous, with lower base \begin{figure*}[t] \includegraphics[width=\columnwidth]{model_full.pdf} - \caption{\label{fig:model_full} Second-order susceptibility qualitatively predicts responses to sine-wave stimuli. \figitem[]{A} Absolute value of the second-order susceptibility, \eqnref{eq:susceptibility}, for both positive and negative frequencies. \susceptf{} was estimated from $N=10^6$ FFT segments of model simulations in the noise-split condition (cell ``2013-01-08-aa'', table~\ref{modelparams}). Dashed white lines mark zero frequencies. Nonlinear responses at \fsum{} are quantified in the upper right and lower left quadrants. Nonlinear responses at \fdiff{} are quantified in the upper left and lower right quadrants. The baseline firing rate of this cell was at $\fbase=120$\,Hz. The position of the orange/red letters corresponds to the beat frequencies used for the stimulation with pure sine waves in the subsequent panels and indicates the sum/difference of those beat frequencies. \figitem{B--E} Black line -- power spectral density of model simulations in response to stimulation with two pure sine waves, \fone{} and \ftwo, in addition to the receiving fish's own EOD (three-fish scenario). The contrast of beat beats is 2\,\%. Colored circles highlight the height of selected peaks in the power spectrum. Grey line -- power spectral density of model in the baseline condition. \figitem{B} The sum of the two beat frequencies match \fbase{}. \figitem{C} The difference of \fone{} and \ftwo{} match \fbase{}. \figitem{D} Only the first beat frequency matches \fbase{}. \figitem{E} None of the two beat frequencies matches \fbase{}.} + \caption{\label{fig:model_full} Second-order susceptibility qualitatively predicts responses to sine-wave stimuli. \figitem[]{A} Absolute value of the second-order susceptibility, \eqnref{chi2}, for both positive and negative frequencies. \susceptf{} was estimated from $N=10^6$ FFT segments of model simulations in the noise-split condition (cell ``2013-01-08-aa'', table~\ref{modelparams}). Dashed white lines mark zero frequencies. Nonlinear responses at \fsum{} are quantified in the upper right and lower left quadrants. Nonlinear responses at \fdiff{} are quantified in the upper left and lower right quadrants. The baseline firing rate of this cell was at $\fbase=120$\,Hz. The position of the orange/red letters corresponds to the beat frequencies used for the stimulation with pure sine waves in the subsequent panels and indicates the sum/difference of those beat frequencies. \figitem{B--E} Black line -- power spectral density of model simulations in response to stimulation with two pure sine waves, \fone{} and \ftwo, in addition to the receiving fish's own EOD (three-fish scenario). The contrast of beat beats is 2\,\%. Colored circles highlight the height of selected peaks in the power spectrum. Grey line -- power spectral density of model in the baseline condition. \figitem{B} The sum of the two beat frequencies match \fbase{}. \figitem{C} The difference of \fone{} and \ftwo{} match \fbase{}. \figitem{D} Only the first beat frequency matches \fbase{}. \figitem{E} None of the two beat frequencies matches \fbase{}.} \end{figure*} \subsection{Second-order susceptibility explains nonlinear responses to pure sinewave stimuli} Using the RAM stimulation we found pronounced nonlinear responses, in particular 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 $f_1$ and $f_2$ 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, $f_2$ was similar to the baseline firing rate, corresponding to the horizontal ridge of the second-order susceptibility (\subfigrefb{fig:lifresponse}{B}). The difference frequency is not covered by the so-far shown part of the second-order susceptibility, in which only the response at the sum of the two stimulus frequencies is addressed. -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 $\chi_2$ matrix is symmetric with respect to the origin. In the upper-right and lower-left quadrants of \susceptf{}, stimulus-evoked responses at the sum $f_1 + f_2$ are shown, whereas in the lower-right and upper-left quadrants responses at the difference $f_2 - f_1$ are shown (\figref{fig:model_full}). The vertical and horizontal lines at $f_1 = r$ and $f_2 = r$ are very pronounced in the upper-right quadrant of \subfigrefb{fig:model_full}{A} and extend into the upper-left quadrant, where they fade out towards more negative $f_1$ frequencies. The peak in the response power-spectrum at $f_2 - f_1$ 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}). +However, the second-order susceptibility \eqnref{chi2} is a spectral measure that is based on Fourier transforms and thus is also defined for negative stimulus frequencies. The full $\chi_2$ matrix is symmetric with respect to the origin. In the upper-right and lower-left quadrants of \susceptf{}, stimulus-evoked responses at the sum $f_1 + f_2$ are shown, whereas in the lower-right and upper-left quadrants responses at the difference $f_2 - f_1$ are shown (\figref{fig:model_full}). The vertical and horizontal lines at $f_1 = r$ and $f_2 = r$ are very pronounced in the upper-right quadrant of \subfigrefb{fig:model_full}{A} and extend into the upper-left quadrant, where they fade out towards more negative $f_1$ frequencies. The peak in the response power-spectrum at $f_2 - f_1$ 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 in the noise-split configuration (\subfigrefb{fig:model_full}{A})? We test this by stimulating the same model with two weak beats (\subfigrefb{fig:model_full}{B--E}). Qualitatively, the second-order susceptibility predicts the presence and absence of peaks in the response spectrum at the sums and differences of the two beat frequencies. However, a quantitative prediction fails. This is because pure sine waves influence a nonlinear system in different ways than a white-noise stimulus. This will be addressed in detail in a future manuscript. @@ -659,7 +661,7 @@ The weakly nonlinear interactions in low-CV P-units could facilitate the detecta \subsection{Conclusions} -We have demonstrated pronounced nonlinear responses in primary electrosensory afferents at weak stimulus amplitudes and sufficiently low intrinsic noise levels. The observed nonlinearities match the expectations from previous theoretical studies \citep{Voronenko2017,Franzen2023}. The resulting nonlinear components introduce spectral components not present in the original stimulus, but may provide an edge in the context of signal detection problems at stimulus amplitudes close to threshold \citep{Schlungbaum2023}. Electrosenory afferents share an evolutionary history with hair cells \citep{Baker2019} and share many response properties with mammalian auditory nerve fibers \citep{Barayeu2023, Joris2004}. Thus, we expect weakly nonlinear responses for near-threshold stimulation in auditory nerve fibers as well. +We have demonstrated pronounced nonlinear responses in primary electrosensory afferents at weak stimulus amplitudes and sufficiently low intrinsic noise levels. The observed nonlinearities match the expectations from previous theoretical studies \citep{Voronenko2017,Franzen2023}. The resulting nonlinear components introduce spectral components not present in the original stimulus, but may provide an edge in the context of signal detection problems at stimulus amplitudes close to threshold \citep{Schlungbaum2023}. Electrosensory afferents share an evolutionary history with hair cells \citep{Baker2019} and share many response properties with mammalian auditory nerve fibers \citep{Barayeu2023, Joris2004}. Thus, we expect weakly nonlinear responses for near-threshold stimulation in auditory nerve fibers as well. \section{Acknowledgements} @@ -696,48 +698,51 @@ The fish were stimulated with band-limited white noise stimuli with a cut-off fr \subsection{Data analysis} Data analysis was done in Python 3 using the packages matplotlib \citep{Hunter2007}, numpy \citep{Walt2011}, scipy \citep{scipy2020}, pandas \citep{Mckinney2010}, nixio \citep{Stoewer2014}, and thunderlab (\url{https://github.com/bendalab/thunderlab}). \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 (median 32\,s). The coefficient of variation (CV) of the interspike intervals (ISI) is their standard deviation relative to their mean: $\rm{CV}_{\rm base} = \sqrt{\langle (ISI- \langle ISI \rangle) ^2 \rangle} / \langle ISI \rangle$. If the baseline was recorded several times in a recording, the measures from the longest recording were taken. +The baseline firing rate $r$ was calculated as the number of spikes divided by the duration of the baseline recording (median 32\,s). The coefficient of variation (CV) of the interspike intervals (ISI) is their standard deviation relative to their mean: $\rm{CV}_{\rm base} = \sqrt{\langle (ISI- \langle ISI \rangle) ^2 \rangle} / \langle ISI \rangle$. If the baseline was recorded several times in a recording, the measures from the longest recording were taken. \paragraph{White noise analysis} \label{response_modulation} -In the stimulus-driven case, the neuronal activity of the recorded cell is modulated around the average firing rate that is similar to \fbase{} and in that way encodes the time-course of the stimulus. Spiking activity -\begin{equation} - \label{eq:spikes} - x_k(t) = \sum_i\delta(t-t_{k,i}) -\end{equation} -is recorded for each stimulus presentation $k$, as a train of times $t_{k,i}$ where action potentials occured. -%If only a single trial was recorded or is used for the analysis, we drop the trial index $k$. +When stimulated with band-limited white noise stimuli, neuronal activity is modulated around the average firing rate that is similar to the baseline firing rate and in that way encodes the time-course of the stimulus. For an estimate of the time-dependent firing rate $r(t)$ we convolved each spike train with normalized Gaussian kernels with standard deviation of 2\,ms and averaged the resulting single-trail firing rates over trials. The response modulation quantifies the variation of $r(t)$ computed as the standard deviation in time $\sigma_{s} = \sqrt{\langle (r(t)-\langle r(t) \rangle_t )^2\rangle_t}$, where $\langle \cdot \rangle_t$ denotes averaging over time. -The single-trial firing rate -\begin{equation} - r_k(t) = x_k(t) * K(t) = \int_{-\infty}^{+\infty} x_k(t') K(t-t') \, \mathrm{d}t' -\end{equation} -was estimated by convolving the spike train with a kernel. We used a Gaussian kernel -\begin{equation} -K(t) = {\scriptstyle \frac{1}{\sigma\sqrt{2\pi}}} e^{-\frac{t^2}{2\sigma^2}} -\end{equation} -with standard deviation $\sigma$ set to 2\,ms. Averaging over $n$ repeated stimulus presentations results in the trial averaged firing rate -\begin{equation} - \label{eq:rate} - r(t) = \left\langle r_k(t) \right\rangle _k = \frac{1}{n} \sum_{k=1}^n r_k(t) -\end{equation} +% Spiking activity +% \begin{equation} +% \label{eq:spikes} +% x_k(t) = \sum_i\delta(t-t_{k,i}) +% \end{equation} +% is recorded for each stimulus presentation $k$, as a train of times $t_{k,i}$ where action potentials occured. +% %If only a single trial was recorded or is used for the analysis, we drop the trial index $k$. -The 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. +% The single-trial firing rate +% \begin{equation} +% r_k(t) = x_k(t) * K(t) = \int_{-\infty}^{+\infty} x_k(t') K(t-t') \, \mathrm{d}t' +% \end{equation} +% was estimated by convolving the spike train with a kernel. We used a Gaussian kernel +% \begin{equation} +% K(t) = {\scriptstyle \frac{1}{\sigma\sqrt{2\pi}}} e^{-\frac{t^2}{2\sigma^2}} +% \end{equation} +% with standard deviation $\sigma$ set to 2\,ms. Averaging over $n$ repeated stimulus presentations results in the trial averaged firing rate +% \begin{equation} +% \label{eq:rate} +% r(t) = \left\langle r_k(t) \right\rangle _k = \frac{1}{n} \sum_{k=1}^n r_k(t) +% \end{equation} + +% The average firing rate during stimulation, $r_s = \langle r(t) \rangle_t$, is given by the temporal average $\langle \cdot \rangle_t$ over the duration of the stimulus of the trial-averaged firing rate. To quantify how strongly a neuron was driven by the stimulus, we computed the response modulation as the standard deviation $\sigma_{s} = \sqrt{\langle (r(t)-r_s)^2\rangle_t}$ of the trial-averaged firing rate. \paragraph{Spectral analysis}\label{susceptibility_methods} -The neuron is driven by the stimulus and thus its spiking response depends on the time course of the stimulus. To characterize the relation between stimulus $s(t)$ and response $x(t)$, we calculated the first- and second-order susceptibilities in the frequency domain. +To characterize the relation between the spiking response evoked by white-noise stimuli, we estimated the first- and second-order susceptibilities in the frequency domain. For this we converted spike times into binary vectors $x_k$ with $\Delta t = 0.5$\,ms wide bins that are set to 2\,kHz where a spike occurred and zero otherwise. Fast Fourier transforms (FFT) $S(\omega)$ and $X(\omega)$ of the stimulus $s_k$ (also down-sampled to a sampling rate of 2\,kHz) and $x_k$, respectively, were computed according to $X(\omega) = \sum_{k=0}^{n} \, x_k e^{- i \omega k}$ for $n=512$ long segments of 256\,ms duration with no overlap, resulting in a spectral resolution of about 4\,Hz. For simplicity we use angular frequencies $\omega=2\pi f$ instead of frequencies $f$. + +In the experimental data the duration of the noise stimuli varied and they were presented once or repeatedly (frozen noise). For the analysis we discarded the responses within the initial 200\,ms of stimulation in each trial. To make the recordings comparable we always used the first 100 segments from as many trials as needed for the following analysis. -Fast fourier transforms (FFT) $\tilde s_T(\omega)$ and $\tilde x_T(\omega)$ of $s(t)$ and $x(t)$, respectively, were computed according to $\tilde x_T(\omega) = \int_{0}^{T} \, x(t) e^{- i \omega t}\,dt$ for $T=0.5$\,s long segments with no overlap, resulting in a spectral resolution of 2\,Hz for the experimental data. For simplicity we use angular frequencies $\omega=2\pi f$ instead of frequencies $f$. In the experimental data, most stimuli had a duration of 10\,s and were chopped into 20 segments. Spectral measures were computed for single trials of neural responses, series of spike times, \eqnref{eq:spikes}. For spectral parameters used in simulations, see below. -\notejb{Update FFT parameter. Describe binary spiketrain normalized by dt and plain FFT. Power spectra scaled by dt/nfft, bispectra by dt**2/nfft.} +In the simulations we generated for each trial a new realization of the noise stimulus. We discarded the first 500\,ms of the response and used the following 10 FFT segments for the analysis. This was repeated until $10^6$ or $10^7$ FFT segments were collected. The power spectrum of the stimulus $s(t)$ was estimated as \begin{equation} - \label{powereq} - S_{ss}(\omega) = \frac{\langle \tilde s_T(\omega) \tilde s_T^*(\omega)\rangle}{T} + \label{powerss} + S_{ss}(\omega) = \frac{\Delta t}{n} \langle S(\omega) S^*(\omega) \rangle \end{equation} -with $\tilde s^*$ being the complex conjugate of $\tilde s$ and $\langle \cdot \rangle$ denoting the average over the segments. The power spectrum of the spike trains $S_{xx}(\omega)$ was estimated accordingly. The cross-spectrum $S_{xs}(\omega)$ between stimulus and evoked spike trains was estimated according to +with $S^*$ being the complex conjugate of $S$ and $\langle \cdot \rangle$ denoting the average over the FFT segments. The power spectrum of the spike trains $S_{xx}(\omega)$ was estimated accordingly. Likewise, the cross-spectrum $S_{xs}(\omega)$ between stimulus and evoked spike trains was estimated according to \begin{equation} - \label{cross} - S_{xs}(\omega) = \frac{\langle \tilde x_T(\omega) \tilde s_T^*(\omega)\rangle}{T} + \label{crossxs} + S_{xs}(\omega) = \frac{\Delta t}{n} \langle X(\omega) S^*(\omega) \rangle \end{equation} The first-order susceptibility (transfer function) @@ -745,43 +750,39 @@ The first-order susceptibility (transfer function) \label{linearencoding_methods} \chi_{1}(\omega) = \frac{S_{xs}(\omega) }{S_{ss}(\omega) } \end{equation} -was then computed from $S_{xs}(\omega)$ and $ S_{ss}(\omega)$. +was then computed from $S_{xs}(\omega)$ and $ S_{ss}(\omega)$. We report $\chi_{1}(\omega)$ in Hz/\%, i.e. firing rate per percent stimulus contrast. Multiplying $\chi_{1}(\omega)$ with the contrast of a sinusoidal stimulus in percent results in the amplitude of the evoked firing rate modulation in Hertz. The second-order cross-spectrum \begin{equation} - \label{eq:crosshigh} - S_{xss} (\omega_{1},\omega_{2}) = \frac{\langle \tilde x_T(\omega_{1} + \omega_{2}) \tilde s_T^*(\omega_{1}) \tilde s_T^*(\omega_{2}) \rangle}{T} + \label{crossxss} + S_{xss}(\omega_{1},\omega_{2}) = \frac{\Delta t^2}{n} \langle X(\omega_{1} + \omega_{2}) S^*(\omega_{1}) S^*(\omega_{2}) \end{equation} -describes nonlinear interactions that generate responses at the sum and difference (at negative $\omega_2$) of two stimulus frequencies $\omega_1$ and $\omega_2$. +quantifies nonlinear interactions that generate responses at the sum and difference (for negative $\omega_2$) evoked by two stimulus frequencies $\omega_1$ and $\omega_2$. The second-order susceptibility \begin{equation} - \label{eq:susceptibility} + \label{chi2} \chi_{2}(\omega_{1}, \omega_{2}) = \frac{S_{xss} (\omega_{1},\omega_{2})}{2 S_{ss} (\omega_{1}) S_{ss} (\omega_{2})} \end{equation} -normalizes the second-order cross-spectrum by the spectral power at the two stimulus frequencies. -% Applying the Fourier transform this can be rewritten resulting in: -% \begin{equation} -% \label{susceptibility} -% \chi_{2}(\omega_{1}, \omega_{2}) = \frac{TN \sum_{n=1}^N \int_{0}^{T} dt\,r_{n}(t) e^{-i(\omega_{1}+\omega_{2})t} \int_{0}^{T}dt'\,s_{n}(t')e^{i \omega_{1}t'} \int_{0}^{T} dt''\,s_{n}(t'')e^{i \omega_{2}t''}}{2 \sum_{n=1}^N \int_{0}^{T} dt\, s_{n}(t)e^{-i \omega_{1}t} \int_{0}^{T} dt'\,s_{n}(t')e^{i \omega_{1}t'} \sum_{n=1}^N \int_{0}^{T} dt\,s_{n}(t)e^{-i \omega_{2}t} \int_{0}^{T} dt'\,s_{n}(t')e^{i \omega_{2}t'}} -% \end{equation} +normalizes the second-order cross-spectrum by the spectral power at the two stimulus frequencies. We report $\chi_{2}(\omega_{1}, \omega_{2})$ in Hz/\%$^2$, i.e. firing rate per percent stimulus contrast squared. Throughout the manuscript we only show the absolute values of the complex-valued second-order susceptibility matrix and ignore the corresponding phases. -\paragraph{Nonlinearity index}\label{projected_method} -We expected to see elevated values in the second-order susceptibility at $\omega_1 + \omega_2 = \fbase$ \citep{Voronenko2017,Franzen2023}. To characterize this in a single number we computed the peakedness of the nonlinearity via +\paragraph{Susceptibility index}\label{projected_method} +We expected to see a sharp ridge in the second-order susceptibility at $\omega_1 + \omega_2 = r$ \citep{Voronenko2017,Franzen2023}. To characterize this in a single number we computed a susceptibility index. First, we projected the absolute values of the second-order susceptibility matrix onto the diagonal by averaging over anti-diagonal elements. In this projection $D(f)$ we took the position of the maximum \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))} + \label{dmax} + f_{\rm peak} = {\rm argmax} D(\fbase{}-50\,\rm{Hz} \leq f \leq r + 50\,\rm{Hz}) \end{equation} -from the second-order susceptibilities estimated from experimental data. $D(f)$ is the projection of the absolute values of the second-order susceptibility matrix onto the diagonal, calculated by taking the mean of the anti-diagonal elements. The peak of $D(f)$ at the neuron's baseline firing rate $\fbase$ was found by finding the maximum of $D(f)$ in the range $\fbase \pm 5$\,Hz. This peak was then normalized by the median of $D(f)$ (\subfigrefb{fig:punit}{G}). Since in most experimentally measured cells the second-order susceptibilities was more or less flat, normalizing by the median is working well. If the same RAM was recorded several times in a cell, each trial resulted in a separate second-order susceptibility matrix. For the population statistics in \figref{fig:dataoverview} the mean of the resulting \nli{} values is used. - -\notejb{computed from FFTs based on 100 segments of 512 samples at a temporal resolution of 0.5\,ms} - -For the model simulations, the second-order susceptibilities where not flat but often showed a broader bump on which the peak at the neuron's baseline firing rate occured (see, for example, \figrefb{fig:modelsusceptcontrasts}, panels \panel[i]{B}, \panel[iii]{C}, and \panel[iv]{D}). Instead of normalizing by the median, we normalized the model simulations by the mean values of the projection $D(f)$ in small ranges close to the left and right of the neuron's baseline firing rate: +within $\pm 50$\,Hz of the neuron's baseline firing rate $r$ as the position of the expected peak. For an estimate of the noise-floor surrounding this peak we averaged over 10\,Hz wide windows 10\,Hz to the left and right of the peak: +\begin{equation} + \label{dref} + D_{\rm ref} = \frac{1}{2}\left(\langle D(f_{\rm peak}-20\,\rm{Hz} \leq f \leq f_{\rm peak}-10\,\rm{Hz})\rangle_f + \langle D(f_{\rm peak}+10\,\rm{Hz} \leq f \leq f_{\rm peak}+20\,\rm{Hz}))\rangle_f\right) +\end{equation} +The size of the peak relative to this reference is then the susceptibility index \begin{equation} - \label{eq:nli_equation2} - \nli{} = \frac{\max D(\fbase{}-5\,\rm{Hz} \leq f \leq \fbase{}+5\,\rm{Hz})}{\frac{1}{2}\left[\mathrm{mean}(D(\fbase{}-20\,\rm{Hz} \leq f \leq \fbase{}-10\,\rm{Hz})) + \mathrm{mean}(D(\fbase{}+10\,\rm{Hz} \leq f \leq \fbase{}+20\,\rm{Hz}))\right]} + \label{siindex} + SI = \frac{D(f_{\rm peak})}{D_{\rm ref}} \end{equation} -Using \eqnref{eq:nli_equation} instead of \eqnref{eq:nli_equation2} for the estimates of the second-order susceptibility based on $N=10$ segments does not change the results in \subfigref{fig:modelsusceptlown}{C \& D}. That is, experimental data based on low numbers of segments can be safely analyzed using the more robust \eqnref{eq:nli_equation}. +Values larger than one indicate a sharp ridge in the susceptibility matrix close to where the stimulus frequencies add up to the baseline firing rate. \begin{figure*}[t] @@ -867,7 +868,7 @@ The eight free parameters of the P-unit model $\beta$, $\tau_m$, $\mu$, $D$, $\t \end{tabular} \end{table*} -In \figrefb{fig:noisesplit} and \figrefb{fig:model_full}, for each of the stimulus and response segments needed for the spectral analysis, \eqnsref{powereq}--\eqref{eq:susceptibility}, a simulation was run. The first second was discarded and the analysis was based on the last second of the data. The resulting spectra thus have a spectral resolution of 1\,Hz. To speed up model simulations for \figrefb{fig:modelsusceptcontrasts} and \figrefb{fig:trialnr} we simulated up to $10^6$ trials of 3\,s duration. In each trial the first 0.5\,s were skipped, yielding 10 FFT segments of $T=0.25$\,s (4\,Hz resolution). +In \figrefb{fig:noisesplit} and \figrefb{fig:model_full}, for each of the stimulus and response segments needed for the spectral analysis, \eqnsref{powerss}--\eqref{chi2}, a simulation was run. The first second was discarded and the analysis was based on the last second of the data. The resulting spectra thus have a spectral resolution of 1\,Hz. To speed up model simulations for \figrefb{fig:modelsusceptcontrasts} and \figrefb{fig:trialnr} we simulated up to $10^6$ trials of 3\,s duration. In each trial the first 0.5\,s were skipped, yielding 10 FFT segments of $T=0.25$\,s (4\,Hz resolution). \subsection{Noise split} @@ -881,7 +882,7 @@ Based on the Furutsu-Novikov theorem \citep{Furutsu1963,Novikov1965,Lindner2022, \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. +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{crossxss}, were computed between the RAM stimulus $s_{\xi}(t)$ and the spike train $x(t)$ it evoked. In this way, the effective signal-to-noise ratio can be increased while maintaining the total noise in the system. %\bibliographystyle{apalike}%alpha}%}%alpha}%apalike} \bibliography{journalsabbrv,references}