Graphics Processing Unit-Accelerated Code for Computing Second-Order Wiener Kernels and Spike-Triggered Covariance

Abstract

Sensory neuroscience seeks to understand and predict how sensory neurons respond to stimuli. Nonlinear components of neural responses are frequently characterized by the second-order Wiener kernel and the closely-related spike-triggered covariance (STC). Recent advances in data acquisition have made it increasingly common and computationally intensive to compute second-order Wiener kernels/STC matrices. In order to speed up this sort of analysis, we developed a graphics processing unit (GPU)-accelerated module that computes the second-order Wiener kernel of a system's response to a stimulus. The generated kernel can be easily transformed for use in standard STC analyses. Our code speeds up such analyses by factors of over 100 relative to current methods that utilize central processing units (CPUs). It works on any modern GPU and may be integrated into many data analysis workflows. This module accelerates data analysis so that more time can be spent exploring parameter space and interpreting data.

Fig 1. Verified algorithm and degree of speed up.

(A) The second order component of a cell’s response is modeled as a linear weighting of pairwise products in the stimulus. The linear weighting matrix is the Volterra kernel, G⁽²⁾. (B) The GPU code correctly extracts the Wiener kernel from the stimulus-response pair. We simulated a model cell that used a Volterra kernel G⁽²⁾ to respond to a Gaussian input with unit variance. We then used the GPU code to estimate the Weiner kernel K⁽²⁾. Differences between the actual and estimated kernel go to 0 as the number of stimulus-response samples increases. In this case, the Wiener kernel equals the Volterra kernel because the response does not depend on higher order functions of the stimulus. (C-F) Performance of the GPU-enabled module versus the reference CPU implementation. Unless otherwise stated for each simulation parameter sweep, the dataset included 64 responses of 2¹⁹ samples, with inputs of 8 spatial dimensions each. The code extracted filters with 64 temporal offsets. All calculation times are given in seconds; code was run on hardware as described in Methods. Numbers above the GPU datapoints indicate the factor speedup of the GPU module relative to the CPU module. (C) Calculation time comparison for a sweep of the number of samples of each response. (D) Calculation time comparison for a sweep of the number of temporal offsets. (E) Calculation time comparison for a sweep of the number of spatial dimensions. (F) Calculation time comparison for a sweep of the number of responses (e.g. recorded cells in an imaging dataset).

Fig 2. Demonstration of input filter extraction from a model spiking cell.

(A) A model spiking neuron. Two linear filters act on the input, are rectified, and then summed. This determines the probability of firing, which is modeled as a Poisson process. Gaussian random inputs were fed into the model and resulting spikes computed with a mean rate of 0.0045 per time bin. (B) We extracted a spike triggered average from the responses (top). The extracted filter is a mix of the two input filters, and does not accurately represent the underlying linear processing. The outer product of the spike triggered average (bottom) will be used to compute forms of the STC matrix. (C) Three forms of the spike triggered covariance (STC) matrix: raw (C₀), STA subtracted (C₁), and with the STA projected out (C₂). See Methods for details. (D) The eigenspectrum of each STC matrix, with significant eigenvalues highlighted. (E) Eigenvectors of each STC matrix that correspond to significant eigenvalues. In the case of C₂, the STA serves as the second filter. (F) Comparison of the actual input filters (purple and green) with linear combinations of the extracted filters (grey); the eigenvectors span the space of linear input filters.