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

updated results to new figures

Browse Source
This commit is contained in:
Jan Benda 2025-05-26 23:52:21 +02:00
parent d56c0c7c7e
commit fcf9cfa431
6 changed files with 44 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
ampullaryexamplecell.py
View File

@ -134,7 +134,8 @@ if __name__ == '__main__':
if k == 0: if k == 0:
axc.set_ylabel('') axc.set_ylabel('')
fig.common_yticks(axs[1, :]) fig.common_yticks(axs[1, :])
fig.tag([axs[0, :]], xoffs=-3, yoffs=1) fig.tag([axs[0, :3]], xoffs=-3, yoffs=1)
fig.tag(axs[0, 3:])
fig.savefig() fig.savefig()
print() print()

5
dataoverview.py
View File

@ -300,6 +300,11 @@ if __name__ == '__main__':
si = '_nmax' si = '_nmax'
si_thresh = 1.8 si_thresh = 1.8
cells = np.unique(list(punit_model['cell']) + list(punit_data['cell']) + list(ampul_data['cell']))
specimen = np.unique(['-'.join(c.split('-')[:3]) for c in cells])
print(f'Total of {len(cells)} cells recorded in {len(specimen)} specimens')
print()
cvmodel = punit_model['cvbase'] cvmodel = punit_model['cvbase']
cvdata = punit_data['cvbase'] cvdata = punit_data['cvbase']
u, p = mannwhitneyu(cvmodel, cvdata) u, p = mannwhitneyu(cvmodel, cvdata)

10
lifsuscept.py
View File

@ -117,13 +117,7 @@ def plot_gain(ax, s, r0, freqs, chi1):
def plot_chi2(ax, s, r0, freqs, chi2): def plot_chi2(ax, s, r0, freqs, chi2):
chi2 = np.abs(chi2) chi2 = np.abs(chi2)
vmax = np.quantile(chi2, 0.996) vmax = np.quantile(chi2, 0.996)
ten = 10**np.floor(np.log10(vmax)) vmax = 300
for fac, delta in zip([1, 2, 3, 4, 6, 8, 10],
[0.5, 1, 1, 2, 3, 4, 5]):
if fac*ten >= vmax:
vmax = fac*ten
ten *= delta
break
pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax, pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax,
rasterized=True) rasterized=True)
ax.set_aspect('equal') ax.set_aspect('equal')
@ -139,7 +133,7 @@ def plot_chi2(ax, s, r0, freqs, chi2):
cb.outline.set_color('none') cb.outline.set_color('none')
cb.outline.set_linewidth(0) cb.outline.set_linewidth(0)
cax.set_ylabel('$|\\chi_2(f_1, f_2)|$') cax.set_ylabel('$|\\chi_2(f_1, f_2)|$')
cax.set_yticks_delta(ten) cax.set_yticks_delta(100)
if __name__ == '__main__': if __name__ == '__main__':

6
modelsusceptcontrasts.py
View File

@ -127,8 +127,10 @@ def plot_summary_contrasts(axs, s, cells):
print(f' mean SI = {np.mean(nli):.2f}, stdev = {np.std(nli):.2f}') print(f' mean SI = {np.mean(nli):.2f}, stdev = {np.std(nli):.2f}')
n = np.sum(nli > nli_thresh) n = np.sum(nli > nli_thresh)
print(f' {n} cells ({100*n/len(nli):.1f}%) have SI > {nli_thresh:.1f}:') print(f' {n} cells ({100*n/len(nli):.1f}%) have SI > {nli_thresh:.1f}:')
for name, cv in cdata[nli > nli_thresh, ['cell', 'cvbase']].row_data(): for name, cv, respmod in cdata[nli > nli_thresh,
print(f' {name:<22} CV={cv:4.2f}') ['cell', 'cvbase', 'respmod2']].row_data():
print(f' {name:<22} CV={cv:4.2f} '
f'response modulation={respmod:4.0f}Hz')
print(f' triangle cells have SI >= {np.min(nli[cdata["triangle"] > 0.5]):.2f}') print(f' triangle cells have SI >= {np.min(nli[cdata["triangle"] > 0.5]):.2f}')
print() print()
print('overall lowest baseline CV:', print('overall lowest baseline CV:',

71
nonlinearbaseline.tex
View File

@ -441,9 +441,9 @@ Here we search for such weakly nonlinear responses in electroreceptors of the tw
Theoretical work on leaky integrate-and-fire and conductance-based models suggests a distinct structure of the second-order response function for neurons with low levels of intrinsic noise driven in the supra-threshold regime with low stimulus amplitudes (\figrefb{fig:lifresponse}, \citealp{Voronenko2017}). Here, we explored a large set of recordings of P-units and ampullary cells of the active and passive electrosensory systems of the brown ghost knifefish \textit{Apteronotus leptorhynchus} to search for such weakly nonlinear responses in real neurons. Additional simulations of LIF-based models of P-unit spiking put the experimental findings into context. We start with demonstrating the basic concepts using example P-units and models. \notejg{relativ viel Wiederholung des ob gesagten, kuerzen?} Theoretical work on leaky integrate-and-fire and conductance-based models suggests a distinct structure of the second-order response function for neurons with low levels of intrinsic noise driven in the supra-threshold regime with low stimulus amplitudes (\figrefb{fig:lifresponse}, \citealp{Voronenko2017}). Here, we explored a large set of recordings of P-units and ampullary cells of the active and passive electrosensory systems of the brown ghost knifefish \textit{Apteronotus leptorhynchus} to search for such weakly nonlinear responses in real neurons. Additional simulations of LIF-based models of P-unit spiking put the experimental findings into context. We start with demonstrating the basic concepts using example P-units and models. \notejg{relativ viel Wiederholung des ob gesagten, kuerzen?}
\subsection{Nonlinear responses in P-units stimulated with two frequencies} \subsection{Nonlinear responses in P-units stimulated with two frequencies}
Without external stimulation, a P-unit is driven by the fish's own EOD alone (with a specific EOD frequency \feod{}) and spontaneously fires action potentials at the baseline rate \fbase{}. Accordingly, the power spectrum of the baseline activity has a peak at \fbase{} (\subfigrefb{fig:motivation}{A}). Superposition of the receiver's EOD with an EOD of another fish with frequency $f_1$ results in a beat, a periodic amplitude modulation of the receiver's EOD. The frequency of the beat is given by the difference frequency $\Delta f_1 = f_1 - \feod$ between the two fish. P-units encode this beat in their firing rate \citep{Bastian1981a, Barayeu2023} and consequently the power spectrum of the response has a peak at the beat frequency (\subfigrefb{fig:motivation}{B}). A second peak at the first harmonic of the beat frequency is indicative of a nonlinear process that here is easily identified by the clipping of the P-unit's firing rate at zero. Pairing the fish with another fish at a higher beat frequency $\Delta f_2 = f_2 - \feod > \Delta f_1$ results in a weaker response with a single peak in the response power spectrum, suggesting a linear response (\subfigrefb{fig:motivation}{C}). The weaker response to this beat can be explained by the beat tuning of the cell \citep{Walz2014}. Note, $\Delta f_2$ has been deliberately chosen to match the recorded P-unit's baseline firing rate. Without external stimulation, a P-unit is driven by the fish's own EOD alone (with a specific EOD frequency \feod{}) and spontaneously fires action potentials at the baseline rate \fbase{}. Accordingly, the power spectrum of the baseline activity has a peak at \fbase{} (\subfigrefb{fig:motivation}{A}). Superposition of the receiver's EOD with an EOD of another fish with frequency $f_1$ results in a beat, a periodic amplitude modulation of the receiver's EOD. The frequency of the beat is given by the difference frequency $\Delta f_1 = f_1 - \feod$ between the two fish. P-units encode this beat in their firing rate \citep{Bastian1981a} and consequently the power spectrum of the response has a peak at the beat frequency (\subfigrefb{fig:motivation}{B}). A second peak at the first harmonic of the beat frequency is indicative of a nonlinear process that here is associated with the clipping of the P-unit's firing rate at zero \citep{Barayeu2023}. Pairing the fish with another fish at a higher beat frequency $\Delta f_2 = f_2 - \feod > \Delta f_1$ results in a weaker response with a single peak in the response power spectrum, suggesting a linear response (\subfigrefb{fig:motivation}{C}). The weaker response to this beat can be explained by the beat tuning of the cell \citep{Walz2014}. Note, $\Delta f_2$ has been deliberately chosen to match the recorded P-unit's baseline firing rate.
When stimulating with both foreign signals simultaneously, additional peaks appear in the response power spectrum at the sum $\Delta f_1 + \Delta f_2$ and the difference frequency $\Delta f_2 - \Delta f_1$ (\subfigrefb{fig:motivation}{D}). Thus, the cellular response is not equal to the sum of the responses to the two beats presented separately. These peaks at the sum and the difference of the two stimulus frequencies are a hallmark of nonlinear interactions that, by definition, are absent in linear systems\notejg{should not occur?}. When stimulating with both foreign signals simultaneously, additional peaks appear in the response power spectrum at the sum $\Delta f_1 + \Delta f_2$ and the difference frequency $\Delta f_2 - \Delta f_1$ (\subfigrefb{fig:motivation}{D}). Thus, the cellular response is not equal to the sum of the responses to the two beats presented separately. These additional peaks at the sum and the difference of the two stimulus frequencies are a hallmark of nonlinear interactions that, by definition, are absent in linear systems\notejg{should not occur?}.
\subsection{Linear and weakly nonlinear regimes} \subsection{Linear and weakly nonlinear regimes}
@ -455,61 +455,59 @@ When stimulating with both foreign signals simultaneously, additional peaks appe
The stimuli used in \figref{fig:motivation} had the same not-small amplitude. Whether this stimulus condition falls into the weakly nonlinear regime as in \citet{Voronenko2017} is not clear. In order to illustrate how the responses to two beat frequencies develop over a range of amplitudes we use a stochastic leaky-integrate-and-fire (LIF) based P-unit model fitted to a specific electrophysiologically measured cell \citep{Barayeu2023}. The stimuli used in \figref{fig:motivation} had the same not-small amplitude. Whether this stimulus condition falls into the weakly nonlinear regime as in \citet{Voronenko2017} is not clear. In order to illustrate how the responses to two beat frequencies develop over a range of amplitudes we use a stochastic leaky-integrate-and-fire (LIF) based P-unit model fitted to a specific electrophysiologically measured cell \citep{Barayeu2023}.
At very low stimulus contrasts (in the example cell less than approximately 1.5\,\% relative to the receiver's EOD amplitude) the spectrum has small peaks only at the beat frequencies (\subfigref{fig:regimes}{A}, green and purple). The amplitudes of these peaks initially increase linearly with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines), an indication of the linear response at lowest stimulus amplitudes. At very low stimulus contrasts (in the example cell less than approximately 1.2\,\% relative to the receiver's EOD amplitude) the spectrum has small peaks only at the beat frequencies (\subfigref{fig:regimes}{A,B}, green and purple). The amplitudes of these peaks initially increase linearly with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines), an indication of the linear response at lowest stimulus amplitudes.
This linear regime is followed by the weakly nonlinear regime (in the example cell between approximately 1.5\,\% and 4\,\% stimulus contrasts). In addition to the peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:regimes}{B}, orange and red). The amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines). Note, that we have chosen $\Delta f_2$ to match the baseline firing rate $f_{base}$ of the neuron. This linear regime is followed by the weakly nonlinear regime (in the example cell between approximately 1.2\,\% and 3.5\,\% stimulus contrasts). In addition to the peaks at the stimulus frequencies, peaks at the sum and the difference of the stimulus frequencies appear in the response spectrum (\subfigref{fig:regimes}{C}, orange and red). The amplitudes of these two peaks initially increase quadratically with stimulus amplitude (\subfigref{fig:regimes}{E}, thin lines). Note, that we have chosen $\Delta f_2$ to match the baseline firing rate $f_{base}$ of the neuron.
At higher stimulus amplitudes, the linear response and the weakly-nonlinear response begin to deviate from their linear and quadratic dependency on amplitude (\subfigrefb{fig:regimes}{C \& E}) and additional peaks appear in the response spectrum (\subfigrefb{fig:regimes}{D}). At high stimulus contrasts, additional nonlinearities in the system, in particular clipping of the firing rate, shape the responses. At higher stimulus amplitudes, the linear response and the weakly-nonlinear response begin to deviate from their linear and quadratic dependency on amplitude (\subfigrefb{fig:regimes}{E}) and additional peaks appear in the response spectrum (\subfigrefb{fig:regimes}{D}). At high stimulus contrasts, additional nonlinearities in the system, in particular clipping of the firing rate, shape the responses.
For this example, we chose very specific stimulus (beat) frequencies. %One of these matching the P-unit's baseline firing rate. For this example, we chose very specific stimulus (beat) frequencies. %One of these matching the P-unit's baseline firing rate.
In the following, however, we are interested in how the nonlinear responses depend on different combinations of stimulus frequencies in the weakly nonlinear regime. For the sake of simplicity we will drop the $\Delta$ notation even though P-unit stimuli are beats. In the following, however, we are interested in how the nonlinear responses depend on different combinations of stimulus frequencies in the weakly nonlinear regime. For the sake of simplicity we will drop the $\Delta$ notation even though P-unit stimuli are beats.
\subsection{Nonlinear signal transmission in low-CV P-units} \subsection{Nonlinear signal transmission in example P-units}
Weakly nonlinear responses are expected in cells with sufficiently low intrinsic noise levels, i.e. low baseline CVs \citep{Voronenko2017}. P-units fire action potentials probabilistically phase-locked to the self-generated EOD \citep{Bastian1981a}. Skipping of EOD cycles leads to the characteristic multimodal ISI distribution with maxima at integer multiples of the EOD period (\subfigrefb{fig:punit}{A}). In this example, the baseline ISI distribution has a CV$_{\text{base}}$ of 0.2, which is at the lower end of the P-unit population \citep{Hladnik2023}. Spectral analysis of the baseline activity shows two major peaks: the first is located at the baseline firing rate \fbase, the second is located at the discharge frequency \feod{} of the electric organ and is flanked by two smaller peaks at $\feod \pm \fbase{}$ (\subfigref{fig:punit}{B}). P-units fire action potentials probabilistically phase-locked to the self-generated EOD \citep{Bastian1981a}. Skipping of EOD cycles leads to the characteristic multimodal ISI distribution with maxima at integer multiples of the EOD period (\subfigrefb{fig:punit}{A}). In this example, the baseline ISI distribution has a CV$_{\text{base}}$ of 0.49, which is at the center of the P-unit population \citep{Hladnik2023}. Spectral analysis of the baseline activity shows two major peaks: the first is located at the baseline firing rate $r$, the second is located at the discharge frequency \feod{} of the electric organ (\subfigref{fig:punit}{B}).
\begin{figure*}[p] \begin{figure*}[p]
\includegraphics[width=\columnwidth]{punitexamplecell} \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 four more example P-units (``2021-06-18-ae'', ``2017-07-18-ai'', ``2018-08-24-ak'', ``2018-08-14-ac'') covering the range of baseline firing rates and CV$_{\text{base}}$s as indicated. The first two cells show an anti-diagonal, but not the full expected triangular structure. The second-order susceptibilities of the other two cells are flat and consequently the SI($r$) values are close to one.} \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'').}
\end{figure*} \end{figure*}
Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:punit}{C}, top trace, red line), are commonly used to characterize stimulus-driven responses of sensory neurons using transfer functions (first-order susceptibility), spike-triggered averages, or stimulus-response coherences. Here, we additionally estimate the second-order susceptibility to quantify nonlinear encoding. P-unit spikes align more or less clearly with fluctuations in the RAM stimulus. A higher stimulus intensity, here a higher contrast of the RAM relative to the EOD amplitude (see methods), entrains the P-unit response more clearly (light and dark purple for low and high contrast stimuli, \subfigrefb{fig:punit}{C}). Linear encoding, quantified by the transfer function, \eqnref{linearencoding_methods}, is similar for the two RAM contrasts in this low-CV P-unit (\subfigrefb{fig:punit}{D}), as expected for a linear system. The first-order susceptibility is low for low frequencies, peaks in the range below 100\,Hz and then falls off again \citep{Benda2005}. Noise stimuli, here random amplitude modulations (RAM) of the EOD (\subfigref{fig:punit}{C}, top trace, red line), are commonly used to characterize stimulus-driven responses of sensory neurons using transfer functions (first-order susceptibility), spike-triggered averages, or stimulus-response coherences. Here, we additionally estimate the second-order susceptibility to quantify nonlinear encoding. P-unit spikes align more or less clearly with fluctuations in the RAM stimulus. A higher stimulus intensity, here a higher contrast of the RAM relative to the EOD amplitude (see methods), entrains the P-unit response more clearly (light and dark 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 \fbase{} \citep{Voronenko2017,Franzen2023}. Only then, additional stimulus-evoked peaks appear in the spectrum of the spiking response that would show up in the second-order susceptibility as a horizontal, a vertical, and an anti-diagonal line (\subfigrefb{fig:lifresponse}{B}, pink triangle in \subfigsref{fig:punit}{E, F}). 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}).
For the low-CV P-unit, we observe a band of slightly elevated second-order susceptibility for the low RAM contrast at \fsumb{} (yellowish anti-diagonal between pink edges, \subfigrefb{fig:punit}{E}). This structure vanishes for the stronger stimulus (\subfigref{fig:punit}{F}). Further, the overall level of the second-order susceptibility is reduced with increasing stimulus strength. To quantify the structural changes in the susceptibility matrices we projected the susceptibility values onto the diagonal by averaging over the anti-diagonals (\subfigrefb{fig:punit}{G}). At low RAM contrast this projected second-order susceptibility indeed has a small peak at \fbase{} (\subfigrefb{fig:punit}{G}, dot on top line). For the higher RAM contrast, however, this peak vanishes and the overall level of the second-order susceptibility is reduced (\subfigrefb{fig:punit}{G}). The reason behind this reduction is that a RAM with a higher contrast is not only a stimulus with an increased amplitude, but also increases the total noise in the system. Increased noise is known to linearize signal transmission \citep{Longtin1993, Chialvo1997, Roddey2000, Voronenko2017} and thus the second-order susceptibility is expected to decrease. 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 projected second-order susceptibility 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.
In contrast, a high-CV P-unit (CV$_{\text{base}}=0.4$) does not exhibit pronounced nonlinearities even at low stimulus contrasts (\figrefb{fig:punithighcv}). Irrespective of the CV, neither of the two example P-units shows the complete expected structure of nonlinear interactions. In other P-units we also observe ridges where the stimulus frequencies add up to the unit's baseline firing rate (\subfigrefb{fig:punit}{H}), but we never observed the expected triangular structure. In most P-units, however, we did not observe any structure in the second-order susceptibility (\subfigrefb{fig:punit}{I}).
\subsection{Ampullary afferents exhibit strong nonlinear interactions} \subsection{Ampullary afferents exhibit strong nonlinear interactions}
\begin{figure*}[p] \begin{figure*}[p]
\includegraphics[width=\columnwidth]{ampullaryexamplecell} \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 four more example afferents (``2010-11-26-an'', ``2010-11-08-aa'', ``2011-02-18-ab'', and ``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{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'').}
\end{figure*} \end{figure*}
Electric fish possess an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogenous electric stimuli. The population of ampullary afferents is much less heterogeneous, and known for the much lower CVs of their baseline ISIs (CV$_{\text{base}}=0.06$ to $0.22$, \citealp{Grewe2017}). Ampullary cells do not phase-lock to the EOD and the ISIs are unimodally distributed (\subfigrefb{fig:ampullary}{A}). As a consequence of the high regularity of their baseline spiking activity, the corresponding power spectrum shows distinct peaks at the baseline firing rate \fbase{} and its harmonics. Since the cells do not respond to the self-generated EOD, there is no peak at \feod{} (\subfigrefb{fig:ampullary}{B}). When driven by a low-contrast noise stimulus (note: this is no longer an AM stimulus, \subfigref{fig:ampullary}{C}), ampullary cells exhibit very pronounced bands in the second-order susceptibility, where \fsum{} is equal to \fbase{} or its harmonic (yellow diagonals in \subfigrefb{fig:ampullary}{E}), implying strong nonlinear response components at these frequency combinations (\subfigrefb{fig:ampullary}{G}, top). With higher stimulus contrasts these bands disappear (\subfigrefb{fig:ampullary}{F}), the projection onto the diagonal loses its distinct peak at \fsum{} and its overall level is reduced (\subfigrefb{fig:ampullary}{G}, bottom). Electric fish possess an additional electrosensory system, the passive or ampullary electrosensory system, that responds to low-frequency exogenous electric stimuli. The population of ampullary afferents is much less heterogeneous, and known for the much lower CVs of their baseline ISIs (CV$_{\text{base}}=0.06$ to $0.22$, \citealp{Grewe2017}). Ampullary cells do not phase-lock to the 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: this is no longer an AM stimulus, \subfigref{fig:ampullary}{C}), ampullary cells exhibit very pronounced bands in the second-order susceptibility, where \fsum{} is equal to \fbase{} or its harmonic (yellow 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 \fsum{} and its overall level is reduced (\subfigrefb{fig:ampullary}{G}, bottom). Some ampullary afferents, however, do not show any such structure in their second-order susceptibility (\subfigrefb{fig:ampullary}{I}).
\subsection{Model-based estimation of the second-order susceptibility} \subsection{Model-based estimation of the second-order susceptibility}
In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagonal of the second-order susceptibility, where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, which is in line with theoretical expectations \citep{Voronenko2017,Franzen2023}. However, a pronounced nonlinear response at frequencies \foneb{} or \ftwob{}, although predicted by theory, cannot be observed. In the following, we investigate how these discrepancies can be understood. In the example recordings shown above (\figsrefb{fig:punit} and \fref{fig:ampullary}), we only observe nonlinear responses on the anti-diagonal of the second-order susceptibility, where the sum of the two stimulus frequencies matches the neuron's baseline firing rate, which is in line with theoretical expectations \citep{Voronenko2017,Franzen2023}. However, a pronounced nonlinear response where one of the stimulus frequencies matches the baseline firing rate (horizontal and vertical ridges), although predicted by theory, cannot be observed. In the following, we investigate how these discrepancies can be understood.
\begin{figure*}[p] \begin{figure*}[p]
\includegraphics[width=\columnwidth]{noisesplit} \includegraphics[width=\columnwidth]{noisesplit}
\caption{\label{fig:noisesplit} Estimation of second-order susceptibilities. \figitem{A} \suscept{} (right) estimated from $N=198$ 256\,ms long FFT segments of an electrophysiological recording of another P-unit (cell ``2017-07-18-ai'', $r=78$\,Hz, CV$_{\text{base}}=0.22$) driven with a RAM stimulus with contrast 5\,\% (left). \figitem[i]{B} \textit{Standard condition} of model simulations with intrinsic noise (bottom) and a RAM stimulus (top). \figitem[ii]{B} \suscept{} estimated from simulations of the cell's LIF model counterpart (cell ``2017-07-18-ai'', table~\ref{modelparams}) based on 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.} \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*} \end{figure*}
%\notejb{Since the model overestimated the sensitivity of the real P-unit, we adjusted the RAM contrast to 0.9\,\%, such that the resulting spike trains had the same CV as the electrophysiological recorded P-unit during the 2.5\,\% contrast stimulation (see table~\ref{modelparams} for model parameters).} \notejb{chi2 scale is higher than in real cell} One simple reason could be the lack of data, i.e. the estimation of the second-order susceptibility is not good enough. Electrophysiological recordings are limited in time, and therefore only a limited number of trials, here repetitions of the same frozen RAM stimulus, are available. As a consequence, the cross-spectra, \eqnref{eq:crosshigh}, are insufficiently averaged and the full structure of the second-order susceptibility might be hidden in finite-data noise. This experimental limitation can be overcome by using a computational model for the P-unit, a stochastic leaky integrate-and-fire model with adaptation current, dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of another example cell 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, i.e. repetitions of the same RAM stimulus, are available. As a consequence, the cross-spectra, \eqnref{eq:crosshigh}, are insufficiently averaged and the full structure of the second-order susceptibility might be hidden in finite-data noise. This experimental limitation can be overcome by using a computational model for the P-unit, a stochastic leaky integrate-and-fire model with adaptation current, dendritic preprocessing, and parameters fitted to the experimentally recorded P-unit (\figrefb{flowchart}) \citep{Barayeu2023}. The model faithfully reproduces the second-order susceptibility of a low-CV cell estimated from the same low number of FFT (fast fourier transform) segments as in the experiment ($\n{}=11$, compare \panel{A} and \panel[ii]{B} in \figrefb{fig:noisesplit}). In model simulations we can increase the number of FFT segments beyond what would be experimentally possible, here to one million (\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 reduces the order of magnitude of the susceptibility estimate until close to one million the estimate levels out at a low levels (\subfigrefb{fig:noisesplit}\,\panel[iv]{B}).
In model simulations we can increase the number of FFT segments beyond what would be experimentally possible, here to one million (\subfigrefb{fig:noisesplit}\,\panel[iii]{B}). In this example, the estimate of the second-order susceptibility indeed improves. It gets less noisy and the diagonal at \fsum{} is emphasized. However, the expected vertical and horizontal lines at \foneb{} and \ftwob{} are still mainly missing. At a lower stimulus contrast of 1\,\% (\subfigrefb{fig:noisesplit}{C}), however, one million FFT segements 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. Again, this cannot be done experimentally, because the problem of insufficient averaging becomes even more severe for weak AMs (low contrast). In the model, however, we know the time course of the intrinsic noise and can use this knowledge to determine the susceptibilities by input-output correlations via the Furutsu-Novikov theorem \citep{Furutsu1963, Novikov1965}. This theorem, in its simplest form, states that the cross-spectrum $S_{x\eta}(\omega)$ of a Gaussian noise $\eta(t)$ driving a nonlinear system and the system's output $x(t)$ is proportional to the linear susceptibility according to $S_{x\eta}(\omega)=\chi(\omega)S_{\eta\eta}(\omega)$. Here $\chi(\omega)$ characterizes the linear response to an infinitely weak signal $s(t)$ in the presence of the background noise $\eta(t)$. Likewise, the nonlinear susceptibility can be determined in an analogous fashion from higher-order input-output cross-spectra (see methods, equations \eqref{eq:crosshigh} and \eqref{eq:susceptibility}) \citep{Egerland2020}. In line with an alternative derivation of the Furutsu-Novikov theorem \citep{Lindner2022}, we can split the total noise and consider a fraction of it as a stimulus. This allows us to calculate the susceptibility from the cross-spectrum between the output and this stimulus fraction of the noise. Adapting this approach to our P-unit model (see methods), we replace the intrinsic noise by an approximately equivalent RAM stimulus $s_{\xi}(t)$ and a weak remaining intrinsic noise $\sqrt{2D \, c_{\rm{noise}}}\;\xi(t)$ with $c_\text{noise} = 0.1$ (see methods, equations \eqref{eq:ram_split}, \eqref{eq:Noise_split_intrinsic}, \eqref{eq:Noise_split_intrinsic_dendrite}, \subfigrefb{fig:noisesplit}\,\panel[i]{C}). We tune the amplitude of the RAM stimulus $s_{\xi}(t)$ such that the output firing rate and variability (CV) are the same as in the baseline activity (i.e. full intrinsic noise $\sqrt{2D}\;\xi(t)$ in the voltage equation but no RAM) and compute the cross-spectra between the RAM part of the noise $s_{\xi}(t)$ and the output spike train. This procedure has two consequences: (i) by means of the cross-spectrum between the output and \signalnoise, which is a large fraction of the noise, the signal-to-noise ratio of the measured susceptibilities is drastically improved; (ii) the total noise in the system has been reduced (by what was before the external RAM stimulus $s(t)$), which makes the system more nonlinear. For both reasons we now see the expected nonlinear features in the second-order susceptibility for a sufficient number of segments (\subfigrefb{fig:noisesplit}\,\panel[iii]{C}), but not for a number of segments comparable to the experiment (\subfigrefb{fig:noisesplit}\,\panel[ii]{C}). In addition to the strong response for \fsumb{}, we now also observe pronounced nonlinear responses at \foneb{} and \ftwob{} (vertical and horizontal lines, \subfigrefb{fig:noisesplit}\,\panel[iii]{C}). 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$ (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) 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}).
Note, that the increased number of segments goes along with a substantial reduction of second-order susceptibility values (\subfigrefb{fig:noisesplit}\,\panel[iii]{C}). Only for more than about $10^5$ segments does the estimate of the second-order susceptibility converge for most of the model cells.
\begin{figure}[p] \begin{figure}[p]
@ -520,17 +518,17 @@ Note, that the increased number of segments goes along with a substantial reduct
\subsection{Weakly nonlinear interactions in many model cells} \subsection{Weakly nonlinear interactions in many model cells}
In the previous section we have shown one example cell for which we find in the corresponding model strong ridges in the second-order susceptibility where one \notejg{same here...} or the sum of two stimulus frequencies match the neuron's baseline firing rate (\figrefb{fig:noisesplit}\,\panel[iii]{C}). Using our 39 P-unit models, we now can explore how many P-unit model neurons show such a triangular structure. In the previous section we have shown one example cell for which we find in the corresponding model strong ridges in the second-order susceptibility where one \notejg{same here...} or the sum of two stimulus frequencies match the neuron's baseline firing rate (\figrefb{fig:noisesplit}\,\panel[iii]{C}). Using our 39 P-unit models, we now can explore how many P-unit model neurons show such a triangular structure.
By just looking at the second-order susceptibilities estimated using the noise-split method (first column of \figrefb{fig:modelsusceptcontrasts}) we can readily identify strong triangular patterns in 11 of the 39 model cells (28\,\%, see \figrefb{fig:modelsusceptcontrasts}\,\panel[i]{A}\&\panel[i]{B} for two examples). In another 5 cells (13\,\%) the triangle is much weaker and sits on top of a smooth bump of elevated second-order susceptibility (\figrefb{fig:modelsusceptcontrasts}\,\panel[i]{C} shows an example). The remaining 23 model cells (59\,\%) show no triangle (see \figrefb{fig:modelsusceptcontrasts}\,\panel[i]{D} for an example). This categorization is supported by the peakedness of the nonlinearity, \nli{} \eqnref{eq:nli_equation2}, which quantifies the height of the ridge where the stimulus frequencies add up to the neuron's baseline firing rate. Values above one indicate an elevated ridge. The absence of such a ridge results in values close to one. Indeed, the cells showing only a weak triangle (orange) arise out of values around one and the cells showing strong triangles (red) have consistently \nli{} values exceeding 1.9 (\figrefb{fig:modelsusceptcontrasts}\,\panel[i]{E}). 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).
The \nli{} correlates with the cell's CV of its baseline interspike intervals ($r=-0.59$, $p<0.001$). The lower the cell's CV, the higher the \nli{} value and thus the stronger the triangular structure of its second-order susceptibility. The model cells with the most distinct triangular pattern in their second-order susceptibility are the ones with the lowest CVs, hinting at low intrinsic noise levels. 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}).
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.
\subsection{Weakly nonlinear interactions vanish for higher stimulus contrasts}
As pointed out in the context of stimulation with sine waves (\figrefb{fig:regimes}), the weakly nonlinear regime can only be observed for sufficiently weak stimulus amplitudes. In the model cells we estimated second-order susceptibilities for RAM stimuli with a contrast of 1, 3, and 10\,\%. For the 10\,\% contrast, $10^5$ FFT segments were sufficient for estimating the susceptibility (\figrefb{fig:trialnr}\,\panel[iv]{E}). The estimate converged similarly as in the noise-split configuration (\figrefb{fig:trialnr}\,\panel[i]{E}). For 3\,\% contrast, one million segments were just sufficient (\figrefb{fig:trialnr}\,\panel[iii]{E}) and for 1\,\% contrast even $10^7$ segments were sometimes not enough for a good estimate (\figrefb{fig:trialnr}\,\panel[ii]{E}). This again demonstrates the detrimental effect of smaller signal-to-noise ratios on the estimate of the second-order susceptibility and the advantage of the noise-split method.
The estimates for 1\,\% contrast (\figrefb{fig:modelsusceptcontrasts}\,\panel[ii]{E}), however, were quite similar to the estimates from the noise-split method, corresponding to a stimulus contrast of 0\,\% ($r=0.92$, $p\ll 0.001$). Thus, RAM stimuli with 1\,\% contrast are sufficiently small to not destroy weakly nonlinear interactions by their linearizing effect. 48\,\% of the model cells have a \nli{} value greater than 1.2. \subsection{Weakly nonlinear interactions vanish for higher stimulus contrasts}
As pointed out above, the weakly nonlinear regime can only be observed for sufficiently weak stimulus amplitudes. In the model cells we estimated second-order susceptibilities for RAM stimuli with a contrast of 1, 3, and 10\,\%. The estimates for 1\,\% contrast (\figrefb{fig:modelsusceptcontrasts}\,\panel[ii]{E}) were quite similar to the estimates from the noise-split method, corresponding to a stimulus contrast of 0\,\% ($r=0.97$, $p\ll 0.001$). Thus, RAM stimuli with 1\,\% contrast are sufficiently small to not destroy weakly nonlinear interactions by their linearizing effect. 51\,\% of the model cells have an SI($r$) value greater than 1.2.
At a RAM contrast of 3\,\% the \nli{} values become smaller (\figrefb{fig:modelsusceptcontrasts}\,\panel[iii]{E}). Only 9 cells (23\,\%) have \nli{} values exceeding 1.2. Finally, at 10\,\% the \nli{} values of all cells drop below 1.2, except for two cells (5\,\%, \figrefb{fig:modelsusceptcontrasts}\,\panel[iv]{E}). The cell shown in \subfigrefb{fig:modelsusceptcontrasts}{A} is one of them. Interestingly, the baseline CVs of these two cells (0.16 and 0.23) are not the lowest one (0.11) of the model cells, but their sensitivities are quite low (small response modulations and thus small output noise). At 10\,\% contrast the \nli{} values are no longer correlated with the ones in the noise-split configuration ($r=0.23$, $p=0.15$). To summarize, the regime of distinct nonlinear interactions at frequencies matching the baseline firing rate extends in this set of P-unit model cells to stimulus contrasts ranging from a few percents to about 10\,\%. At a RAM contrast of 3\,\% the SI($r$) values become smaller (\figrefb{fig:modelsusceptcontrasts}\,\panel[iii]{E}). Only 7 cells (18\,\%) have SI($r$) values exceeding 1.2. Finally, at 10\,\% the SI($r$) values of all cells drop below 1.2, except for three cells (8\,\%, \figrefb{fig:modelsusceptcontrasts}\,\panel[iv]{E}). The cell shown in \subfigrefb{fig:modelsusceptcontrasts}{A} is one of them. At 10\,\% contrast the SI($r$) values are no longer correlated with the ones in the noise-split configuration ($r=0.32$, $p=0.05$). To summarize, the regime of distinct nonlinear interactions at frequencies matching the baseline firing rate extends in this set of P-unit model cells to stimulus contrasts ranging from a few percents to about 10\,\%.
\begin{figure}[tp] \begin{figure}[tp]
\includegraphics[width=\columnwidth]{modelsusceptlown} \includegraphics[width=\columnwidth]{modelsusceptlown}
@ -538,18 +536,11 @@ At a RAM contrast of 3\,\% the \nli{} values become smaller (\figrefb{fig:models
\end{figure} \end{figure}
\subsection{Weakly nonlinear interactions can be deduced from limited data} \subsection{Weakly nonlinear interactions can be deduced from limited data}
Estimating second-order susceptibilities reliably requires large numbers (millions) of FFT segments (\figrefb{fig:trialnr}). Electrophysiological measurements, however, suffer from limited recording durations and hence limited numbers of available FFT segments and estimating weakly nonlinear interactions from just ten segments appears futile. The question arises, to what extend such limited-data estimates are still informative? Estimating second-order susceptibilities reliably requires large numbers (millions) of FFT segments (\figrefb{fig:noisesplit}). Electrophysiological measurements, however, suffer from limited recording durations and hence limited numbers of available FFT segments and estimating weakly nonlinear interactions from just a few hundred segments appears futile. The question arises, to what extend such limited-data estimates are still informative?
The second-order susceptibility matrices that are based on only 10 sgements look flat and noisy, lacking the triangular structure \subfigref{fig:modelsusceptcontrasts}{B}. The anti-diagonal ridge, however, when the sum of the stimulus frequencies matches the neuron's baseline firing rate seems to be present whenever the converged estimate shows a clear triangular structure (compare \subfigref{fig:modelsusceptcontrasts}{B} and \subfigref{fig:modelsusceptcontrasts}{a}).
%A comparison of well converged estimates based on millions of segments with estimates based on just ten segments as a worst-case scenario (\subfigrefb{fig:modelsusceptlown}{A\&B}) seems hopeless on a first glance. These estimates using just ten segments look flat and noisy, no triangular structure is visible. However, the anti-diagonal ridge where the stimulus frequencies add up to the neuron's baseline firing rate seems to be present when the converged estimate shows a clear triangular structure.
The \nli{} characterizes the ridgeness of the second-oder susceptibility plane. Comparing it when based on 10 versus one or ten million segments for all 39 model cells (\subfigrefb{fig:modelsusceptlown}{C}) supports this impression. As we have seen in the context of \subfigref{fig:modelsusceptcontrasts}{E} for converged estimates, values greater than one indicate triangular structures in the second-order susceptibility. The \nli{} values based on just ten segments correlate quite well with the ones from the converged estimates for contrasts 1\,\% and 3\,\% ($r=0.9$, $p<0.001$). At a contrast of 10\,\% this correlation is weaker ($r=0.53$, $p<0.001$), because there are only two cells left with \nli{} values greater than one. Despite the good correlations, care has to be taken to set a threshold on the \nli{} values for deciding whether a triangular structure would emerge for a much higher number of segments. Because the estimates are noisier, there could be false positives for a too low threshold. Setting the threshold to 1.6 avoids false positives for the price of a few false negatives.
Trying to predict whether there is a triangular structure in the noise-split configuration is more difficult (\subfigrefb{fig:modelsusceptlown}{D}). The correlations are weaker and at a stimulus contrast of 10\,\% the correlation is no longer significant. One false positive arises at a contrast of 1\,\%, and false positives are absent for higher contrasts. False negatives increase with increasing contrast and at a contrast of 10\,\% all \nli{} values based on 10 segments are so low that just one triangular pattern can be predicted. This makes sense given the absent correlation between \nli{} values estimated at 10\,\% stimulus contrast and the noise-split configuration described above.
Overall, observing \nli{} values greater than at least 1.6, even for a number of FFT segments as low as ten, seems to be a reliable indication for a triangular structure in the second-order susceptibility at the corresponding stimulus contrast. Small stimulus contrasts of 1\,\% are less informative, because of the bad signal-to-noise ratio. Intermediate stimulus contrasts around 3\,\% seem to be optimal, because there, most cells still have a triangular structure in their susceptibility and the signal-to-noise ratio is better. At RAM stimulus contrasts of 10\,\% or higher the signal-to-noise ratio is even better, but only few cells remain with weak triangularly shaped susceptibilities that might be missed as a false positives. Note that increasing the number of segments used for estimating the susceptibilities to 100 or 1000 improves the situation only marginally (not shown). The second-order susceptibility matrices that are based on only 100 sgements look flat and noisy, lacking the triangular structure (\subfigref{fig:modelsusceptlown}{B}). The anti-diagonal ridge, however, where the sum of the stimulus frequencies matches the neuron's baseline firing rate, seems to be present whenever the converged estimate shows a clear triangular structure (compare \subfigref{fig:modelsusceptlown}{B} and \subfigref{fig:modelsusceptlown}{A}). The SI($r$) characterizes the height of the ridge in the second-oder susceptibility plane at the neuron's baseline firing rate $r$. Comparing it when based on 100 versus one or ten million segments for all 39 model cells (\subfigrefb{fig:modelsusceptlown}{C}) supports this impression. As we have seen in the context of \subfigref{fig:modelsusceptcontrasts}{E} for converged estimates, values greater than one indicate triangular structures in the second-order susceptibility. The SI($r$) values based on just one hundred segments correlate quite well with the ones from the converged estimates for contrasts 1\,\% and 3\,\% ($r=0.9$, $p<0.001$). At a contrast of 10\,\% this correlation is weaker ($r=0.38$, $p<0.05$), because there are only three cells left with SI($r$) values greater than 1.2. Despite the good correlations, care has to be taken to set a threshold on the SI($r$) values for deciding whether a triangular structure would emerge for a much higher number of segments. Because at low number of segments the estimates are noisier, there could be false positives for a too low threshold. Setting the threshold to 1.8 avoids false positives for the price of a few false negatives.
Overall, observing SI($r$) values greater than about 1.8, even for a number of FFT segments as low as one hundred, seems to be a reliable indication for a triangular structure in the second-order susceptibility at the corresponding stimulus contrast. Small stimulus contrasts of 1\,\% are less informative, because of their bad signal-to-noise ratio. Intermediate stimulus contrasts around 3\,\% seem to be optimal, because there, most cells still have a triangular structure in their susceptibility and the signal-to-noise ratio is better. At RAM stimulus contrasts of 10\,\% or higher the signal-to-noise ratio is even better, but only few cells remain with weak triangularly shaped susceptibilities that might be missed as a false positives.
\begin{figure*}[tp] \begin{figure*}[tp]
\includegraphics[width=\columnwidth]{dataoverview} \includegraphics[width=\columnwidth]{dataoverview}
@ -601,12 +592,10 @@ Overall, observing \nli{} values greater than at least 1.6, even for a number of
\end{figure*} \end{figure*}
\subsection{Low CVs and weak stimuli are associated with distinct nonlinearity in recorded electroreceptive neurons} \subsection{Low CVs and weak stimuli are associated with distinct nonlinearity in recorded electroreceptive neurons}
Now we are prepared to evaluate our pool of 221 P-units and 47 ampullary afferents recorded in 71 specimen. Now we are prepared to evaluate our pool of 39 P-unit model cells, 159 P-units, and 30 ampullary afferents recorded in 75 specimen of \textit{Apteronotus leptorhynchus}.
\notejb{We need to state the number of trials} \notejb{We need to state the number of trials}
\notejb{We need to know which contrasts} \notejb{We need to know which contrasts}
% All the statements about nonlinear encoding in p-type and ampullary electroreceptor afferents based on single-cell examples shown above are supported by the analysis of our pool of 221 P-units and 47 ampullary afferents recorded in 71 specimen. For comparison across cells we quantify the structure of the second-order susceptibilities by the SI($r$), \eqnref{eq:nli_equation}, that quantifies the expected ridge of the second-order susceptibility at the baseline firing rate (e.g. \subfigref{fig:punit}{G}).
\notejb{\nli{} explanation is needed earlier}
For comparison across cells we summarize the structure of the second-order susceptibilities in a single number, the peakedness of the projected nonlinearity at \fbase{} (\nli{}) \eqnref{eq:nli_equation}, that characterizes the size of the expected peak of the projections of a \suscept{} matrix onto its diagonal at the baseline firing rate (e.g. \subfigref{fig:punit}{G}). \nli{} assumes high values when the peak at \fbase{} is pronounced relative to the median of projections onto the diagonal and is small when there is no distinct peak.
Three P-units stand out with \nli{} values exceeding two, but additional \notejb{XXX} cells have \nli{} values greater than 1.8. Based on our insights from the P-unit models these would have the expected triangular structure in their susceptibilities when estimated with a sufficiently high number of segments. However, we can only speculate how many of the cells with lower \nli{} values are false negatives. \notejb{because also many have been measured at too strong contrasts}. Three P-units stand out with \nli{} values exceeding two, but additional \notejb{XXX} cells have \nli{} values greater than 1.8. Based on our insights from the P-unit models these would have the expected triangular structure in their susceptibilities when estimated with a sufficiently high number of segments. However, we can only speculate how many of the cells with lower \nli{} values are false negatives. \notejb{because also many have been measured at too strong contrasts}.
The \nli{} values of the P-unit population correlate weakly with the CV of the baseline ISIs. Cells with lower baseline CVs tend to have more pronounced peaks in their projections than those that have high baseline CVs (\subfigrefb{fig:dataoverview}{A}). This negative correlation is more pronounced against the CV measured during stimulation (\subfigrefb{fig:dataoverview}{C}). The \nli{} values of the P-unit population correlate weakly with the CV of the baseline ISIs. Cells with lower baseline CVs tend to have more pronounced peaks in their projections than those that have high baseline CVs (\subfigrefb{fig:dataoverview}{A}). This negative correlation is more pronounced against the CV measured during stimulation (\subfigrefb{fig:dataoverview}{C}).

2
punitexamplecell.py
View File

@ -299,7 +299,7 @@ if __name__ == '__main__':
if k == 1: if k == 1:
axc.set_ylabel('') axc.set_ylabel('')
fig.common_yticks(axs[1, :]) fig.common_yticks(axs[1, :])
fig.tag([axs[0, :]], xoffs=-3, yoffs=1) fig.tag([axs[0, :2], axs[0, 2:]], xoffs=-3, yoffs=1)
fig.savefig() fig.savefig()
print() print()