Characterization and closed-loop control of infrared thalamocortical stimulation produces spatially constrained single-unit responses

Abstract

Deep brain stimulation (DBS) is a powerful tool for the treatment of circuitopathy-related neurological and psychiatric diseases and disorders such as Parkinson's disease and obsessive-compulsive disorder, as well as a critical research tool for perturbing neural circuits and exploring neuroprostheses. Electrically mediated DBS, however, is limited by the spread of stimulus currents into tissue unrelated to disease course and treatment, potentially causing undesirable patient side effects. In this work, we utilize infrared neural stimulation (INS), an optical neuromodulation technique that uses near to midinfrared light to drive graded excitatory and inhibitory responses in nerves and neurons, to facilitate an optical and spatially constrained DBS paradigm. INS has been shown to provide spatially constrained responses in cortical neurons and, unlike other optical techniques, does not require genetic modification of the neural target. We show that INS produces graded, biophysically relevant single-unit responses with robust information transfer in rat thalamocortical circuits. Importantly, we show that cortical spread of activation from thalamic INS produces more spatially constrained response profiles than conventional electrical stimulation. Owing to observed spatial precision of INS, we used deep reinforcement learning (RL) for closed-loop control of thalamocortical circuits, creating real-time representations of stimulus-response dynamics while driving cortical neurons to precise firing patterns. Our data suggest that INS can serve as a targeted and dynamic stimulation paradigm for both open and closed-loop DBS.

Keywords: Biological; Health; and Medical Sciences: Neuroscience; closed-Loop DBS; deep brain stimulation; deep reinforcement learning; infrared neural stimulation.

Fig. 1.

Implantation and EEG-MLR procedures. A) Rodents were implanted with fiber optic optrodes into the MGB and 16 channel microwire arrays into auditory cortex. Placement of microwire array was confirmed by tonic single-unit responses evoked from 80 dB filtered Gaussian noise stimuli during implantation. B) Schematic of the four-channel EEG-MLR recording preparation. C) Schematic of the rodent auditory thalamocortical circuit. Stimulation optrodes were placed in the MGB with excitatory thalamocortical projections to layers 3–4 of primary and secondary auditory cortices. Microwire array recording electrodes were placed in layers 3–4 of primary auditory cortex confirmed during surgery by low-latency single-unit activity. D) Histological confirmation of stimulation optrode placement in the MGB. E) EEG-MLR prepostsurgical ratios show small changes in wave P1, N1, and P2 correlates of auditory thalamocortical function in amplitude and latency due to passive presence of device at 65 or 85 dB-SPL click stimuli. While changes in amplitudes and latencies were observed, differences did not rise to level of significance (P > 0.05). Rodent implantation and EEG diagrams were created using BioRender (www.biorender.com) under publication license.

Fig. 2.

A) Example INS-evoked PSTHs. BARS estimates of 1.8 mJ (solid line) and 0.92 mJ (dotted line) per pulse show higher energy pulses drive higher firing, lower latency responses. B) Bayesian hierarchical linear regression models of cortical dose-response profiles elicited by INS describe the effects of varying INS parameters on evoked firing rates. Distributions of regression parameters as estimated from observed data are given for applied laser energy (top left), laser interstimulus interval (top right), and laser energy-interstimulus interval (bottom left) interactions. Assessment of statistical significance was made using the Bayesian convention of HDI estimation. Distributions of regression parameters are shown, with parameters being significant if 95% of parameter value distribution does not contain null value 0 (vertical dotted line). Regressions show that increases in applied energy significantly increase maximum cortical firing rates with a maximum a priori estimate of 0.58 increase in log firing rate in response to increases in log energy (95% HDI does not overlap 0). The width of the 95% HDI of the energy parameter (0.27–0.88) suggests that while cortical firing rates increase with increases in laser energy generally, total increase in maximum firing rates per unit change in INS energy is dependent on the physiology of the neuron. Slight decreases in firing rate with increased ISIs were observed (MAP = −0.055), but not significant (95% HDI overlaps 0). Laser energy and ISI interactions also did not significantly change evoked cortical firing rates (95% HDI overlaps 0). Basal firing rates (bottom right) of neurons were significantly above 0 (95% HDI does not overlap 0, MAP estimate = 2.2). C) Evoked single-unit spike train MI increases as INS energy increases (standard error presented as shaded error bars).

Fig. 3.

Examples of auditory cortex firing classes evoked from INS stimuli. A) Evoked cortical firing activity was classified into onset, onset-sustained, sustained, and offset classes. Responses were classified as onset if INS induced a significant increase in firing activity above spontaneous activity compared to the 200 ms prestimulus window with a return to spontaneous firing rates before cessation of INS stimuli (top left). A response was classified as sustained if INS elicited a significant increase in firing rate above spontaneous rate that maintained firing rates of atleast $\frac{\mathrm{onset}\mathrm{rate}}{\mathrm{sustained}\mathrm{rate}}<3$ through the duration of the stimulus (top right). Responses were classified as onset-sustained if INS elicited a significant increase above spontaneous rate with sustained activity above a firing rate of $\frac{\mathrm{onset}\mathrm{rate}}{\mathrm{sustained}\mathrm{rate}}\ge 3$ (bottom left). A response was classified as offset if firing rates significantly increased from spontaneous rates after cessation of INS stimulus (bottom right). Onset, sustained, and onset-sustained responses received an additional classification of inhibition if poststimulus firing rates fell below 95% of mean prestimulus spontaneous rates. B) PCA decomposition was performed to assess if neural classes formed distinct clusters. Dimensionality reduction into the top three components of largest explained variance shows class specific clustering shapes with overlap between classes, suggesting that evoked responses are not uniquely stereotyped into fixed classes, but exist across a continuum. This is exemplified in (B), top left where PCA clustering of onset and offset classes into top two coefficients of explained variance show clusters with defined shapes and strong overlap of offset clusters into a subsection of the onset space. The decomposition of all classes into top three components of explained variance is shown in (B) top right, bottom left, and bottom right, respectively.

Fig. 4.

JPSTH analysis reveals INS thalamocortical recruitment is spatially constrained. A) JPSTH methods quantify the correlation dynamics of pairs of neurons. For each INS condition, PSTHs for two neurons are time-aligned on an orthogonal axis. A JPSTH matrix is created by counting coincidences of spikes for both PSTHs across all time normalized to the variances of both PSTHs, creating a time-varying covariance response to INS stimuli, where events occurring across the lower left to upper right (the main diagonal) correspond to correlations close in time. Then, the stimulus normalized JPSTH matrix is calculated by creating a JPSTH with one neuron shifted in time by one stimulus trial, which is subtracted bin-by-bin from the original JPSTH to remove physiologically induced correlations arising after the primary INS stimulation train. The normalized JPSTH represents long-term correlated functional connectivity after accounting for direct INS-induced responses. JPSTHs were calculated for all INS responsive neurons on the recording array. Cross-correlograms were calculated by summing coincident firing in diagonal bins oriented in the direction of the main diagonal of the JPSTH, normalized by bin length. JPSTH matrices were smoothed by a two-dimensional Gaussian smoothing kernel with variance of three for visualization. All reported correlations and JPSTHs were calculated from raw, unsmoothed covariance histograms. Statistical significance of correlations was assessed via shuffled permutation testing. B) Bar plots of maximum correlations vs. distance. Distances of 0 correspond to neurons recorded on the same electrode. INS-induced response correlations (left) show that spread of activation within cortical layers III/IV was constrained to a lateral spread of $\le 1,500\mu$ m with 90% of responses constrained to $\le 1,000\mu$ m. Laser energies $<1$ mJ per pulse (right) further limited correlated lateral spread to $\le 1,250\mu$ m with 90% of responses $\le 750\mu \mathrm{m}$. C) Pairwise normalized JPSTH correlations (left) measuring poststimulation induced connectivity show layers III/IV cortical lateral spreads limited to $\le 1,250\mu$ m with 90% of responses constrained to $\le 625\mu \mathrm{m}$. Laser energies ≤ 1 mJ per pulse (right) further show poststimulation induced connectivity limited to $\le 1,000\mu$ m with 90% of responses constrained to $\le 625\mu \mathrm{m}$.

Fig. 5.

SpikerNet, a deep RL-based closed-loop control system. A) Schematic of SpikerNet operation, which utilizes TD3 RL. The state is representative of a response as recorded from the electrode environment. The agent is the set of all safe stimulation parameters. B) SpikerNet is able to find arbitrary neural firing patterns through repeated iterations of stimulation through the environment. C) SpikerNet partakes in search and targeting behavior to find target responses and to learn stimulation parameters that best drive the neural environment to target state. D) Example evoked responses during SpikerNet search and learning show a wide variety of firing classes are evoked during algorithm search. While fits were calculated around the window of evoked activity (200–300 ms), more complex multipeaked and offset responses were observed (trials 12, 22, and 26).