All data analyses were performed offline using custom-written codes in MATLAB software (MathWorks). Spike times were defined as the times at which the recorded neural signal crossed a given threshold value from below. A binary sequence was constructed from the spike times by discretizing time into bins of 0.1-ms width and setting the content of a given bin to 10,000 if a spike occurred within it or 0 otherwise. The time-varying FR was obtained by low pass filtering the binary sequence using a second-order Butterworth filter with cutoff frequencies 0.2, 0.35, 0.75, 1.5, 2.5, and 3.5 Hz for envelope frequencies 0.05, 0.1, 0.2, 0.5, 0.75, and 1 Hz, respectively, as done previously (39). Neural sensitivity for a given sinusoidal frequency was then measured as the ratio between the amplitude of the sinusoidal envelope stimulus as extracted by the dipole and the amplitude of the neural FR modulation. We also quantified neural sensitivity to the adaptation stimulus using linear systems identification techniquesG(f)=Prs(f)Pss(f)(1)where Prs(f) is the cross-spectral density between the adaptation stimulus and the binary sequence, and Pss(f) is the adaptation stimulus power spectral density as a function of frequency f. The cross-spectral density was computed in MATLAB using nfft = 1024 × 260, noverlap = nfft/2 with multitaper techniques with eight Slepian tapers using the functions “cpsd” and “pwelch” in MATLAB, respectively. Behavioral responses were measured as changes in the animal’s EOD frequency as done previously (12, 13, 39, 42). Specifically, the sequence of inverse intervals between consecutive EOD zero crossings was used to compute the instantaneous EOD frequency, which was then interpolated to match the sampling rate of the envelope stimulus and further low pass filtered (second-order Butterworth filter with 0.05-Hz cutoff frequency). Behavioral sensitivity for a given sinusoidal frequency was then measured as the ratio between the amplitude of the sinusoidal envelope stimulus as extracted by the dipole and the amplitude of the sinusoidal EOD frequency modulation. The exponents were obtained by fitting a power law to the behavioral and neural sensitivities as a function of frequency, respectively. The whiteness index was measured by taking the normalized area under the power spectral density curve using a trapezoidal method and dividing by the maximum normalized area to achieve a value between 0 and 1 as done previously (12, 39). The matching index between the behavioral sensitivity and the adapting stimulus was computed using 1 − |αstim − αbehavior| as done previously (12), where αstim and αbehavior are the best-fit power law exponents of the adapting stimulus (i.e., either −2 or 0) and of the behavioral sensitivity, respectively. We note that, as not all neurons could be held during the full 2-hour-long adaptation stimulus presentation, we always compared the average whiteness index values between the first and last three repetitions of the adaptation stimulus. Spectrograms were computed by estimating the power spectral density of the spiking response for every 100 s during stimulation using multitaper techniques with eight Slepian functions. Power spectra were then stacked for visualization. The time-dependent whiteness index was computed as described above for every 100 s during stimulation.

Theory posits that the response power spectral density Prr(f) is related to the gain |G(f)|, the stimulus power spectral density Pss(f), and the power spectral density of the trial-to-trial variability in the neural response P0(f) by the following equation (17, 44)Prr(f)P0(f)+G(f)2Pss(f)(2)

Hence, we considered both the contributions of the variability as well as the tuning function |G(f)| to predict the neural response to stimulation. Trial-to-trial variability was estimated from responses to at least three repeated presentations of the adaptation stimulus each lasting 100 s (17).