diff --git a/modelsusceptlown.py b/modelsusceptlown.py index 9ac9081..2747110 100644 --- a/modelsusceptlown.py +++ b/modelsusceptlown.py @@ -90,10 +90,10 @@ def plot_si_diags(ax, s, data, alphax, alphay, xthresh, ythresh, cell_name): ax.set_xlim(0, 9) n = datax[0, 'nsegs'] if alphax == 0: - ax.set_xlabel(f'SI, $c=0$, $N=10^{np.log10(n):.0f}$') + ax.set_xlabel(f'SI($r$), $c=0$, $N=10^{np.log10(n):.0f}$') else: - ax.set_xlabel(f'SI, $N=10^{np.log10(n):.0f}$') - ax.set_ylabel('SI, $N=100$') + ax.set_xlabel(f'SI($r$), $N=10^{np.log10(n):.0f}$') + ax.set_ylabel('SI($r$), $N=100$') ax.set_xticks_delta(4) ax.set_yticks_delta(4) ax.set_minor_xticks_delta(1) diff --git a/nonlinearbaseline.tex b/nonlinearbaseline.tex index 40fcfc4..c801144 100644 --- a/nonlinearbaseline.tex +++ b/nonlinearbaseline.tex @@ -266,14 +266,14 @@ We here analyze nonlinear responses in two types of primary electroreceptor affe \begin{figure*}[t] \includegraphics[width=\columnwidth]{lifsuscept} - \caption{\label{fig:lifresponse} First- and second-order response functions of the leaky integrate-and-fire model. \figitem{A} Magnitude of the first-order (linear) response function $|\chi_1(f)|$, also known as the ``gain'' function, quantifies the response amplitude relative to the stimulus amplitude, both measured at the same stimulus frequency. \figitem{B} Magnitude of the second-order (non-linear) response function $|\chi_2(f_1, f_2)|$ quantifies the response at the sum of two stimulus frequencies. For linear systems, the second-order response function is zero, because linear systems do not create new frequencies and thus there is no response at the sum of the two frequencies. The plots show the analytical solutions from \citet{Lindner2001} and \citet{Voronenko2017} with $\mu = 1.1$ and $D = 0.001$. Note that the leaky integrate-and-fire model is formulated without dimensions, frequencies are given in multiples of the inverse membrane time constant.} + \caption{\label{fig:lifsuscept} First- and second-order response functions of the leaky integrate-and-fire model. \figitem{A} Magnitude of the first-order (linear) response function $|\chi_1(f)|$, also known as the ``gain'' function, quantifies the response amplitude relative to the stimulus amplitude, both measured at the same stimulus frequency. \figitem{B} Magnitude of the second-order (non-linear) response function $|\chi_2(f_1, f_2)|$ quantifies the response at the sum of two stimulus frequencies. Because frequencies can also be negative, the sum is actually a difference in the upper left and lower right quadrants, since $f_2 + f_1 = f_2 - |f_1|$ for $f_1 < 0$. For linear systems, the second-order response function is zero, because linear systems do not create new frequencies and thus there is no response at the sum of the two frequencies. The plots show the analytical solutions from \citet{Lindner2001} and \citet{Voronenko2017} with $\mu = 1.1$ and $D = 0.001$. Note that the leaky integrate-and-fire model is formulated without dimensions, frequencies are given in multiples of the inverse membrane time constant.} \end{figure*} -We like to think about signal encoding in terms of linear relations with unique mapping of a given input value to a certain output of the system under consideration. Indeed, such linear methods, for example the transfer function or first-order susceptibility shown in figure~\ref{fig:lifresponse}, have been widely and successfully applied to describe and predict neuronal responses and are an invaluable tools to characterize neural systems \citep{Eggermont1983,Borst1999}. Nonlinear mechanisms, on the other hand, are key on different levels of neural processing. Deciding for one action over another is a nonlinear process on the systemic level. On the cellular level, spiking neurons are inherently nonlinear. Whether an action potential is elicited depends on the membrane potential to exceed a threshold \citep{Hodgkin1952, Koch1995}. Because of such nonlinearities, understanding and predicting neuronal responses to sensory stimuli is in general a difficult task. +We like to think about signal encoding in terms of linear relations with unique mapping of a given input value to a certain output of the system under consideration. Indeed, such linear methods, for example the transfer function or first-order susceptibility shown in \subfigref{fig:lifsuscept}{A}, have been widely and successfully applied to describe and predict neuronal responses and are an invaluable tools to characterize neural systems \citep{Eggermont1983,Borst1999}. Nonlinear mechanisms, on the other hand, are key on different levels of neural processing. Deciding for one action over another is a nonlinear process on the systemic level. On the cellular level, spiking neurons are inherently nonlinear. Whether an action potential is elicited depends on the membrane potential to exceed a threshold \citep{Hodgkin1952, Koch1995}. Because of such nonlinearities, understanding and predicting neuronal responses to sensory stimuli is in general a difficult task. The transfer function that describes the linear properties of a system, is the first-order term of a Volterra series. Higher-order terms successively approximate nonlinear features of a system \citep{Rieke1999}. Second-order kernels have been used in the time domain to predict visual responses in catfish \citep{Marmarelis1972}. In the frequency domain, second-order kernels are known as second-order response functions or susceptibilities. Nonlinear interactions of two stimulus frequencies generate peaks in the response spectrum at the sum and the difference of the two. Including higher-order terms of the Volterra series, the nonlinear nature of mammalian visual systems \citep{Victor1977, Schanze1997}, auditory responses in the Torus semicircularis of frogs \citep{Aertsen1981}, locking in chinchilla auditory nerve fibers \citep{Temchin2005}, spider mechanoreceptors \citep{French2001}, and bursting responses in paddlefish \citep{Neiman2011} have been demonstrated. -Noise in nonlinear systems, however, linearizes the system's response properties \citep{Yu1989, Chialvo1997}. Also, in the limit to small stimuli, nonlinear systems can be well described by linear response theory \citep{Roddey2000, Doiron2004, Rocha2007, Sharafi2013}. With increasing stimulus amplitude, the contribution of the second-order kernel of the Volterra series becomes more relevant. For these weakly nonlinear responses analytical expressions for the second-order susceptibility have been derived for leaky-integrate-and-fire (LIF) \citep{Voronenko2017} and theta model neurons \citep{Franzen2023}. In the suprathreshold regime, in which the LIF generates a baseline firing rate in the absence of an external stimulus, the linear response function has a peak at the baseline firing rate and its harmonics (\subfigrefb{fig:lifresponse}{A}) and the second-order susceptibility shows very distinct ridges of elevated nonlinear responses, exactly where one of two stimulus frequencies equals or both frequencies add up to the neuron's baseline firing rate (\subfigrefb{fig:lifresponse}{B}). In experimental data, such structures in the second-order susceptibility have not been reported yet. +Noise in nonlinear systems, however, linearizes the system's response properties \citep{Yu1989, Chialvo1997}. Also, in the limit to small stimuli, nonlinear systems can be well described by linear response theory \citep{Roddey2000, Doiron2004, Rocha2007, Sharafi2013}. With increasing stimulus amplitude, the contribution of the second-order kernel of the Volterra series becomes more relevant. For these weakly nonlinear responses analytical expressions for the second-order susceptibility have been derived for leaky-integrate-and-fire (LIF) \citep{Voronenko2017} and theta model neurons \citep{Franzen2023}. In the suprathreshold regime, in which the LIF generates a baseline firing rate in the absence of an external stimulus, the linear response function has a peak at the baseline firing rate and its harmonics (\subfigrefb{fig:lifsuscept}{A}) and the second-order susceptibility shows very distinct ridges of elevated nonlinear responses, exactly where one of two stimulus frequencies equals or both frequencies add up to the neuron's baseline firing rate (\subfigrefb{fig:lifsuscept}{B}). In experimental data, such structures in the second-order susceptibility have not been reported yet. Here we search for such weakly nonlinear responses in electroreceptors of the two electrosensory systems of the wave-type electric fish \textit{Apteronotus leptorhynchus}, i.e. the tuberous (active) and the ampullary (passive) electrosensory system. The p-type electroreceptor afferents of the active system (P-units) are driven by the fish's high-frequency, quasi-sinusoidal electric organ discharges (EOD) and encode disturbances of it \citep{Bastian1981a}. The electroreceptors of the passive system are tuned to lower-frequency exogeneous electric fields such as caused by muscle activity of prey \citep{Kalmijn1974}. As different animals have different EOD-frequencies, being exposed to stimuli of multiple distinct frequencies is part of the animal's everyday life \citep{Benda2020,Henninger2020} and weakly nonlinear interactions may occur in the electrosensory periphery. In communication contexts \citep{Walz2014, Henninger2018} the EODs of interacting fish superimpose and lead to periodic amplitude modulations (AMs or beats) of the receiver's EOD. Nonlinear mechanisms in P-units, enable encoding of AMs in their time-dependent firing rates \citep{Bastian1981a, Walz2014, Middleton2006, Barayeu2023}. When multiple animals interact, the EOD interferences induce second-order amplitude modulations referred to as envelopes \citep{Yu2005, Fotowat2013, Stamper2012Envelope} and saturation nonlinearities allow also for the encoding of these in the electrosensory periphery \citep{Savard2011}. Field observations have shown that courting males were able to react to the extremely weak signals of distant intruding males despite the strong foreground EOD of the nearby female \citep{Henninger2018}. Weakly nonlinear interactions at particular combinations of signals can be of immediate relevance in such settings as they could boost detectability of the faint signals \citep{Schlungbaum2023}. @@ -323,11 +323,11 @@ P-units fire action potentials probabilistically phase-locked to the self-genera 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{chi2}, quantifies for each combination of two stimulus frequencies $f_1$ and $f_2$ the amplitude and phase of the stimulus-evoked response at the sum $f_1 + f_2$ (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 $f_1$ and $f_2$ or their sum $f_1 + f_2$ 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 $f_1$ and $f_2$ the amplitude and phase of the stimulus-evoked response at the sum $f_1 + f_2$ (and also the difference, \subfigrefb{fig:lifsuscept}{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 $f_1$ and $f_2$ or their sum $f_1 + f_2$ 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:lifsuscept}{B}). For the example P-unit, we observe a ridge of elevated second-order susceptibility for the low RAM contrast at $f_1 + f_2 = r$ (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. -Overall we observed in 17\,\% of the 159 P-units ridges where the stimulus frequencies add up to the unit's baseline firing rate. Two more examples are shown in \subfigref{fig:punit}{H}. However, we never observed the full triangular structure expected from theory (\subfigrefb{fig:lifresponse}{B}). In all other P-units, we did not observe any structure in the second-order susceptibility (\subfigrefb{fig:punit}{I}). +Overall we observed in 17\,\% of the 159 P-units ridges in the second-order susceptibility where the stimulus frequencies add up to the unit's baseline firing rate. Two more examples are shown in \subfigref{fig:punit}{H}. However, we never observed the full triangular structure expected from theory (\subfigrefb{fig:lifsuscept}{B}). In all other P-units, we did not observe any structure in the second-order susceptibility (\subfigrefb{fig:punit}{I}). \subsection{Ampullary afferents exhibit strong nonlinear interactions} @@ -380,7 +380,7 @@ At a RAM contrast of 3\,\% the SI($r$) values become smaller (\figrefb{fig:model \begin{figure}[tp] \includegraphics[width=\columnwidth]{modelsusceptlown} - \caption{\label{fig:modelsusceptlown}Inferring the triangular structure of the second-order susceptibility from limited data. \figitem{A} Reliably estimating the structure of the second-order susceptibility requires a high number of FFT segments $N$ in the order of one or even ten millions. As an example, susceptibilities of the model cell ``2012-12-21-ak-invivo-1'' (baseline firing rate of 157\,Hz, CV=0.15) are shown for the noise-split configuration ($c=0$\,\%) and RAM stimulus contrasts of $c=1$, $3$, and $10$\,\% as indicated. For contrasts below $10$\,\% this cell shows a nice triangular pattern in its susceptibilities, quite similar to the introductory example of a LIF in \figrefb{fig:lifresponse}. \figitem{B} However, with limited data of $N=100$ trials the susceptibility estimates are noisy and show much less structure, except for the anti-diagonal at the cell's baseline firing rate. The SI($r$) quantifies the height of this ridge where the two stimulus frequencies add up to the neuron's baseline firing rate. \figitem{C} Correlations between the estimates of SI($r$) based on 100 FFT segments (density to the right) with the converged ones based on one or ten million segments at a given stimulus contrast for all $n=39$ model cells. The black circle marks the model cell shown in \panel{A} and \panel{B}. The black diagonal line is the identity line and the dashed line is a linear regression. The correlation coefficient and corresponding significance level are indicated in the top left corner. The thin vertical line is a threshold at 1.2, the thin horizontal line a threshold at 1.8. The number of cells within each of the resulting four quadrants denote the false positives (top left), true positives (top right), true negatives (bottom left), and false negatives (bottom right) for predicting a triangular structure in the converged susceptibility estimate from the estimates based on only 100 segments.} + \caption{\label{fig:modelsusceptlown}Inferring the triangular structure of the second-order susceptibility from limited data. \figitem{A} Reliably estimating the structure of the second-order susceptibility requires a high number of FFT segments $N$ in the order of one or even ten millions. As an example, susceptibilities of the model cell ``2012-12-21-ak-invivo-1'' (baseline firing rate of 157\,Hz, CV=0.15) are shown for the noise-split configuration ($c=0$\,\%) and RAM stimulus contrasts of $c=1$, $3$, and $10$\,\% as indicated. For contrasts below $10$\,\% this cell shows a nice triangular pattern in its susceptibilities, quite similar to the introductory example of a LIF in \figrefb{fig:lifsuscept}. \figitem{B} However, with limited data of $N=100$ trials the susceptibility estimates are noisy and show much less structure, except for the anti-diagonal at the cell's baseline firing rate. The SI($r$) quantifies the height of this ridge where the two stimulus frequencies add up to the neuron's baseline firing rate. \figitem{C} Correlations between the estimates of SI($r$) based on 100 FFT segments (density to the right) with the converged ones based on one or ten million segments at a given stimulus contrast for all $n=39$ model cells. The black circle marks the model cell shown in \panel{A} and \panel{B}. The black diagonal line is the identity line and the dashed line is a linear regression. The correlation coefficient and corresponding significance level are indicated in the top left corner. The thin vertical line is a threshold at 1.2, the thin horizontal line a threshold at 1.8. The number of cells within each of the resulting four quadrants denote the false positives (top left), true positives (top right), true negatives (bottom left), and false negatives (bottom right) for predicting a triangular structure in the converged susceptibility estimate from the estimates based on only 100 segments.} \end{figure} \subsection{Weakly nonlinear interactions can be deduced from limited data} @@ -458,7 +458,7 @@ The population of ampullary cells is more homogeneous, with generally lower base \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 sine wave stimuli \notejg{with finite amplitudes, needed?} 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. +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 sine wave stimuli \notejg{with finite amplitudes, needed?} 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:lifsuscept}{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{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 $|\chi_2|(f_1, f_2)$, 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}). @@ -551,9 +551,9 @@ When stimulated with band-limited white noise stimuli, neuronal activity is modu 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 numerically according to \begin{equation} \label{fourier} - X(\omega) = \sum_{k=0}^{n-1} \, x_k e^{- i \omega k} + X(\omega) = \sum_{k=0}^{N-1} \, x_k e^{- i \omega k} \end{equation} -for $n=512$ long segments of $T=n \Delta t = 256$\,ms duration with no overlap, resulting in a spectral resolution of about 4\,Hz. Note, that for a real Fourier integral a factor $\Delta t$ is missing. For simplicity we use angular frequencies $\omega=2\pi f$ instead of frequencies $f$. +for $N=512$ long segments of $T=N \Delta t = 256$\,ms duration with no overlap, resulting in a spectral resolution of about 4\,Hz. Note, that for a real Fourier integral a factor $\Delta t$ is missing. 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. @@ -562,14 +562,14 @@ In the simulations we generated for each trial a new realization of the noise st The power spectrum of the stimulus $s(t)$ was estimated as \begin{equation} \label{powerss} - S_{ss}(\omega) = \frac{\Delta t}{n} \langle S(\omega) S^*(\omega) \rangle + S_{ss}(\omega) = \frac{\Delta t}{N} \langle S(\omega) S^*(\omega) \rangle \end{equation} -with $S^*$ being the complex conjugate of $S$ and $\langle \cdot \rangle$ denoting the average over the FFT segments. The factor in front is $\Delta t^2$ from the missing integration factors of the two Fourier transforms, \eqnref{fourier}, divided by $T=n \Delta t$, needed to make this a proper power spectral density. +with $S^*$ being the complex conjugate of $S$ and $\langle \cdot \rangle$ denoting the average over the FFT segments. The factor in front is $\Delta t^2$ from the missing integration factors of the two Fourier transforms, \eqnref{fourier}, divided by $T=N \Delta t$, needed to make this a proper power spectral density. 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{crossxs} - S_{xs}(\omega) = \frac{\Delta t}{n} \langle X(\omega) S^*(\omega) \rangle + S_{xs}(\omega) = \frac{\Delta t}{N} \langle X(\omega) S^*(\omega) \rangle \end{equation} The first-order susceptibility (transfer function) @@ -582,9 +582,9 @@ was then computed from $S_{xs}(\omega)$ and $ S_{ss}(\omega)$. We report $\chi_{ The second-order cross-spectrum \begin{equation} \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}) + 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} -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$. +quantifies nonlinear interactions that generate responses at the sum and difference (for negative $\omega_1$ or $\omega_2$) evoked by two stimulus frequencies $\omega_1$ and $\omega_2$. The second-order susceptibility \begin{equation} \label{chi2}