Inferring entire spiking activity from local field potentials

Extracellular recordings are typically analysed by separating them into two distinct signals: local field potentials (LFPs) and spikes. Previous studies have shown that spikes, in the form of single-unit activity (SUA) or multiunit activity (MUA), can be inferred solely from LFPs with moderately good accuracy. SUA and MUA are typically extracted via threshold-based technique which may not be reliable when the recordings exhibit a low signal-to-noise ratio (SNR). Another type of spiking activity, referred to as entire spiking activity (ESA), can be extracted by a threshold-less, fast, and automated technique and has led to better performance in several tasks. However, its relationship with the LFPs has not been investigated. In this study, we aim to address this issue by inferring ESA from LFPs intracortically recorded from the motor cortex area of three monkeys performing different tasks. Results from long-term recording sessions and across subjects revealed that ESA can be inferred from LFPs with good accuracy. On average, the inference performance of ESA was consistently and significantly higher than those of SUA and MUA. In addition, local motor potential (LMP) was found to be the most predictive feature. The overall results indicate that LFPs contain substantial information about spiking activity, particularly ESA. This could be useful for understanding LFP-spike relationship and for the development of LFP-based BMIs.

noise and takes into account the spike contribution from small neurons. A number of studies have demonstrated that ESA can achieve better accuracy and reliability in measuring evoked responses (e.g. receptive field) in the visual cortex of monkeys while receiving various visual stimuli 14,15,[20][21][22] . ESA has also been shown to yield more accurate and robust decoding of hand kinematics compared to SUA and MUA from three monkeys performing different tasks 23 .
Understanding the relationship between the LFP and spikes plays a critical role in addressing many issues in neuroscience research, such as neuronal functional organisation and connectivity, neuronal communication between different brain areas, cell assembly formation, neural coding and information processing [24][25][26][27] . In addition, it is also relevant for brain-machine interface (BMI), for example, measuring behavioural task-related information within different neural signals 28,29 and extracting indirect spiking information features from LFPs in biofeedback based BMI 30 . As LFPs are thought to represent mainly the input to local neuronal networks, while spiking activity represents the output from local neuronal networks 31 , it is conceivable to relate the two signals based on system identification-based inference approach. Previous studies have shown that SUA 30,32,33 and MUA 28,34 can be inferred solely from LFPs with moderately good accuracy. So far, however, there has been no study that investigates the relationship between the LFPs and ESA.
In light of this, our present study aims to address the above-mentioned issue by using extracellular recordings from the motor cortex area of three macaque monkeys while performing two different tasks: point-to-point task and instructed delay reach-to-grasp task. We assess quantitatively how well ESA can be inferred from LFPs with multivariate multiple linear regression (MLR) method. We then analyse which features within LFPs are highly predictive of ESA. Furthermore, we compare the inference accuracy of ESA with those of SUA and MUA. We also investigate the impact of the number of LFP channels on the inference performance. Finally, we examine LFP channel importance for the inference performance.

Results
LFP feature informativeness for ESA inference. We evaluated and compared the informativeness of six different LFP features: the smoothed time-domain amplitude of LFP known as local motor potential (LMP) and average power spectra within five frequency bands (theta, delta, alpha, beta, and gamma). The informativeness of LFP features for the inference of ESA was measured using two approaches. First, we quantified and compared the inference performance of LFP features (in terms of average CC) with an MLR model that was fit independently (separately) on each LFP feature ( p = 96 ). Second, we took the absolute value of coefficients (weights) of an MLR model that was fit simultaneously on all LFP features combined ( p = 576 ). We then compared the average coefficients associated with each LFP feature. Figure 1a-c compare the average CC of six LFP features from monkey I, L, and N, respectively. Across 960 cases (10 blocks × 96 ESA channels) from each subject, LMP was found to consistently yield the highest inference performance (largest average CC) while power spectra in the alpha band (shortly referred to as alpha) consistently yielded the lowest inference performance (a-c) Boxplot comparison of average CC across LFP features from monkey I, L, and N, respectively. Asterisks indicate LFP features whose inference performance differed significantly from that of LMP (***p < 0.001 ). The horizontal lines and circles within the boxes indicate the median and mean, respectively. The boxes represent interquartile range (IQR) from 25th to 75th percentiles. The whisker extends 1.5 times the IQR. (d-f) Bar plot comparison of average coefficients (i.e. weights) across LFP features from monkey I, L, and N, respectively.
Comparison of inference performance: ESA vs SUA vs MUA. We next compared the inference of ESA from LFP features to the inference of other types of spiking (SUA and MUA) from LFP features. We found the average of ESA as follows: 0.70 ± 0.02, 0.75 ± 0.02, and 0.55 ± 0.02 for monkey I, L, and N, respectively. The average CC of SUA (MUA) were 0.33 ± 0.01, 0.59 ± 0.02, and 0.41 ± 0.01 ( 0.56 ± 0.02, 0.76 ± 0.02, and 0.41 ± 0.02) for monkey I, L, and N, respectively. The average CC of ESA was statistically significantly larger than that of SUA (MUA) in three (two) subjects as can be seen from Fig. 2a-c. The higher inference performance of ESA was consistently observed across long-term recording sessions from monkey I as shown in Fig. 2d. The average CC of ESA, SUA, and MUA were 0.65 ± 0.01, 0.32 ± 0.01, and 0.48 ± 0.01, respectively. The relative inference performance of ESA was 2.03 and 1.35 times higher compared to those of SUA and MUA, respectively. Figure 2e displays boxplot of ESA inference that was statistically significant different than those of SUA and MUA (***p < 0.001 ). Snippet examples of actual and inferred ESA from monkey I (channel 5), L (channel 45), and N (channel 62) are illustrated in Fig. 2f-h. Overall, the order of spiking activity from highest to lowest inference performance was ESA> MUA > SUA.
The above results were quantified in terms of average CC, which is a translation and scale invariant metric that only assesses the shape similarity between actual and inferred signals. In the presence of a large but constant inference error, the CC value can be high (implying falsely good accuracy). Therefore, as an additional metric for the inference performance, we also used RMSE that represents an average magnitude of inference error between actual and inferred signals. Similar to that of CC metric, the inference performance of ESA across three subjects and long-term recording sessions were consistently better (smaller RMSE) than those of SUA and MUA (see Fig. 3). The average RMSE of ESA, SUA, and MUA for monkey I were 0.70 ± 0.01, 0.94 ± 0.01, and 0.82 ± 0.01, respectively. The average RMSE of ESA, SUA, and MUA for monkey L were 0.64 ± 0.02, 0.79 ± 0.01, and 0.62 ± 0.02, respectively. The average RMSE of ESA, SUA, and MUA for monkey N were 0.81 ± 0.01, 0.90 ± 0.01, and 0.89 ± 0.01, respectively. Figure 3d shows inference performance comparison among different types of spiking activity over long-term recording sessions. The average RMSE of ESA, SUA, and MUA were 0.76 ± 0.01, 0.95 ± 0.01, and 0.87 ± 0.01, respectively. We found statistically significant difference in average RMSE among ESA, SUA, and MUA (see Fig. 3e).
We were also interested in determining whether these findings could hold true when using different algorithms. We repeated the same procedures but using two different deep learning algorithms: multilayer perceptron (MLP) and long short-term memory (LSTM). More detailed description of MLP and LSTM used in this study along with their hyperparameter configurations can be found in Supplementary Information. Comparison of inference performance among ESA, SUA, and MUA using MLP and LSTM can be seen in Supplementary Figure 3 and Supplementary Figure 4, respectively. We observed similar trend where on average the inference performance of ESA was higher than those of SUA and MUA across subjects and long-term recording sessions. We then compared ESA inference performance under MLR, MLP, and LSTM algorithms. We found that the inference performance of MLR ( 0.70 ± 0.02, 0.75 ± 0.02, and 0.55 ± 0.02 ) was lower than those of MLP ( 0.73 ± 0.02, 0.75 ± 0.02, and 0.62 ± 0.02) and LSTM ( 0.73 ± 0.02, 0.76 ± 0.01, and 0.63 ± 0.02) across three subjects (monkey I, L, and N, respectively). Relative to MLR, the deep learning algorithms yielded an average performance improvement of 4.15%, 0.53%, and 12.97% for monkey I, L, and N, respectively. There was statistically significant difference in inference performance among MLR, MLP, and LSTM within single recording session across three subjects. However, in the case of long-term recording sessions from monkey I, the difference in inference performance was not statistically significant as shown in Supplementary Figure 5 www.nature.com/scientificreports/ Impact of number of LFP channels on inference performance. Next, we investigated the impact of different number of LFP channels on inference performance for different types of spiking activity. Performance comparison among ESA, SUA, and MUA across three subjects are plotted in Fig. 4a-c. The results across subjects showed that the inference performance of ESA was always higher than those of SUA and MUA regardless of the number of LFP channels. In addition, the inference performance of ESA, SUA, and MUA improved with the increasing number of LFP channels. Initially, the inference performance improved quickly but after a certain point, it only improved marginally (i.e. reaching plateau). However, the rate at which the inference performance reached a plateau was different among ESA, SUA, and MUA. Further analysis revealed that the inference performance of ESA reached a plateau quicker (in fewer number of LFP channels) than those of SUA and MUA. Averaging across three subjects, the inference performance of ESA, SUA, and MUA reached 90% of their maximum performance when using 35, 42, and 42 LFP channels, respectively. Furthermore, we examined whether there is relationship between LFP interchannel correlation and the rate at which the inference performance reached a plateau. Comparison of LFP interchannel correlation across three subjects is plotted in Fig. 4d-f. The average LFP interchannel correlation for monkey I, L, and N were 0.63, 0.73, and 0.41, respectively. We found that the higher the LFP interchannel correlation, the faster the rate at which the inference performance of ESA, SUA, and MUA reached their plateau. The inference of ESA, SUA, and MUA reached their 90% of maximum performance when using 25, 30, and 35 LFP channels, respectively, for monkey I (10, 25, and 20 LFP channels, respectively, for monkey L; 70, 70, and 70 LFP channels, respectively, for monkey N).
LFP channel importance for ESA inference. Next, we assessed the LFP channel importance for the inference of all 96 ESA channels. The LFP channel importance was quantified using two metrics: average CC www.nature.com/scientificreports/ and average coefficient (see Methods section). It was then sorted based on the distance between all possible pairs of electrodes associated with the ESA channel output and LFP channel input. Figure 5a,c,e show scatter plots of average CC across interelectrode distances for monkey I, L, and N, respectively. We found statistically significant decreasing trends of average CC with the increase of interelectrode distance in two out of three subjects ( p < 0.001 for monkey I and N). Despite the decreasing trend, the shortest interelectrode distance did not necessarily yield the highest average CC. Similarly, the farthest the interelectrode distance also did not necessarily yield the lowest average CC. The average (maximum) CC for monkey I, L, and N was 0.40 (0.73), 0.38 (0.81), and 0.12 (0.66), respectively. Figure 5b,d,f show examples of LFP channel importance score (average CC) mapped onto 10-by-10 heatmap grids for ESA inference of channel 5 (monkey I), channel 45 (monkey L), and channel 62 (monkey N), respectively. When using average coefficient for quantifying LFP channel importance, we also observed decreasing trends (negative slope) in two out of three subjects. However, only one case (monkey N) was found to be statistically significant ( p < 0.001 ) as can be seen from Fig. 6a,c,e. Representative examples of LFP channel importance score for ESA inference of channel 5 (monkey I), channel 45 (monkey L), and channel 62 (monkey N) are plotted as heatmaps of 10 × 10 grid corresponding to the Utah array configuration, as shown in Fig. 6b,d,f.
Lastly, we also examined LFP channel importance score for the inference of SUA and MUA using the same approach as for ESA. Since each channel can have more than one SUA, we averaged the LFP channel importance score for SUA inference from the same channel. The results of LFP channel importance score for SUA inference using average CC and average coefficient metrics are visualised in Supplementary Figures 7 and 8, respectively. As for MUA inference, LFP channel importance score for MUA inference in terms of average CC and average coefficient are visualised in Supplementary Figures 9 and 10, respectively. From a total of six cases (3 subjects and two metrics), we found decreasing trends (negative slope) in four cases (three cases being statistically significant) for SUA and five cases (three cases being statistically significant) for MUA.

Discussion
The present study investigates the relationship between local field potential (LFP) and entire spiking activity (ESA) by asking whether we can infer ESA solely from LFPs. In doing so, we firstly examined which feature within LFPs carries the highest information (i.e. most predictive) about ESA. Our experimental results revealed that LMP emerged as the most predictive feature. The power from high-frequency band (gamma) also showed high inference performance, especially when using 300 Hz cut-off frequency (see Supplementary Figure 2), although its inference performance was lower than that of LMP. On the other hand, the power from intermediate frequency bands carried only a little information about ESA. These findings were consistently observed across subjects performing different tasks. A considerable number of previous studies have reported that power modulation within high-frequency band ( > 30Hz) were highly informative, whereas the intermediate frequency bands were little informative about spiking activity (SUA and MUA) 26,[32][33][34][35] . It is believed that the high-frequency band of LFPs are thought to reflect synchronised spiking activity arising from local neuronal population 32 . Our study confirms www.nature.com/scientificreports/ and extends the previous findings by showing that apart from power within high-frequency LFPs, LMP also contained substantial information about spiking activity. LMP even yielded the highest inference performance in both CC and RMSE metrics across subjects. Even though LMP has been frequently used as a feature for decoding behavioural tasks [36][37][38] , the use of LMP for inferring spiking activity has not been investigated. A rather similar feature to LMP is low-frequency LFP (lf-LFP) which is obtained by low-pass filtering the broadband LFPs 28,30 .
The frequency band of lf-LFP is more clearly separated than that of LMP as low-pass filter results in better roll-off and stopband attenuation. On the other hand, LMP is obtained by a moving average filter which is simple, fast and yields less random white noise while maintaining the smoothing of LFP signals well 39 .
Given the highly informativeness of LMP, one may ask what the biophysical origin of LMP is. Unfortunately, the answer to this question remains to be established. It has been speculated that LMP is related to firing rate modulation of neuronal population near the recording electrode 40 or distant neuronal population activity that is connected to the recording site 26,41 . It is also possible that LMP reflects evoked LFPs in the motor cortex in response to movement-related tasks 32,40 . A prior study suggested that LMP plays an important role not only in the execution but also in the preparation of movement (e.g. anticipation) 42 . This aligns well with our our results in which LMP showed the highest inference performance on movement experiment with and without preparatory delay interval (dataset II and I, respectively).
Next, we evaluated and compared the inference performance of ESA with those of SUA and MUA using LMP feature and multivariate multiple linear regression (MLR). Results across recording sessions and subjects showed that the inference performance of ESA was consistently and significantly higher than those of SUA and MUA. The same trends were also observed when using two deep learning based inference algorithms (MLP and LSTM). Previous studies have only conducted the inference of SUA and MUA from LFPs 28,30,[32][33][34] . Therefore, to the best of our knowledge, our study is the first that evaluates the inference of ESA from LFPs and systematically compares its performance to those of SUA and MUA. LFP has been reported to have spatial reach within around 200-400 μm 21,43,44 with other studies suggesting a broader spatial reach up to a few millimetres 45,46 . Our results may indicate that LFP, which has relatively broad spatial reach, is more related to population spiking activity than single-unit activity. ESA which reflects population spiking activity contains richer spiking information compared to SUA and MUA. Both SUA and MUA are typically obtained by using a threshold-based technique which is prone to false-positive or false-negative spike detection. If we set the threshold value too high, we could miss true but lower amplitude spikes; on the contrary, if we set the threshold too low, we may detect some background noises as spikes. In contrast to SUA and MUA, ESA does not use thresholding, instead, it employs full-rectification and low pass filter which are rather insensitive to random high-frequency noise and are able to preserve full information of spiking activity 15 . Moreover, thresholding favours large (pyramidal) neurons with high amplitude spikes and other neurons very close to the electrode tips 14,15 . The small neurons (with low amplitude spikes) or slightly farther neurons may not be detected. In contrast to that, ESA integrates the contribution from all neurons in the vicinity of the recording electrode tip.
We next investigated the impact of different number of LFP channels on the inference performance of spiking activity. The inference performance of ESA was found to be consistently higher than those of SUA and MUA regardless of the number of LFP channels. Moreover, ESA reached an inference performance plateau quicker (i.e. in fewer number of LFP channels) than SUA and MUA. This may suggest that ESA exhibits larger spatial coverage (encompassing smaller and farther neurons), which in turn contains more spiking information and higher interchannel correlation 23 . We also found that the higher LFP interchannel correlation (information redundancy), the quicker the inference performance of spiking activity reaches its plateau.
To gain more insights into the relationship between LFP and spiking activity, we quantified the LFP channel importance for ESA inference using two metrics, which are average CC and average coefficient. We sorted the average LFP channel importance scores according to the interelectrode distance in an ascending order. We found that four out of six cases (across three subjects and two metrics) exhibited decreasing linear trend (negative slope) between the LFP channel importance score and interelectrode distance with three cases being statistically significant. Similar trends were also observed in the case of SUA and MUA inferences. These results are in good agreement with the previous findings 26,34 which showed that the inference performance degraded with the increase of interelectrode distance. This trend could be related to contribution of neuronal population giving rise to LFP signal which has been found to decay with the increasing electrode distance 5,7 . In addition, previous study reported that phase synchronisation between spiking activity and LFP (termed spike-field coherence) decreased with the increasing interelectrode distance 47 .
In this study, LFP was obtained using a low-pass filter with lower cut-off frequency (100 Hz) than the one typically used (300 Hz) to remove possible contamination by spike waveforms from nearby neurons 48 . Since slow modulation of spiking activity as low as 10 Hz can contribute to the LFPs 8,32,49 , this raised the question of how the cut-off frequency impacts the inference performance. We therefore performed the same procedure and compared the inference performance using different cut-off frequencies {10, 50, 100, and 300 Hz}. Results across subjects revealed that there was no significant difference in inference performance, indicating negligible contamination by nearby spike waveforms (see Supplementary Figure 6). This aligns well with previous studies which suggested that the contribution of slow modulation of spiking activity is negligible compared to other sources such as synaptic activity and membrane potential oscillations 24,50 .
Our study uses Utah array with fixed electrode depth (1 mm for dataset I and 1.5 mm for dataset II) which records neural signals from likely the same layer (presumably layer 5 51,52 ). Since neural signals have been shown to exhibit different properties 53,54 and contain different amount of task-related information 55,56 across depths or layers, it is unclear whether the same findings can be observed using electrode with recording tips on varying depths (also known as laminar or linear electrode). This can be a potentially interesting avenue for future research.
In summary, we have shown that LFPs, in the form of LMPs, carry substantial information about spiking activity, particularly ESA. Since spiking activity has been widely used as an input signal for BMI, our finding www.nature.com/scientificreports/ corroborates the increasingly accumulating evidence that LFPs can be used as an alternative input signal. This is especially relevant for chronic recordings where the spike signals have been found to be unstable or even degrading over long periods of time. In this case, LFP-based BMI can be implemented using biomimetic approach for decoding behavioural parameters 38,57,58 . Alternatively, LFP can also be used to infer spiking activity which is then applied to BMI decoding using biofeedback (operant conditioning) approach 13,30 .

Methods
Neural recording and behavioural task. Electrophysiological recordings were obtained from two public neural datasets, herein referred to as dataset I 59 and dataset II 60 . These datasets were recorded from the motor cortex area of three Rhesus macaque monkeys (Macaca mulatta) while performing predefined tasks with 96-channel silicon-based intracortical microelectrode (Utah) array. The data acquisition along with the behavioural task associated with each dataset are briefly described as follows.
Dataset I. This dataset was recorded from a male monkey (indicated by I) while performing a point-to-point task. The monkey had to to reach randomly drawn circular targets uniformly distributed around an 8-by-8 square grid. A sequence of new random targets was presented immediately and continuously after target acquisition without an inter-trial interval. The recordings were made by using a Utah array (platinum contact, 400 k nominal impedance, 400 μm interelectrode spacing, 1 mm electrode length, and 4 mm × 4 mm area) referenced to a silver wire placed under the dura (several cm away from the electrodes). The recordings were pre-amplified and filtered with a 4th-order 7.5 kHz low-pass filter with a roll-off of 24 dB per octave. The recordings were then digitised with 16-bit resolution at 24.4 kHz sampling rate, which are hereinafter called as raw neural signals. More detailed description of the experimental setup is described elsewhere 61 . For our experiments and analyses, we used a total of 26 recording sessions spanning 7.3 months between the first (I20160627_01) and last (I20170131_02) sessions, with duration ranging from 6 to 13.6 min (average of 8.88 ± 1.96 min).
Dataset II. This dataset was recorded from two monkeys (female and male) indicated as L and N, respectively, while performing an instructed delayed reach-to-grasp task. The monkey had to grasp an object with one of grip types (side grip or precision grip) and displace it with either a high or low pulling force. In each trial, the monkey had to perform one of four possible trial types (from a combination of grip type and force type) randomly drawn from an equiprobable distribution. Before initiating the movement, the monkey had to wait for 1000 ms (preparatory delay). The trial was successful if the monkey could reach, grasp, pull and hold the object for 500 ms within the position window. The recordings were made by using a Utah array (iridium oxide contact, 50 k average impedance at 1 kHz, 400 μm interelectrode spacing, 1.5 mm electrode length, and 4 mm × 4 mm area) referenced to two wires connected to the connector pedestal. The recordings were amplified and filtered with a 1st-order 0.3 Hz high-pass filter and a 3rd-order 7.5 kHz Butterworth low-pass filter. The recordings were then digitised with 16-bit resolution at 30 kHz sampling rate which are then called as raw neural signals. This dataset contains two recording sessions: L101210-001 (11.49 min) for monkey L and I140703-001 (16.43 min) for monkey N. A detailed description of the experimental setup and the corresponding behavioural task is provided elsewhere 62 . www.nature.com/scientificreports/ Signal processing and feature extraction. Signal processing steps from the raw neural signals to obtain LFP, ESA, MUA, and SUA along with their associated features are illustrated in Fig. 7a-d. We briefly describe the signal processing and feature extraction steps as follows.
Local field potential (LFP). LFP was obtained by low-pass filtering the raw neural signal with a 4th-order Butterworth filter at 100 Hz and then downsampling it to 1 kHz. A cut-off frequency of 100 Hz was selected to eliminate possible contamination from multiunit spiking activity within higher LFP band ( 100 − 300 Hz). The filtering was performed in forward and backward directions to avoid any phase shift. We extracted six different LFP features consisting of average amplitude feature called local motor potential (LMP) 40 and average spectral power features from five different frequency bands: delta (0.5-4 Hz), theta (4-8 Hz), alpha (8-12 Hz), beta (12-30 Hz) and gamma (30-100 Hz). The LMP was computed using time-domain moving average filter with 256 ms rectangular window. The spectral power feature was computed using short-time Fourier transform (STFT) with a 256 ms Hanning window and then averaged across frequency bins within each band. All the features were extracted through an overlapping fashion (206 ms overlap) to yield a sample every 50 ms.
Entire spiking activity (ESA). ESA was obtained by high-pass filtering the raw neural signal with 1st-order Butterworth filter at 300 Hz. The filtered signal was then full-wave rectified, low-pass filtered with 1st-order Butterworth filter at 12 Hz and downsampled to 1 kHz. All the filtering processes were performed in both forward and (a,c,e) Scatter plot of LFP importance score over inter-electrode distance (μm) from monkey I, L, and N, respectively. Red solid lines represent linear regression lines used to test whether or not there is a significant linear trend between inter-electrode distance and LFP channel importance score. Asterisks indicate that there is a significant linear trend (two-tailed one-sample t test; **p < 0.01 , ***p < 0.001 ). (b,d,f) Examples of heatmap of LFP channel importance score for ESA inference from monkey I (channel 5), monkey L (channel 45), and monkey N (channel 62), respectively. The importance score is mapped onto a 10 × 10 grid spatially corresponding to Utah electrode array configuration. white numbers inside the grids denote the ESA channel being inferred. White boxes on the grid represent unused (unconnected) electrodes. The larger the average CC, the more important is the channel for the inference. www.nature.com/scientificreports/ backward directions. We extracted one ESA feature using time-domain moving average filter with a rectangular window of 256 ms width and 206 ms overlap (similar to that of LMP feature).
Multiunit activity (MUA). Both datasets comprise the raw neural signals and pre-processed spikes (detected and sorted spikes). The spike waveforms were extracted by (i) filtering the raw neural signals with Butterworth filter (4th-order bandpass filter from 500 Hz to 5 kHz for monkey I; 4th-order high-pass filter at 250 Hz for monkey L; 2nd-order bandpass filter from 250 Hz to 5 kHz for monkey N), (ii) storing snippets of the filtered signals that crossed certain threshold values (48/64, 48, and 38 samples for monkey I, L, and N, respectively). MUA was defined as all the detected spikes and represented by their spike times. Details of the spike detection are described elsewhere by Makin et al. 61 for dataset I and Brochier et al. 62 for dataset II. We extracted spike rate feature by averaging the number of spikes within 256 ms window size (overlapped by 206 ms) to obtain feature sample every 50 ms. We only included MUA with spike rate exceeding 0.5 Hz for our experiments and analyses, which yielded 91, 96, and 96 units for monkey I, L, and N, respectively (see Table 1).
Single-unit activity (SUA). SUA was obtained by aligning the extracted spike waveforms, reducing the dimensionality to a few principal components, and then sorting them into distinct putative single units via certain algorithms (operator defined templates for dataset I; K-Means and Valley Seeking for dataset II). Details of the spike sorting procedure for each dataset can be found in other studies 61,62 . We computed spike rate from SUA Figure 6. LFP channel importance score (quantified in terms of average coefficient) for ESA inference across subjects. (a,c,e) Scatter plot of LFP importance score over inter-electrode distance (μm) from monkey I, L, and N, respectively. Red solid lines represent linear regression lines used to test whether or not there is a significant linear trend between inter-electrode distance and LFP channel importance score. Asterisks indicate that there is a significant linear trend (two-tailed one-sample t test; ***p < 0.001 ). (b,d,f) Examples of heatmap of LFP channel importance score for ESA inference from monkey I (channel 5), monkey L (channel 45), and monkey N (channel 62), respectively. The importance score is mapped onto a 10 × 10 grid spatially corresponding to Utah electrode array configuration. Black numbers inside the grids denote the ESA channel being inferred. White boxes on the grid represent unused (unconnected) electrodes. The larger the average coefficient, the more important is the channel for the inference. where y ik is the k-th response for the i-th observation; b 0k is the regression intercept for the k-th response; b jk is the j-th predictor's regression slope (also called coefficient) for the k-th response; x ij is the j-th predictor for the i-th observation; e ik is multivariate Gaussian residual (error) term that accounts for all other factors influencing the response variables other than the predictors. Equation (1) can be written in matrix form as The regression coefficients ( B ) were estimated using ordinary least squares (OLS), that is, by minimising the sum of squared error (i.e. difference between the observed responses and predicted responses). The solution for OLS problem is given below.
where B denotes the (p + 1) × m coefficient matrix, X represents the n × (p + 1) design matrix (LFP features), Y represents the n × m response matrix (ESA, SUA, or MUA). The MLR model was implemented using Scikitlearn 63 (v0.22.1) machine learning library in Python programming language.
Performance evaluation and metrics. Each dataset was divided into 10 non-overlapping contiguous blocks of equal size which were then categorised into three sets: training set (8 concatenated blocks), validation set (1 block) and testing set (1 block). We trained (fit) an MLR model on the training set and evaluated it on the testing set. We repeated the performance evaluation on 10 different blocks of testing set. The model's input (LFP features) and output (ESA, SUA, or MUA features) were standardised (i.e. z-transformed) to have zero mean and unit variance. The model's performance evaluation was assessed using Pearson's correlation coefficients (CC) metric. CC measures the linear correlation between the actual and inferred spiking activity. In addition, we also evaluated the model performance using another metric called root mean square error (RMSE). It is a measure of the average magnitude of the inference error. Both CC 28,30,35,64,65 and RMSE 28,35 metrics have been used in several prior studies. CC and RMSE are defined as follows: (1) www.nature.com/scientificreports/ where y t and ŷ t represents the actual and inferred spiking activity at time step t, respectively and N is the total number of observations (i.e. samples).

Impact of number of LFP channels on inference performance.
To examine the impact of number of LFP channels on the inference performance, we selected randomly p distinct LFP channels where p = {1, 5, 10, 15, . . . , 90, 95} from a total of 96 channels. We then used these p LFP channels to train an MLR model and evaluate its performance for inferring all channels or units of spiking activity (ESA, SUA, and MUA). This procedure was repeated for 30 iterations to obtain the average inference performance along with its confidence interval. We also measured LFP interchannel correlation by computing Pearson's correlation coefficient from all possible pair combination of LFP channels, yielding 96 × 96 correlation matrix.
LFP channel importance for inference performance. The LFP channel importance for ESA inference performance was measured using two approaches. First, we quantified the LFP channel importance in terms of average CC with an MLR model that was fit independently (separately) on each LFP channel ( p = 1 ). Second, we took the absolute value of coefficients (weights) of an MLR model that was fit simultaneously on all LFP channels ( p = 96 ). We then averaged the coefficients associated with each LFP channel. The larger the average CC and average coefficients, the more important is that channel for ESA inference. We calculated the distance between the electrode location of ESA channels (output) and the electrode location of LFP channels (input). The LFP channel importance was then sorted based on the interelectrode distance in an ascending order. The interelectrode distance ranged from 0 to 4561 μm. We used the slope of linear regression to examine whether there is a significant linear trend between the interelectrode distance and LFP channel importance score (twotailed one-sample t test).

Statistical analysis.
For each session, the mean and standard error of the mean (SEM) of the inference performance were evaluated across the number of units of spiking activity (ESA, SUA, and MUA) on 10 different blocks within the testing set. To test for significant effects between a pair of different LFP features or different spiking activity, we used a paired two-tailed t-test whenever the difference between the pairs follows normal distribution, otherwise we used Wilcoxon signed-rank test. The significance level ( α ) was set to 0.05. When the results are visualised using boxplot, the middle horizontal line and circle mark inside each box represent the median and mean, respectively. The coloured solid box represents interquartile range (IQR) from 25th to 75th percentiles. The whisker extends 1.5 times the IQR. All the analyses were conducted in Python (v3.6.10).