Abstract
Microstimulation can modulate the activity of individual neurons to affect behavior, but the effects of stimulation on neuronal spiking are complex and remain poorly understood. This is especially challenging in the human brain where the response properties of individual neurons are sparse and heterogeneous. Here we use microelectrode arrays in the human anterior temporal lobe in 6 participants (3 female) to examine the spiking responses of individual neurons to microstimulation delivered through multiple distinct stimulation sites. We demonstrate that individual neurons can be driven with excitation or inhibition using different stimulation sites, which suggests an approach for providing direct control of spiking activity at the single-neuron level. Spiking responses are inhibitory in neurons that are close to the site of stimulation, while excitatory responses are more spatially distributed. Together, our data demonstrate that spiking responses of individual neurons can be reliably identified and manipulated in the human cortex.
SIGNIFICANCE STATEMENT One of the major limitations in our ability to interface directly with the human brain is that the effects of stimulation on the activity of individual neurons remain poorly understood. This study examines the spiking responses of neurons in the human temporal cortex in response to pulses of microstimulation. This study finds that individual neurons can either be excited or inhibited depending on the site of stimulation. These data suggest an approach for modulating the spiking activity of individual neurons in the human brain.
Introduction
Clinicians and researchers have been using direct electrical stimulation of the human brain for decades to treat neurologic disease and to advance our understanding of human brain function (Borchers et al., 2011; Lozano et al., 2019). However, the precise effects of electrical stimulation on neural activity, particularly at the level of individual neurons, still remain relatively unknown. The use of smaller stimulation currents delivered through microelectrodes has the potential of providing more targeted stimulation that can modulate the activity of individual neurons (Butovas and Schwarz, 2003; M. Histed et al., 2009). Microstimulation could therefore, in principle, afford both clinicians and researchers selective control over neuronal firing, which could in turn enable direct investigations into the causal role neuronal activity may play in human behavior and our ability to interface directly with the human brain to address neurologic disorders (Lebedev and Nicolelis, 2006; Schiller and Tehovnik, 2008; Lewis et al., 2015).
The effects of stimulation on neuronal spiking responses, however, are complex and may be related to the location of the stimulating electrode with respect to neurons and axons of passage, which in turn are dictated by the underlying neural architecture of the cortical regions being investigated. This complexity highlights a central challenge for using targeted microstimulation to manipulate neuronal activity in humans, which is that in many brain regions there is considerable sparseness and heterogeneity in the response properties of neighboring neurons, even within the same patch of cortex (Jang et al., 2017; Wittig et al., 2018; Vaz et al., 2020). High-density recordings have been used over the past decade to acquire and analyze larger populations of spiking neurons in animals, which has advanced our understanding of how these populations of neurons contribute to perception (Zhang et al., 2011; Pagan et al., 2013), cognition (Kaufman et al., 2015; Peixoto et al., 2021), and action (Kao et al., 2015; Willett et al., 2021). Despite these advances, though, the ability to use microstimulation to manipulate neuronal activity within these sparse and heterogeneous populations codes has not scaled in proportion (Carrillo-Reid et al., 2017).
To address this challenge, therefore, it is necessary to first characterize the physiological effects of microstimulation at the level of individual neurons. Here we characterize the effects of microstimulation from multiple locations across a patch of human temporal lobe, an association cortex characterized by neuronal heterogeneity. Using a semichronic microelectrode array (MEA) custom-built to allow microstimulation and recording across distributed sites, we capture the spiking responses of individual neurons in 6 participants as we deliver biphasic pulses of electrical stimulation through multiple separate microelectrode contacts. We show, for the first time, in the human brain physiological effects of stimulation in the same neurons from multiple distinct stimulation sites. We demonstrate that a single neuron can be driven with excitation or inhibition using different stimulation sites, which suggests a possible approach toward achieving selective control at the single-neuron level. The responses to microstimulation are dominated by inhibition in regions that are spatially proximal to the site of stimulation, but we also find a smaller yet significant number of excitatory responses that are more spatially distributed. In both cases, the responses to stimulation are modulated by stimulation amplitude. Together, our data demonstrate that neuronal spiking responses to microstimulation in the human brain are complex yet can be reliably observed and characterized. Our data therefore provide a step toward future approaches that may seek to manipulate the activity of individual neurons to test causal contributions to a population code, and ultimately toward establishing precise communication of information with the human brain.
Materials and Methods
Participants
Six participants (34.6 ± 3.75 [mean ± SEM] years; range 24-46 years; 3 female) with drug-resistant epilepsy underwent a surgical procedure in which platinum recording contacts were implanted subdurally on the cortical surface as well as deep within the brain parenchyma. In each case, the clinical team determined the placement of the contacts to localize epileptogenic regions. In all the participants investigated here, the clinical region of investigation was the temporal lobes.
For research purposes, we additionally placed two 64-channel MEAs (3.2 × 3.2 mm, Cereplex SI; Blackrock Microsystems) ∼2 cm apart in the middle temporal gyrus of the anterior temporal lobe and ∼2 cm from the temporal pole in addition to the subdural contacts (see Fig. 1A). We implanted MEAs only in participants with a presurgical evaluation, indicating clear seizure localization in the temporal lobe, and the implant site in the ATL was chosen to fall within the expected resection area. Each MEA was placed in an area of cortex that appeared normal both on the preoperative MRI and on visual inspection. Four of the 6 participants received a surgical resection, which included the tissue where the MEAs were implanted. One participant received a surgical resection that involved regions posterior to the implanted MEA in the temporal lobe and had a considerable improvement in seizures following surgery. The other participant did not have seizure activity that was captured during the monitoring period, which suggested a safe and effective surgical resection would be possible, and did not experience any change in seizure frequency following removal of the monitoring electrodes. Neither participant experienced any noted change in cognitive function.
Data were collected at the Clinical Center at the National Institutes of Health. The Institutional Review Board approved the research protocol (11-N-0051), and informed consent was obtained from the participants and their guardians explicitly for the placement of the MEAs. All analyses were performed using custom-built MATLAB code. Data are reported as mean ± SEM unless otherwise specified.
Microelectrode recordings
We captured spiking activity and local field potential signals from the MEAs. Each MEA contains an 8 × 8 grid of microelectrodes, with each electrode spaced 400 μm apart and extending 1.5 mm into the cortical surface. Postoperative paraffin blocks of the resected tissue demonstrated that the electrodes extended approximately halfway into the 3-mm-thick gray matter. We digitally recorded signals at 30 kHz from all 128 microelectrodes in each participant using the Cereplex I and a Cerebus data acquisition system (Blackrock Microsystems) with 16-bit precision and a range of ±8 mV. We recorded signals while applying microstimulation to combinations of up to 15 preselected electrodes (see Experimental design). Digitized neural data, external power, and stimulation current were communicated via two silicon-encased wires that exited the scalp in a bundle with other clinical wires far from the primary incision site to minimize the risk of infection.
To extract unit spiking activity, we rereferenced each microelectrode's signal offline by subtracting the mean signal of all the electrodes in the MEA, and then used a second-order Butterworth filter with zero phase distortion to bandpass the signal between 0.3 and 3 kHz. Because microstimulation can create large noise artifacts in the neural data with decays that typically returned to baseline within 5 ms, we zeroed out all neural data from 5 ms before to 5 ms after each stimulation event before sorting the data for single-unit activity. Any residual effects of stimulation following artifact removal did not impact the spike-sorting process, which relies on extracting discrete large deflections in the filtered trace and then projecting potential spike waveforms into a lower dimensional space. Before spike sorting, raw data were filtered between 300 and 3000 Hz, which also mitigates the effects of the stimulation artifact. Using a spike sorting software package (Plexon Offline Sorter), we identified spike waveforms by manually setting a negative or positive voltage threshold depending on the direction of putative action potentials. The voltage threshold was set to include noise signals used in calculating unit isolation quality (see below). Waveforms (duration, 1.067 ms; 32 samples per waveform) that crossed the voltage threshold were stored for spike sorting. Spike clusters were manually identified by viewing the first two principal components, and the difference in peak-to-trough voltage (voltage vs time) of the waveforms. We manually drew a boundary around clusters of waveforms that were differentiable from noise throughout the experimental session. Across all experimental sessions in all participants, we identified a total of 717 putative single units across all sessions (average of 119 ± 34 units per participant). The average baseline spike rate across all units was 0.78 ± 0.55 Hz (median 0.28 Hz). In order to reliably detect stimulation effects, we only retained units with a minimum firing rate of 0.5 Hz, leaving a total of 258 units (43 ± 11 units per participant) and 3870 unique combinations of units and stimulation sites.
Single-unit recording quality measures
Because of variability in the signal quality across recordings and the subjective nature of spike sorting, we quantified the quality of each unit by calculating an isolation score and signal-to-noise ratio (SNR) (Joshua et al., 2007). The isolation score quantifies the distance between the spike and noise clusters in a 32-dimensional space, where each dimension corresponds to a sample in the spike waveform. The spike cluster consisted of all waveforms that were classified as belonging to that unit, and the noise cluster consisted of all waveforms that crossed the threshold that were not classified as belonging to any unit. The isolation score is normalized to be between 0 and 1, and serves as a measure to compare the isolation quality of all units across all experimental sessions and participants. Across participants, the mean isolation score for all retained units was 0.84 ± 0.17 (median 0.91).
In addition to isolation quality, we computed the SNR for each unit as follows:
Experimental design
We used an external microstimulator (CereStim96, Blackrock Microsystems) to apply microstimulation through combinations of up to 15 preselected microelectrodes (see Fig. 1A, inset). The Cerestim96 connects to an externalized pigtail lead using a custom cable, and the system is electrically isolated with a floating ground tied to the participant. During each stimulation experiment, we captured the continuously digitized neural data as well as a digital reference signal marking the precise time of stimulation onset and offset output by the external stimulator. For all stimulation, we applied monopolar biphasic stimulation through each microelectrode such that the metal Cereplex SI digitizer case acts as the cathodic current sink. Stimulation currents are hardware limited to 100 μA. This current level was determined to be safe for long-term microstimulation through these microelectrodes based on current density through the electrode tips, extensive testing by Blackrock, and previous studies (Merrill et al., 2005).
An epileptologist was present for all stimulation experiments to monitor the intracranial EEG recordings during stimulation. Before each experimental session, we cleared the microstimulation sites by simultaneously stimulating all 15 microelectrodes using a train of stimulation pulses at 100 Hz for 200 ms (21 pulses) at the lowest stimulation current of 5 μA per microelectrode. If no afterdischarges or clinical responses were noted, we repeated stimulation after incrementally increasing the stimulation current. We continued this process until we reached the maximum current level of 100 μA per microelectrode. The epileptologist never observed any afterdischarges or clinical responses during the microstimulation experiments, even in cases when all 15 microelectrodes were stimulated with the maximum current simultaneously (total current = 100 μA × 15 microelectrodes = 1.5 mA).
We designed a custom-built MATLAB graphical user interface to run the microstimulation experiments while providing an easy way to pause or stop stimulation based on clinical considerations. Our graphical user interface allows for control of stimulation amplitude, pulse train (frequency and length), and microelectrodes stimulated. During each experimental session, we applied stimulation to a single microelectrode at a time during each trial (see Fig. 1B). Hence, we consider each stimulation event at each microelectrode as a single trial. The custom graphical user interface sequentially cycled through each of the 15 microstimulation sites, with an intertrial interval of at least 1 s between successive sites (intertrial interval of 1.5 s in 1 participant, 1.25 s, 1s in the remaining participants). Because we serially delivered microstimulation to each of the 15 stimulation sites in order before returning back to the first site, the interval between two separate stimulation trials in the same microelectrode was at least 15 s. We sequentially stimulated each of 15 stimulation sites in this manner until we collected a minimum of 20 trials per site in each experiment. Each experimental session lasted ∼30 min. We performed the experiments during the monitoring period between day 2 and day 10 after implantation. Participants were awake and in their monitoring room and were instructed to rest with their eyes open during each experimental session. No participant reported any sensation during the experimental sessions, or gave any indication that they knew when stimulation was delivered. Each participant completed one experimental session, and 2 of the participants completed a second session.
For each trial at each individual microelectrode, we provided microstimulation using charge balanced biphasic pulses (300 μs up, 55 μs flat, 300 μs down). We used a maximum current amplitude of 100 μA per phase (200 μA peak-to-peak) based on the hardware and safety limits, although in a subset of participants we also examined the responses to stimulation for a range of amplitudes between 5 and 100 μA. In these participants, we completed stimulation at one amplitude across sites, gave the participants a short 1 min rest period before repeating stimulation across sites at the next amplitude. On each trial for all experimental sessions, we delivered stimulation pulses at 100 Hz (50 Hz in one participant, P2). Depending on the experimental session, on each trial, we delivered stimulation for a total duration of one pulse, 100 ms (11 pulses), or 200 ms (21 pulses).
Statistical analysis
During each trial, we captured single-unit spiking activity from each unit and estimated the instantaneous spike rate by calculating the total number of spikes in sliding 100 ms windows (10 ms step sizes, 90% overlap). We did not include the first 10 ms after the end of stimulation in each trial for analysis because of the stimulation artifact. We designated the last 400 ms of neural data in each trial, by which point spiking rates had often returned to prestimulation levels, as the baseline. We used the mean and SD of the spike rate across all of the baseline periods to generate a z-scored spike rate for each unit. For each unit, because we use this period from all trials, the pooled baseline activity reflects the baseline activity present on average even when stimulation is applied to the different microelectrodes. Because we sequentially and repeatedly delivered stimulation to each of the stimulation sites, we are able to separately assess the average response of any individual unit to stimulation at a single site across multiple trials.
In order to reliably identify and classify response types, we used a nonparametric test on each 100 ms sliding window in the spike raster to flag statistically significant changes in z-scored spike rate across trials relative to that units baseline rate (two-sample rank sum test). This method flagged a large number of units and stimulation sites with potential stimulation effects at specific time points. However, on close examination of each raster, we found that using a single threshold for this statistical test missed several inhibitory responses (i.e., no flags) in units with a low baseline firing rate that dropped to zero for several hundred milliseconds following stimulation, and conversely several apparent false alarms where isolated time points were flagged as excitatory or inhibitory several hundred milliseconds after stimulation with no obvious corollary in the raster. We systematically varied the α value of the test and the minimum number of consecutive significant time windows required for flagging but found no single combination of these two parameters that reliably flagged all three response types with a low false alarm rate.
We therefore used four separate combinations of α values and number of consecutive windows that reliably identified each of the clearly observable effects in the rasters with few false alarms. First, we considered a unit as exhibiting an excitatory response to a stimulation site if a single window within the first 250 ms following stimulation exhibited an increase in spike rate that was flagged with the nonparametric rank sum test with an α threshold of 0.000667 (α of 0.01 divided by 15 time windows). This test for single windows identifies strong, brief changes in spike rate. Second, we identified strong and brief inhibitory responses using a similar approach. In this case, however, the decreases in spike rate are limited by the floor (a spike rate of zero). We therefore used an α threshold of 0.00333 for the nonparametric test on each window (α of 0.05 divided by 15 time windows, which is a threshold that is 5 times as liberal as the threshold for excitatory responses). We also considered cases when the responses to stimulation were more sustained but less strong. In these cases, we identified an inhibitory response if 10 consecutive sliding windows each exhibited a decrease in spike rate with an α that exceeded a threshold of 0.1. We only considered these cases if all 10 windows exceeded this threshold. Similarly, we identified excitatory responses when 10 consecutive windows each exhibited an increase in spike rate with an α that exceeded a threshold of 0.02 (maintaining the same ratio of being 5 times more conservative for the excitatory responses). For both excitatory and inhibitory responses, if 10 or more consecutive windows satisfied these thresholds, all of the windows were considered as part of a significant deviation from baseline spike rate. We considered any unit that satisfied one or more of these four requirements as exhibiting a putative response, and the response type was assigned accordingly, with units satisfying both excitatory and inhibitory requirements being labeled as complex (see Figs. 1C, 2).
Two human judges (DY and JW) independently inspected the rasters from all spiking responses that were flagged using this method to verify that rasters with visually obvious response types were being properly identified. In an attempt to see whether a simpler set of statistical requirements might yield comparable classification of response types, we compared several algorithms that just used one or two α values, and did not explicitly account for brief versus sustained effects. We compared a lenient threshold (one or more time bins surpassing an α value of 0.05), a moderate threshold (lenient, with correction for multiple comparisons in time using false discovery rate), a strict threshold (lenient, with Bonferroni-corrected α of 0.00022, which is an α of 0.05 divided by 15 time points and 15 stimulation sites), and a dual threshold that yielded the closest match to our four threshold technique (one or more time bins surpassing an α value of 0.01 for excitatory effects, 0.20 for inhibitory effects; both with FDR correction for multiple comparisons in time). The three methods with single α values either had a large number of excitatory false alarms or a low number of missed inhibitory responses. The dual threshold method, which accounted for different sensitivities to increases versus decreases, was the closest match to our method, although a few individual example responses were clearly mislabeled. Therefore, we continued to classify all significant responses using the four separate combinations of α values and number of consecutive windows.
To verify that stimulation can generate different spiking response types, we used a clustering approach. For every unique combination of unit and stimulation site, we extracted the z-scored time series of spiking response. Each response can be considered as a point in n-dimensional space, where n is the number of 100 ms sliding windows used to capture the spiking response in each trial. We computed the distance between every pair of responses in this high dimensional space, and then used multidimensional scaling to project these distances into the two lowest dimensions. Each point in this lower dimensional space therefore represents the response of a single unique combination of unit and stimulation site. We visualized where our identified responses fall in this lower dimensional space (Fig. 1G). In this visualization, each response type appears to cluster to a different region of the low-dimensional space. We separately used a k-means clustering procedure to divide the responses that were projected in the lower dimensional space into three clusters in an unsupervised manner (Fig. 1H). The manual clustering and the unsupervised clustering appear similar. To statistically test whether the two clusterings are more similar than would be expected by chance, we assigned labels to each identified response type. In the manual cluster, we assigned each response type a value of 1, 0, or −1 depending on whether the response was manually classified as excitatory, no response, or inhibitory, respectively. We labeled complex responses as 2 for this comparison, although the unsupervised clusters were divided into only three groups. Hence, the collection of all unique combinations of unit and stimulation site constitutes a vector of assigned values. We then assigned labels to the response types that were clustered using k-means clustering to generate a vector of assigned values. The response of every unique combination of unit and stimulation site occupies the same index in both vectors. We computed the correlation between the two vectors to get a true measure of how similar the clustering is between our manual approach and the unsupervised approach. We then shuffled the response labels 1000 times to generate a distribution of permuted correlations, and compared the true correlation to this distribution to generate an empiric p value of how statistically likely we would observe the true correlation by chance.
For each responsive unit, we characterized the dynamics of the inhibitory, excitatory, and complex responses. For each unique unit and stimulation site combination, we computed the trial-average spike rates using a 100 ms sliding window with 10 ms step size (90% overlap). We then fit the average response with a decaying exponential of the form: y = (y1 − y0) × (1 − exp−(t − toffset)/tau) + y0, where y is the observed spike rate, t is the center of the sliding window, y0 is set to the baseline spike rate, toffset is the starting point of the fit (defined below), and y1 and tau are free parameters of the fit corresponding to the peak spike rate and time constant. We set the starting point of each fit, toffset, to the last time bin containing the maximum spike rate for excitatory responses, and the minimum spike rate for inhibitory responses. For most excitatory responses, the maximum spike rate was observed at a single time point. However, it was common for inhibitory responses to have sustained periods with zero spikes during the inhibitory period; and in those cases, we set toffset as the last bin to have zero spikes.
To examine the spatial relationship between response likelihood, response magnitude, and the radial distance between stimulation and recording sites, we first identified the spatial location of recorded units and their responses to stimulation sites at each spatial location. For all units and stimulation sites on the same array, we calculated radial distances in units of interelectrode distances (400 μm). In some cases, we observed responses on one array following stimulation on the other array (see Fig. 5). In those cases, we assigned the distance between stimulation site and the unit response as 8 mm. After pooling responses from all participants, we calculated the average percentage of responding units at each radial distance. We fit these data with a 1D Gaussian of the form: y = (y1 − y0) × exp−x/σ2 + y0), where y is the observed percentage of responding units, x is the radial distance in units of microelectrodes, y0 is a constant set to the percentage of responses observed when stimulation is on a different array, and y1 and σ are fit parameters corresponding to the peak and spread of the Gaussian fit. We deemed a fit significant if σ was significantly different from zero (95% CIs). For excitatory responses, the fit was not significant, and we report the mean value across all radial distances. We used a similar approach to characterize the magnitude of responses with respect to radial distance. In this case, we computed fits across the values from each response (rather than averaging the proportion of responses at each location). A small proportion of responses were identified, even when the stimulation site was on a different array that was >2 cm away (see Fig. 5), and we did not include these responses in our fits.
In 3 participants, we tested six different stimulation amplitudes, spanning 5-100 μA, to determine whether individual units showed systematic changes in response amplitude with stimulation amplitude. We refer to this as a rate code of stimulation amplitude. Across these participants, we identified 184 unique combinations of units and stimulation sites with an identified response to the maximum stimulation current. For each of these responsive combinations, we defined the time window of rate coding as the period of significant inhibitory or excitatory responses to the maximum stimulation current, and then fit a neurometric function to the change in spike rate across stimulation amplitudes. We used a sigmoidal neurometric function analogous to the psychometric function used in previous microstimulation studies: y = y0 + ((y1 − y0)/(1 + exp(x − x0)/dx)), where y is the observed neural response, x is the stimulation amplitude, y0 and y1 are constants set to the mean observed response to minimum and maximum stimulation levels, respectively, and x0 and dx are free parameters of the model indicating the amplitude of maximum sensitivity and scale of that sensitivity, respectively. Therefore, each fit resulted in an x0 value describing the amplitude of maximal change in spike rate and dx, which describes the width of the maximal change in spike rate over stimulation amplitude. We considered a fit valid if the 95% CI for dx did not include zero, and the 95% CI for x0 did not span the entire rate of stimulation amplitudes (i.e., did not include both 5 and 100 μA). We then computed an optimized rate code by iteratively searching for the ideal time window in which to analyze rates for each combination of unit response and stimulation site. The percentage of valid neurometric fits was similar using both methods (standard rate and optimized rate), as was the stimulation amplitude causing a maximal change in neurons response (x0). However, on average, the fit quality using an optimized rate code was significantly better than that of the optimized time code or standard rate code, and so we present the results of the optimized fits.
Data and code availability
The data and code that support the findings of this study are available from the corresponding author on request and are also available for public download at https://research.ninds.nih.gov/zaghloul-lab/downloads.
Results
We examined neuronal spiking in the human temporal lobe in response to microstimulation in 6 participants (34.6 ± 3.75 [mean ± SEM] years; range 24-46 years; 3 female) who were being monitored for seizures after placement of intracranial electrodes. In addition to the clinical electrodes, we placed two 8 × 8 MEAs in the anterior temporal lobe in each of these participants. We custom-built the arrays such that a total of 15 microelectrodes distributed across fixed locations in both arrays could be used to deliver microstimulation while we simultaneously recorded neuronal spiking activity from single units captured across the remaining microelectrodes (Fig. 1A; see Materials and Methods). Across participants, we recorded spiking activity from 119 ± 34 units per participant (43 ± 11 that met 0.5 Hz spiking criteria; range 19-88) and focused our analyses on the responses of these units to microstimulation.
Microstimulation evokes selective responses in neural firing. A, Schematic illustrating implant of two MEAs in the middle temporal gyrus of the left anterior temporal lobe. Each array contains an 8 × 8 grid of microelectrodes. Seven or 8 microelectrodes dispersed across each array can be used to deliver microstimulation. B, Microstimulation occurs approximately once every second while sequentially cycling through each stimulation site (top). In this example, stimulation amplitude is 100 μA. Stimulation artifacts and spiking activity are visible in the microelectrode recordings, although the size of the stimulation artifact in the example is limited by the scale bar (bottom). C, Example spiking responses of individual units to stimulation. In each example raster, trial responses of a unit to stimulation site that evoked an inhibitory, excitatory, or complex response (effective) are shown on top, and responses of the same unit to a different stimulation site that did not evoke a response (ineffective) are shown on bottom (inset location of each unit and each stimulation site in the array). Bars represent time windows with significant increases (red) or decreases (blue) in spike rate compared with baseline (see Materials and Methods). The timing of spiking responses is variable from trial to trial. For the excitatory response shown here (middle), the spiking response exhibits a temporal jitter of 29.4 ± 12.3 ms. D, Proportion of unique combinations of stimulation site and recorded units demonstrating a response to microstimulation in each participant. E, Distribution across all recorded units of the number of stimulation sites evoking a response in each unit. F, Among the unique combinations of recorded units and stimulation site demonstrating responses, the proportion that are inhibitory, excitatory, or complex across all participants (average 64% inhibitory, 31% excitatory, and 5.2% complex). G, Projection of the stimulation responses into a lower dimensional space based on the time series of the spiking response. The manually identified response types separately cluster into groups of exhibitory, inhibitory, and complex responses. H, Unsupervised clustering of the same spiking responses demonstrates a similar clustering of responses into three groups.
In each trial during the experimental session, we applied biphasic current waveforms through a single stimulating microelectrode. Over the course of the entire session, we sequentially and repeatedly cycled through each of the 15 sites with an interstimulation interval of at least 1 s (Fig. 1B; see Materials and Methods). Hence, for each stimulation site, we were able to quickly record multiple trials in which we captured the spiking response of all identified units to microstimulation. In individual units, we observed that stimulation at one stimulation site resulted in a transient decrease in spike rate following stimulation, often apparent as a complete absence of spikes for a few hundred milliseconds (inhibitory response). In those same units, stimulation from a nearby site caused no apparent change in spike rate (Fig. 1C; see Materials and Methods). Similarly, in other individual units, we observed that stimulation through one site resulted in a transient increase in spiking activity (excitatory response) while stimulation at a different site caused no change in spike rate. Less frequently, we found units that exhibited a combination of a transient increase and decrease that was consistent across stimulation trials (complex response).
We examined the responses of all units from all participants to each stimulation site and classified each response as putatively excitatory, inhibitory, or complex (see Materials and Methods). Individual units that responded to multiple sites could show inhibitory responses to one site and excitatory or complex to another (examples in Fig. 2), indicating that the response type was not a property of the recorded unit, but instead of the unique combination of unit and stimulation site. In total, we analyzed 3870 unique combinations of recorded units and stimulation sites (258 total units above our baseline spike rate criteria of x 15 stimulation sites each). While the vast majority of unique unit and stimulation site combinations showed no effect of stimulation, 13.6 ± 1.2% of all possible combinations of units and stimulation sites across participants exhibited a putative and identifiable response (Fig. 1D). On average across all units, each unit exhibited a response to 1.9 ± 0.16 of the 15 possible stimulation sites (range 0-15; Fig. 1D). Of these combinations exhibiting identified responses to stimulation, 58.1 ± 9.1% exhibited an inhibitory response, 38.1 ± 8.8% exhibited an excitatory response, and 3.8 ± 1.1% exhibited a complex response across participants (Fig. 1F).
Example units from 5 different participants. Each unit responded in at least two different ways (excitatory, inhibitory, or complex) to different stimulation sites. Each unit is shown with a nonresponsive stimulation site, and responsive stimulation sites of multiple types of response. A, Grids represent the location of the unit (white) relative to the stimulation sites: ineffective (black), excitatory (red), inhibitory (blue), and complex (magenta). B, Spike rasters aligned to stimulation onset. Continuous spike rates (z-scored) below, with periods of significant deviations based on our response criteria indicated with solid color-coded lines.
Our approach for classifying responses provides evidence that there are indeed significant spiking responses to stimulation, but it is possible that the different response types we identified may have artificially emerged because of the different choice of thresholds used for response identification. Based on the time series of spiking responses, we therefore projected response of each unique combination of unit and stimulation site into a lower dimensional space. We visually observed that the responses we had labeled as excitatory, inhibitory, and complex appeared to cluster to different regions (Fig. 1G). We then discarded the labels we had manually assigned and used k-means clustering to group the responses into different categories in an unsupervised manner. We found a similar clustering of responses (Fig. 1H). We compared the clusters we had manually identified to those that emerged through the unsupervised approach and found that the two different approaches for grouping the responses were significantly more similar to one another than expected by chance (ρ = 0.46, p < 0.001, permutation procedure; see Materials and Methods).
We noticed substantial variation in the time course of stimulation responses within and across response types. Inhibitory responses, for example, appear to involve a silent period during which no spiking occurs, followed by a slow rise back to baseline spiking levels rather than an immediate return (Fig. 3A). We fit this return to baseline using an exponential (time constant, tau) once the silent period ended (toffset). Similarly, excitatory responses involve a large increase in spiking activity followed by a decay back to baseline that was also well fit by an exponential from the peak spike rate (Fig. 3B). Complex responses have a separate peak and dip in their spike rate, and thus recovered to baseline twice with two separate exponential fits (Fig. 3C). Inhibitory responses on average have smaller peak magnitudes, a delayed peak, and slower time course of recovery to baseline compared with excitatory and complex responses (Fig. 3D).
Time course and magnitude of stimulation responses. A, Two example inhibitory responses shown as a histogram of trial-average spike rates (100 ms windows with 10 ms steps) before and after stimulation (broken x axis). Spike rates decrease and then return to baseline, fit with an exponential time constant (τ) following an inhibitory period (toffset). Inset, Location of the units (black) and stimulation sites (blue) on the array. B, Two example excitatory responses. Spike rates exhibit an initial increase and then return back to baseline, fit with an exponential decay function. C, Two example complex responses, each fit with both an excitatory decay and an inhibitory decay. D, Average time course of z-scored spike rate across all excitatory, inhibitory, and complex responses to stimulation. Shaded error bars indicate 2 SDs. E, Distribution of z-scored peak spike rates for excitatory and inhibitory responses. F, Distribution of time constants for excitatory and inhibitory responses. G, The distribution of offset times for excitatory and inhibitory responses. E-G, For each parameter, the distribution for excitatory and inhibitory responses was significantly different (KS test; for p values, see inset).
We directly compared the peak spiking rate, time constant, and offset time from these fits across the different response types (Figs. 3E–G, 4). Consistent with observations from the average time course, excitatory responses have higher peak spike rates, whereas inhibitory responses recover more slowly and have later offset times than excitatory responses (Kolmogorov–Smirnov nonparametric test, p < 0.01, for exact p values, see insets). These data suggest that, beyond the differences in polarity between response types, the different response types also exhibit different dynamics.
We also examined the dynamics of the less frequently observed complex responses. Complex responses exhibit faster recovery than either excitatory or inhibitory responses (p = 0.003, p < 0.001, respectively, KS test; Fig. 4). Complex responses also have significantly shorter excitatory offset times and significantly longer inhibitory offset times compared with excitatory and inhibitory responses (p = 0.013, p = 0.010, respectively, KS test). Despite these differences, however, the parameters for the complex responses are still within the distribution of the other response types, suggesting that the complex responses may simply reflect a summation of fast excitation and slow inhibition pathways driven by a single stimulation site as opposed to a distinct mechanism of stimulation.
Complex response time course and magnitude. A, Distributions of fit constants (peak rate, offset, and time constants) for inhibitory and excitatory responses. B, Distributions of fit constants (peak rate, offset, and time constants) for complex response types. A KS test was used to determine whether the complex distributions were significantly different from the analogous noncomplex distributions (for p values associated with each test, see insets).
For every example unit presented thus far, we included a schematic relating the location of the unit and stimulation sites that either effectively drove a response or were ineffective (see Figs. 1, 2, insets). In each array, we counted the total number of recorded units from each individual electrode on a MEA and the proportion that responded to stimulation when delivered at each site. In individual examples, we noticed that the highest percentage of responding units was close to the stimulating electrode whether stimulation was delivered through an electrode on one corner of the array or the opposite corner of the array (Fig. 5A). This suggests that response likelihood of a unit may be related to its proximity to the electrode delivering stimulation. To explore this spatial relation systematically across participants, we computed the radial distance between each stimulation site and each recorded unit (Fig. 5B,C). Less frequently, we observed a significant response in one array following stimulation in the other array (24.45 ± 3.03% of all significant responses; 20.88% of significant inhibitory responses and 43.08% of significant excitatory responses). However, in the majority of cases, the responses were local. Across all arrays and all recorded units, the likelihood of observing an inhibitory response is ∼35% for units that are located within one electrode contact (400 μm) of the stimulation site, and falls to a baseline level of ∼4% (Gaussian fit, space constant of 4.4 interelectrode lengths, corresponding to 1.7 mm; Fig. 5D; see Materials and Methods). In contrast, the likelihood of observing an excitatory response is only between 5% and 10% across all recorded units, and appears independent of the radial distance from the stimulation site. The likelihood of observing a complex response also decreases with distance from the stimulation site, but the overall frequency of these responses was lower than the excitatory or inhibitory responses.
Spatial relation between stimulation and responses. A, Left, Number of isolated units per recording electrode in one example array. Dashed locations indicate stimulation electrodes. Right panels, The percentage of units responding to stimulation at the top right and bottom left corners of the array (yellow squares). B, Aggregate unit counts across all participants. For stimulation sites and units that were on the same array (left), the stimulation site was considered the center of a larger array and unit locations were tallied in this translated space. In order to smooth the 2D image (middle), we averaged 8 rotated copies of the images on the left (rotational step of 90 degrees, repeated after flipping the image across the x axis). For stimulation sites and units on different arrays (right), unit locations were tallied based on the absolute location of the unit on the array. C, Proportion of the aggregate unit counts showing each response type. Responses are shown to stimulation sites on the same array (middle) and when stimulation and recording sites are on different arrays (right). D, Percentage of recorded units across all participants showing each response type as a function of radial distance from the stimulating electrode. A Gaussian (red line) was fit to each distribution to estimate the spatial spread, σ, and the peak value closest to stimulation sites. Dashed red lines indicate 95% prediction intervals. Inhibitory and complex responses were well fit by the Gaussian. Excitatory responses are summarized by a single mean value that does not vary with radial distance. E, Magnitude of stimulation response as a function of radial distance from the stimulation site. Inhibitory magnitude is quantified using offset time. Excitatory and complex response magnitudes are quantified using peak spike rate.
In addition to quantifying how likelihoods of each response type are related to the distance of stimulation, we examined whether the response magnitude of engaged units varied systematically with radial distance from the stimulation site. For excitatory and complex responses, we used the z-scored peak spike rate as a measure of response magnitude. For inhibitory responses, because spiking activity drops to zero following stimulation, we instead used the time required for spiking activity of that unit to return to baseline levels (offset time, toffset). For all response types, the magnitude of the response decreases with distance from stimulation site (space constant of 2.65, 1.97, and 7.00 interelectrode lengths for inhibitory, excitatory, and complex responses, respectively; Fig. 5E).
Our data demonstrate that microstimulation can evoke both excitatory and inhibitory responses in individual units. An important question to address is whether stimulation amplitude can be used to systematically modulate the magnitude of the response. Previous studies examining microstimulation in primate somatosensory cortex have demonstrated that sensory percepts can vary with increasing stimulation amplitude, suggesting a relation between stimulation amplitude and neural response (Bartlett et al., 2005; Murphey and Maunsell, 2008; Fridman et al., 2010; Kim et al., 2015). In 3 participants, we examined this question with extended microsimulation data collection that included five amplitudes that were lower than our standard (maximum) stimulation level of 100 μA. In individual units, both inhibitory and excitatory responses appear modulated by the stimulation amplitude (Fig. 6A,C). For each unit, we used response to 100 μA to identify the time window of significant changes in spike rate (i.e., our method of identifying responsive units), and then fit a neurometric curve to the changes in spike rate during this time window across stimulation amplitudes (Fig. 6B,D; see Materials and Methods); ∼30% of the responsive combinations of units and stimulation sites showed systematic changes in spike rate that were well fit by the sigmoidal neurometric function.
Stimulation amplitude affects the magnitude and timing of spiking unit responses. A, Rasters for a single inhibitory unit responding to the same stimulation site at multiple current levels. B, Best-fit neurometric function relating input current to spike rate for this unit. C, D, Rasters and best-fit neurometric function for a single example excitatory unit. E, Neurometric functions relating stimulation amplitude to spike rate calculated across all inhibitory (top) and excitatory (bottom) responses. F, Proportion of unique combinations of recorded units and stimulation sites exhibiting responses as a function of radial distance (between unit and stimulation site) and stimulation amplitude. Distances are binned into four nonoverlapping groups to yield approximately equal numbers of units in each. G, Magnitude of responses as a function of radial distance (between unit and stimulation site) and stimulation amplitude.
To estimate population-level effects, we first pooled spike rates across all inhibitory and excitatory responses and fit a population-level neurometric function to each group (Fig. 6E). Both inhibitory and excitatory responses exhibit their maximal change in spike rate response around 72-79 μA and over a range of ∼17-24 μA (fit parameters were statistically indistinguishable for the two groups, overlapping 95% CIs). Finally, we examined how the relation between stimulation amplitude and response magnitude is affected by the spatial relation between stimulation site and each responsive unit. With few exceptions, the most frequently observed and largest responses occur in units that are closest to the site of microstimulation and with the strongest stimulation amplitudes (Fig. 6F,G).
Discussion
Our data demonstrate that microstimulation in the human brain can drive reliable responses in a sparse subset of neurons up to a few millimeters away from the source of current. By characterizing the effects of microstimulation from several different locations, our data also demonstrate that a single neuron can be driven with excitation or inhibition using different stimulation sites. This observed repertoire of distinct stimulation effects within the same neurons suggests that future approaches may be developed capable of providing direct control of spiking activity at the single-neuron level.
Direct electrical stimulation of the human brain for diagnostic and therapeutic purposes, as well as for research, has enjoyed a long and successful history (Borchers et al., 2011; Lozano et al., 2019). However, it still remains fundamentally unknown how stimulation precisely affects neural activity. This has therefore imposed a significant limitation on how we can effectively deploy electrical stimulation to interface with and modulate brain function. This gap in knowledge is particularly acute at the level of individual neuronal responses. While there have been several studies examining the aggregate effects of large-scale macro-electrode stimulation on physiological responses, precise control of individual neuronal firing has remained elusive.
In most studies examining the effects of microstimulation on behavior, the location of one or more stimulation sites is selected based on well-documented anatomic topographies of neural function in sensory and motor cortices where neighboring neurons share similar response properties (Pantev et al., 1995; Brewer et al., 2002; Feldman and Brecht, 2005; Harrison et al., 2012). Studies using this approach in animal studies have generated psychometric functions that relate stimulation parameters to perception, movement, and other behaviors (Tehovnik et al., 2006; M. H. Histed et al., 2013). Similarly, in humans, microstimulation has been used to elicit sensations of touch (Flesher et al., 2016), sound (Fu, 2005), movement planning (Desmurget et al., 2009), and phosphenes (Dobelle and Mladejovsky, 1974; Bak et al., 1990). Although these studies often do not directly measure how microstimulation affects neural activity, they interpret behavioral effects as being driven by increased spiking activity in neurons immediately surrounding the electrode tip (Tehovnik et al., 2006).
In a more limited number of studies, microstimulation has been used to directly examine the effects of stimulation on neural activity. These studies, almost exclusively performed in animals, have mapped the spatial relationship between individual stimulating electrodes and neural spiking responses (Stoney et al., 1968; McIlwain, 1982; Butovas and Schwarz, 2003; Logothetis et al., 2010). Responses are generally sparse (M. Histed et al., 2009), and can be either excitatory or inhibitory (Butovas and Schwarz, 2003). The traditional interpretation of these data has been that larger stimulation currents result in the recruitment of larger numbers of neurons because of an expanding stimulation field (Seidemann et al., 2002; Tolias et al., 2005; Tehovnik et al., 2006). Recent work has challenged this view, however, and has instead demonstrated that larger stimulation currents increase the proportion of affected units within a fixed spatial distance from the stimulation site (M. Histed et al., 2009). Hence, the effects of stimulation on neuronal spiking responses are likely complex, which has contributed to the difficulty in developing approaches for modulating spiking activity at the level of single neurons.
Indeed, the ultimate goal of stimulation includes selective control over neuronal firing. The introduction of optogenetic manipulation, for example, has afforded animal researchers a powerful tool for exploring and modulating neural circuits at the level of individual neurons, and consequently a powerful approach for investigating the relation between neural activity and behavior (Lee et al., 2020). Selective control of neuronal firing in the human brain could similarly open powerful pathways for both clinical and research innovations. However, a key challenge for using stimulation to modulate neuronal activity at this spatial scale, and consequently function, has been that interpreting and ultimately manipulating neural coding requires considering spiking activity at the population level rather than just within individual neurons.
To address this challenge, one would ideally like to concurrently record a large sample of single units while stimulating from multiple individual sites. In animal studies, the availability of advanced tools, such as optical imaging, makes this approach feasible (Carrillo-Reid et al., 2017). The options for adopting a similar approach in the human brain, however, are more limited. Here, we use MEAs implanted in the human cortex. MEAs have been used extensively for decoding neural signals in brain machine interface applications (Collinger et al., 2013; Aflalo et al., 2015; Wodlinger et al., 2015), but relatively little work has been completed examining how stimulation through such MEAs can drive the responses of individual neurons. To circumvent some of these limitations, we commissioned the development of a new MEA design that afforded us the opportunity to simultaneously stimulate and record neural signals through mulitple microelectrode channels, and thus examine the effects of microstimulation across a population of spiking neurons. Our data demonstrate why this was critical, as we found that the likelihood of any one stimulation site affecting a single unit is quite low, particularly if that unit is more than a few millimeters from the site of stimulation. Our ability to find rare but reliable stimulation effects therefore depended on our use of an array of recording sites with multiple stimulation sites distributed throughout. Hence, the approach we adopted recapitulates some of the advantages offered by optical imaging, namely, enabling us to simultaneously record activity from a relatively large number of neurons in response to stimulation.
The use of direct electrical stimulation by both clinicians and researchers often assumes that stimulation drives widespread and stereotyped responses across the underlying neural populations (Borchers et al., 2011; Lozano et al., 2019). Our finding that microstimulation drives heterogeneous responses among individual units appears to contradict this assumption, and suggests that aggregate behavioral and physiological responses to direct electrical stimulation belie an underlying complexity manifest at the level of individual neurons. Consistent with clinical assumptions, we do find that the dominant response type is inhibitory, and the inhibitory effects of stimulation appear to spread radially, affecting a relatively larger proportion of neurons through a spreading electric field. However, we also find excitatory responses that are selective and more spatially dispersed, consistent with previous microstimulation studies that have argued that excitatory effects reflect direct activation of axonal fibers immediately adjacent to the stimulation electrode (M. Histed et al., 2009).
In both cases, neural responses are modulated by stimulation amplitude, suggesting a neural basis for previously reported psychometric dependencies on electrical stimulation current level (Bartlett et al., 2005; Murphey and Maunsell, 2008; Fridman et al., 2010; Kim et al., 2015). This is consistent with previous modeling work that has suggested that increasing stimulation amplitude leads to increasing axonal activation, but that this effect also decreases with distance (Kumaravelu et al., 2022). Although some prior studies have suggested that the relation between response and stimulation amplitude is unclear (Logothetis et al., 2010), our observation that there is a reliable relation is consistent with effects observed in both nonhuman primates and in rodents (M. Histed et al., 2009; Sombeck et al., 2022). Hence, the use of animal models continues to be valuable for understanding the effects of stimulation in the human brain. Indeed, even at the level of microstimulation, human studies have demonstrated increasing perceptual sensations with increasing stimulation amplitude (Flesher et al., 2016; Fifer et al., 2022), likely related to the underlying effects on neural activity. Together, our data and analyses provide a key stepping stone in using microstimulation to identify and precisely control the spiking activity of individual neurons across several millimeters of human cortex.
Footnotes
This work was supported by the Intramural Research Program of the National Institute of Neurological Disorders and Stroke. This work was also supported in part by the DARPA Restoring Active Memory program. The views, opinions, and/or findings contained in this material are those of the authors and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government. We thank all patients who have selflessly volunteered their time to participate in this study.
The authors declare no competing financial interests.
- Correspondence should be addressed to Kareem A. Zaghloul at kareem.zaghloul{at}nih.gov or John H. Wittig Jr at john.wittig{at}gmail.com
