Abstract
Aberrant neural oscillations hallmark numerous brain disorders. Here, we first report a method to track the phase of neural oscillations in realtime via endpointcorrected Hilbert transform (ecHT) that mitigates the characteristic Gibbs distortion. We then used ecHT to show that the aberrant neural oscillation that hallmarks essential tremor (ET) syndrome, the most common adult movement disorder, can be transiently suppressed via transcranial electrical stimulation of the cerebellum phaselocked to the tremor. The tremor suppression is sustained shortly after the end of the stimulation and can be phenomenologically predicted. Finally, we use featurebased statisticallearning and neurophysiologicalmodelling to show that the suppression of ET is mechanistically attributed to a disruption of the temporal coherence of the aberrant oscillations in the olivocerebellar loop, thus establishing its causal role. The suppression of aberrant neural oscillation via phaselocked driven disruption of temporal coherence may in the future represent a powerful neuromodulatory strategy to treat brain disorders.
Introduction
Synchronous oscillatory firing in large populations of neurons has diverse functional roles in the central nervous system (CNS), including regulation of global functional states, endowing connectivity during development, and providing spatiotemporal reference frames for processing of sensory input^{1,2}. Aberrant synchronous oscillations have been associated with numerous brain disorders^{3,4}. A palpable manifestation of such aberrant oscillation is pathological tremor in essential tremor (ET) syndrome, the most prevalent movement disorder affecting 0.4% of the general population^{5}. While the biomolecular origin of ET remains elusive, rendering pharmacological interventions unspecific and often inefficient^{6}, its systemslevel origin, i.e. oscillatory activity in the corticocerebellothalamocortical (CCTC) network, is well established^{7}. Invasive systemslevel interventions such as lesioning and highfrequency deep brain stimulation (DBS) can successfully treat medication refractory ET^{6,8}, but their widescale application is limited due to the need for brain surgery. However, such aberrant oscillations fundamentally require a delicate cascade of coherent activities across the network components. We here explored whether such a cascade of coherent activities in the CCTC under ET can be disrupted noninvasively by perturbing the synchronous activity of the cerebellum via stimulation that is phaselocked to the tremor oscillation. To phaselock the stimulation to the tremor oscillation, we first present a strategy to mitigate the Gibbs phenomenon distortion^{9} from the Hilbert transformation^{10} to compute the instantaneous phase of an oscillatory signal in realtime, a strategy that we called endpoint corrected Hilbert transform (ecHT). We then demonstrate that if transcranial alternating current stimulation (tACS) of the cerebellum is phaselocked to ET movement it can suppress its amplitude. Finally, we show that the suppression of ET amplitude is attributed to a disruption of the cascade of coherent activities in the olivocerebellar loop.
Results
Realtime computation of instantaneous phase via endpoint corrected Hilbert transform
To enable phaselocking of stimulation to oscillatory activity, we first developed a strategy to compute in realtime the instantaneous phase of oscillatory signals. Traditionally, the instantaneous phase and envelope amplitude, of a bandlimited, timevarying oscillatory signal are computed from a complexified version of the signal, known as the analytic signal, in which the real part is the unmodified signal and the imaginary part is the signal’s Hilbert transform^{10}. The discrete analytic signal is most accurately and efficiently computed in the frequency domain^{11}. However, the Gibbs phenomenon^{9} has made it impossible to accurately compute the instantaneous phase and amplitude at the ends of finitelength analytic signals^{12}. We hypothesised that by applying a causal bandpass filter to the frequency domain representation of the analytic signal we would mitigate the Gibbs phenomenon by establishing a continuity between the two ends of the signal and remove the distortion selectively from the end part of the signal—aka ecHT. See Methods for a detailed description of the ecHT.
To assess whether the ecHT strategy could effectively mitigate the Gibbs phenomenon at the endpoint of the analytic signal, we computed the Hilbert transform of a test signal, i.e. a finitelength discrete cosine waveform, and quantified the error at the endpoint. Figure 1a, b show the Fourier spectra and the Hilbert transforms without the endpoint correction when the signal completed and did not complete full cycles within the sampled time interval, respectively. At the endpoint of the signal without ecHT, the maximal phase error was 179° (mean error 47 ± 50° standard deviation, st.d.), and the maximal amplitude error was 191% (76 ± 69%), Fig. 1c. Fig. 1d shows the same as Fig. 1b but with the endpoint correction. At the endpoint, the ecHT strategy reduced the phase error by at least an order of magnitude (maximal error 12°; mean error 9 ± 2° st.d.) and the amplitude error by at least two orders of magnitude (8%; 4 ± 2%), Fig. 1e. The effects of the filter bandwidth and filter order are shown in Fig. 1f, g, respectively.
Cerebellar stimulation phaselocked to essential tremor movement
Next, we deployed the ecHT to test whether stimulation of the cerebellum phaselocked to the tremor movement can perturb ET in a cohort of 11 human participants with ET (see Supplementary Table 1 for demographic details). We measured the tremor movement of the hand, computed its instantaneous phase in realtime, generated eight different stimulating currents – sinusoidal at six different phase lags (0°, 60°, 120°, 180°, 240°, 300°), a control sinusoidal at the tremor frequency without phaselocking, and a sham, and applied them transcranially to the ipsilateral cerebellum via scalp electrodes (mean current amplitude 2.7 ± 1 st.d. mA). Fig. 2a shows a schematic of the phaselocked stimulation concept, Fig. 2b shows a schematic of the electrode configuration and the theoretical distribution of the electric fields in the brain, computed using finite element method (FEM) modelling. Supplementary Movie 1 shows a representative video. We applied each stimulation condition in a block of 60 s during which the participants maintained a tremor evoking posture. Each block consisted of a 30 s stimulation period (including 5 s of rampup and 5 s of rampdown) and 15 s stimulationfree periods before and after. We repeated the stimulation conditions four times in a doubleblinded random order with a 30 s rest interval between conditions and 5–10 min rest interval between sessions of eight stimulation conditions (see Fig. 2c for a schematic of the study design and Methods).
To assess whether the stimulating currents were delivered at accurate phaselag, we computed, offline using Hilbert transform, the lag between the instantaneous phase of the stimulation waveforms and the instantaneous phase of the tremor movement waveforms. We found that during the phaselocked stimulation, the phaselag distribution of each condition was narrow and different from the other conditions throughout the stimulation period (Fig. 2d(i)) and during the first and second halve periods (Fig. 2d(ii)), (p < 10^{−8} for all periods; Fisher test; see Supplementary Table 2 for full statistics). The difference between the measured phaselag and the set phaselag was small, i.e. 3 ± 11° (mean ± st.d), across all the phaselocked conditions. The mean resultant vector length (quantifying the circular spread)^{13}, was close to one, i.e. 0.98 ± 0.01, across all the conditions, and did not differ between conditions throughout the stimulation period (Fig. 2e(i)), and during the first and second halve periods (Fig. 2e(ii); p > 0.95 for all periods; oneway ANOVA, see Supplementary Table 3 for full statistics). The mean resultant vector length was slightly larger at stimulation blocks with higher tremor amplitude (Fig. 2f(i)) and was slightly smaller at stimulation blocks with higher tremor amplitude st.d. (Fig. 2f(ii)) higher tremor frequency (Fig. 2g(i)) and higher tremor frequency st.d. (Fig. 2g(ii)). In contrast, during the sinusoidal stimulation without phaselocking, the phaselag distribution was not different from a uniform distribution (Fig. 2d(iii); p > 0.4 for all periods; Omnibus test). The mean resultant vector length was small, i.e. 0.19 ± 0.071, and did not differ from sham stimulation (p = 0.37, paired Wilcoxon signedrank test), indicating that the stimulation did not entrain the tremor phase (Fig. 2d, e and Supplementary Table 3). Across all stimulation conditions, the mean resultant vector length was not different in trials in which participants reported sensation underneath the electrodes and trials in which no sensation was reported (p = 0.3, Paired signrank test).
Phasedependent suppression of essential tremor amplitude
After establishing that the stimulating currents were delivered at the desired phase lags, we assessed whether they affected the tremor amplitude. To quantify the stimulation effect relative to the baseline period and relative to the effect of sham stimulation, we computed, for each participant, the zscore of the tremor amplitude relative to the mean and the st.d. of the tremor amplitude during baseline in each stimulation condition, and then subtracted the median zscore of the tremor amplitude during sham stimulation (there was no significant difference in the tremor frequency and amplitude during baseline between conditions, see Supplementary Table 1 for full statistical details). To examine the temporal dynamics of the effect we quantified the zscore values during the first half and second half of the stimulation period, as well as during the poststimulation period.
We found that the stimulation at the tremor frequency without phaselocking resulted in a tremor amplitude reduction, yet not statistically significant (Fig. 3a). A significant tremor amplitude change (reduction or increase) occurred in only a small number of participants (Fig. 3b and Supplementary Table 4). Across these subsets of participants, the change was statistically significant only in those showing a reduction and only during the first half of the stimulation (Fig. 3c, d). The corresponding percentage reduction during the first half period of the stimulation was −10.8 ± 3.0% (mean ± st.d.) relative to baseline. In contrast, stimulation that was phaselocked to the tremor movement resulted in a significant reduction in the tremor amplitude, that increased throughout the stimulation period and sustained during the poststimulation period (Fig. 3e; see Supplementary Fig. 1 for zscore values expressed relative to stimulation without phaselocking). The number of participants who showed a significant reduction in the tremor amplitude was significant during the second half of the stimulation and the poststimulation period, while the number of participants who showed a significant increase in the tremor amplitude was not significant throughout (Fig. 3f and Supplementary Table 4; p value threshold of amplitude change was Bonferroni corrected for six phaselocked conditions). Across these subsets of participants, the reduction/increase in the tremor amplitude was statistically significant throughout (Fig. 3g, h). The corresponding percentage reduction (and increase) in tremor amplitude during the first half period of the stimulation, second half period of the stimulation, and after the stimulation period, was −18.1 ± 2.5% (8.3 ± 4.5%), −15.2 ± 2.2% (1.6 ± 2.0%), and −12.0 ± 2.3% (6.5 ± 3.3%), respectively, relative to baseline. The change in tremor amplitude was not different between sessions (p = 0.64, ANOVA; p = 0.32, linear mixed effect model with sessions as a predictor variable). Across all stimulation conditions, the zscore tremor amplitude was not different in trials in which participants reported sensation underneath the electrodes and trials in which no sensation was reported (p = 0.54, paired ttest).
Comparing the phaselocked conditions, we found that the reduction in tremor amplitude was close to significance (not corrected) only at a phaselag of 0° (Fig. 4a) but the number of participants who showed a significant reduction in tremor amplitude was not significant (Fig. 4b). However, if the phase lags of individual participants were expressed relative to the phase lag that resulted in the largest reduction of their tremor amplitude, the reduction in tremor amplitude and the number of participants who showed a significant reduction, were statistically significant, indicating a narrow range of efficacious phase that can vary between participants (Fig. 4c, d, see Supplementary Table 5 for complete statistical details). The corresponding percentage reduction during the second half period of the stimulation at 0° phaselag was −21.5 ± 4.2% relative to baseline.
To test whether the effect of the stimulation on the tremor amplitude is reproducible, we repeated the experiment in a subset of participants (n = 6, including participants 1,2,3,6, and 11 who showed a reduction in the tremor amplitude and participant 9 who did not; see Supplementary Table 1 for demographic and clinical details during the repeated experiment) and analysed the data in the same way as in the original experiment. We found that in the repetition experiment the stimulation currents were delivered at the same phaselag accuracy as in the original experiment (Supplementary Table 6). As before, stimulation at the tremor frequency without phaselocking resulted in a tremor amplitude reduction, yet not statistically significant (Fig. 4e), however stimulation currents that were phaselocked to the tremor movement resulted in a significant reduction in the tremor amplitude that was sustained during the poststimulation period (Fig. 4f). The participants who showed a significant reduction in the tremor amplitude during the stimulation period in the original experiment also showed a significant reduction in the tremor amplitude in the repetition experiment (see Supplementary Table 7 for full statistics). The zscore reduction in the tremor amplitude across those participants was not different from the original experiment (Fig. 4g). Comparing the phaselocked conditions, we found that across the cohort the reduction in the tremor amplitude was smaller at phaselags of 0° and 300° (Fig. 4h, see also Supplementary Table 8 for full statistics). Within individual participants the phaselag values that reduced the tremor amplitude were consistent in only 20% of the cases.
Prediction of participants’ response from distinct features of the tremor movement
Next, we sought to explore whether the variability in the participants’ response to the stimulation can be attributed to certain characteristics of their ET condition. We divided the participants into two groups, i.e. a ‘responder’ group (n = 7, including participants 1, 2, 3, 6, 8, 9, and 11) and a ‘nonresponder’ group (n = 4, participants 4, 5, 7, and 10). A participant was defined a ‘responder’ if his/her tremor amplitude decreased in at least one of the tested stimulation phases relative to sham and did not increase in any of the tested stimulation phases relative to sham, and a ‘nonresponders’ if his/her tremor amplitude increased in at least one of the tested stimulation phases relative to sham or did not change in any of the tested stimulation phases relative to sham. We first assessed whether certain clinical or demographic characteristics can distinguish between responder and nonresponder groups but found only nonsignificant trends of younger age (p = 0.07, Wilcoxon ranksum test) and higher tremor frequency (p = 0.08) in responders (see Supplementary Table 1 for full statistical details). In addition, we did not find a difference between the groups in the amplitude of the applied currents (p = 0.8).
We then explored whether certain characteristics of the tremor movement can distinguish between the two groups. We deployed a featurebased statistical learning strategy^{14} to extract 7873 different timeseries features from a 10 s segment of the tremor movement before the onset of the stimulation in all the trials with phaselocked stimulation (301 trials in total, including 28 trials per participant except participant 3 in which only 21 trials were recorded); exemplary tremor traces are shown in Fig. 5a. We then used the features and a support vector machine (SVM) with a linear kernel to classify the tremor trials according to the subjects’ responsiveness to a phaselocked stimulation. We found that using all the features, the tremor trials could be classified according to the participants’ response with an accuracy of 97% (Fscore of 96). However, even a small number of features was sufficient for high accuracy classification, using the top 1, 5, 10, and 40 features with highest singlefeature classification accuracy, the tremor trials could be classified with an accuracy of 83%, 81%, 86%, and 92% (Fscore of 82, 80, 85, and 91), respectively (Fig. 5b).
We then used a hierarchical cluster tree approach to search for the most informative features among the 40 features with the highest classification accuracy (Fig. 5c; feature values of individual participants did not differ between trials, p > 0.5; ANOVA). We identified 14 clusters of correlated features and extracted the corresponding features at the centre of those clusters—the list of the most informative features is given in Supplementary Table 9 and the magnitude probability density plots of exemplary features are shown in Fig. 5d (the classification accuracy plateaued at ~14 features, Fig. 5e). The extracted features revealed that the tremor movement in responders was smaller (Fig. 5dii), had a more sinusoidal like regularity (Fig. 5diii and Fig. 5div), and had a higher amplitude symmetry relative to zero (Fig. 5di). The Euclidean distance between feature centroids of the responders class and nonresponders class was 0.55 (feature centroid of a class was computed by averaging the features across the corresponding samples). The feature centroids of individual participants who responded to the stimulation located at a distance <0.5 to the feature centroid of the responders class and had a longer distance to the feature centroid of the nonresponders class (exception was participant 8; Fig. 5f; distance of responders to responders’ class, mean 0.35 ± 0.2 st.d.; responders to nonresponders class, 0.6 ± 0.25; nonresponders and responders class, 0.65 ± 0.15; nonresponders and nonresponders class, 0.35 ± 0.15).
To test whether these features of the tremor movement can potentially help to predict the response of participants to the stimulation, we repeated the experiment in a new cohort of seven human participants with ET. We analysed the data in the same way as in the original cohort and extracted the same 14 features from the 10 s tremor movement before the stimulation onset (see Supplementary Table 10 for demographic details, see Supplementary Table 11 for phaselocking and Supplementary Table 12 tremor amplitude statistics). We found that three participants (i.e. participants 2,3, and 7) responded to the stimulation based on the aforementioned responding criterion. The feature centroids of these participants, but not the rest of the cohort, were located at ≤0.5 distance to the feature centroid of the responders class from the original cohort and had a longer distance to the feature centroid of the nonresponders class from that cohort (Fig. 5g) indicating a consistency in the relationship between the features of the tremor movement and the response to the stimulation.
Suppression of essential tremor amplitude is underpinned by disruption of temporal coherence of movement
After establishing that participants who responded to stimulation had distinct characteristics of tremor movement during baseline, we next sought to explore whether the change in tremor amplitude during stimulation was associated with a change in other characteristics of tremor movement. We divided all the tremor trials with phaselocked stimulation (again 301 trials in total) into three datasets according to the change in tremor amplitude during stimulation relative to sham, i.e. trials with a decrease in tremor amplitude (‘decrease’; 58 trials from 11 subjects), trials with an increase in tremor amplitude (‘increase’; 51 trials from 10 subjects; participant 6, did not show an increase in tremor amplitude in any phaselocked condition), and trials without a change in tremor amplitude (‘nochange’; 192 trials from 11 subjects).
We then deployed the same featurebased statistical learning strategy^{14} to test whether the characteristics of the tremor movement can distinguish between the stimulation and baseline periods in these three datasets. We extracted the same 7873 features as before from a 10 s segment of the tremor movement before the onset of the stimulation and from a corresponding 10 s segment during the middle of the stimulation; exemplary tremor traces with tremor amplitude ‘decrease’ and ‘increase’ are shown in Fig. 6a, b, respectively. We then used the features and the same SVM as before to classify the tremor trials according to the period class, i.e. ‘baseline’, or ‘stimulation’. We found that the ‘decrease’ dataset had a higher probability of classification with high accuracy compared to the ‘increase’ and the ‘nochange’ datasets (Fig. 6c; ‘decrease’ vs. ‘increase’, p = 0.01; ‘decrease’ vs. ‘nochange’, p = 0.008; ‘increase’ vs. ‘nochange’, p = 0.45; and against a null distribution, generated by assigning random values to the feature), ‘decrease’, p = 0.005; ‘increase’, p = 0.34; ‘nochange’, p = 0.58; pairwise Kolmogorov–Smirnov test).
Focusing on the ‘decrease’ dataset, we found that using all the features, the tremor trials during stimulation and baseline could be classified with an accuracy of 79% (Fscore of 79). However, the classification accuracy was dominated by only a few features, using the top 1, 5, 10, and 40 features with highest singlefeature classification accuracy, the tremor trials could be classified with an accuracy of 78%, 79%, 79%, and 80% (Fscore of 78, 81, 81, and 81, respectively; Fig. 6d). We then used, as before, the hierarchical cluster tree approach with a between feature correlation threshold of 0.2 to search for the most informative features among the 40 features with the highest classification accuracy (Fig. 6e). We identified nine clusters of correlated features and extracted the corresponding features at the centre of those clusters—the list of the most informative features is given in Supplementary Table 13 and the magnitude probability density plots of the central features with the highest probability are shown in Fig. 6f. We found that the classification was dominated by two timeseries features, i.e. the ‘information gain’ feature, which estimates how easy it is to predict a data point in the time series from the preceding data points, and the ‘quadratic fit of power spectrum cumulative sum’ feature, which characterises the power spectrum of the time series. The increase in ‘quadratic fit of power spectrum cumulative sum’ during stimulation can be simply attributed to the drop in the spectral peak at the tremor’s frequency. In contrast, the increase in ‘information gain’ during stimulation revealed a loss of linear dependency between consecutive data points of the tremor movement, i.e. a loss of temporal coherence.
To specifically test whether the change in the tremor amplitude was associated with a change in temporal coherence, we computed the change in the magnitude squared coherence during the stimulation period relative to the baseline period in the ‘decrease’ and the ‘increase’ datasets as well as in a dataset consisting of all the trials with sham stimulation (‘sham’). We found that the temporal coherence in the tremor frequencyband decreased in the ‘decrease’ dataset and increased in the ‘increase’ dataset during the stimulation, however, it did not change in the ‘sham’ dataset (Fig. 6g). The change in the tremor amplitude in the ‘decrease’ dataset, but not in the ‘increase’ dataset, was correlated with the change in the tremor temporal coherence. The change in the tremor amplitude in the ‘sham’ dataset was also positively correlated with the change in the tremor temporal coherence, however, with a smaller slope of the linear regression (Fig. 6h; combined dataset, line yintercept c = 0.2, line slope m = 1.2, R^{2} = 0.32; ‘decrease’ dataset, c = −1.4, m = 1.35, R^{2} = 0.49; ‘increase’ dataset, c = 0.94, m = 0.58, R^{2} = 0.004; ‘sham’ dataset, c = −0.3, m = 0.78, R^{2} = 0.32; Pearson correlation; see Supplementary Fig. 2 for a correlation analysis of trials during stimulation without phaselocking). The change in temporal coherence in the ‘decrease’ dataset was correlated with the onset of the stimulation and was maintained during the duration of the stimulation (Fig. 6i).
To explore the possible mechanism by which the disruption of the temporal coherence could result in a suppression of the tremor amplitude, we simulated the CCTC network under ET condition^{15} and phaselocked cerebellar stimulation. We found that the mechanism might be related to the suppression of the aberrant complex spikes in the Purkinje cells (PCs) of the cerebellum due to synchronisation of the hyperpolarizing phase of the stimulating with the onset of the complex spikes. See ‘Neurophysiological model’ in Supplementary Information.
Discussion
In this paper we presented the ecHT strategy to compute the instantaneous phase of oscillatory signals in realtime and validated it using both simulation and measurements with pathologic oscillatory brain activity, i.e. ET. The ecHT strategy is based on the application of a causal bandpass filter to the DFT of the analytic signal to mitigate the distortion, known as the Gibbs phenomenon, from its end. Other frequencydomain and timedomain filters have been previously proposed to mitigate the Gibbs phenomenon from finite signals with a discontinuity^{16} but these filters restore the DFT only away from the discontinuity itself^{17}. There have also been reports of restoring the endpoint of the analytic signal using recursive models, such as autoregression^{18} or polynomial fitting^{19} to forward predict the physiological signal so that the last acquired data points are shifted from the window edge before the computation of the Hilbert transform. Recursive models have been recently tested for phaselocking brain stimulation^{18,20,21}, showing in some cases large st.d. (e.g. ~54°)^{21} and dependency on the coherence of the signal^{18}. Ultimately, the high runtime complexity of recursive models (e.g. autoregression has a runtime complexity of O(n^{3}) for n samples, governed by the parameter estimation operation^{22}) limit their use in applications that require realtime computation using conventional, and/or portable digital hardware.
In comparison, the ecHT is a simple, yet powerful method to accurately compute the Hilbert transform in realtime to track the instantaneous phase and envelope amplitude of an oscillatory signal. The ecHT maintains the same runtime complexity as the original Hilbert transform (i.e. O(nlog(n)) for n samples), allowing implementation in simple and portable hardware. Future studies may be able to improve the accuracy of the ecHT by adjusting, online, the central frequency of the bandpass filter to the instantaneous frequency of the signal, computed e.g. via a time derivative of the instantaneous phase. Given the widespread use of the Hilbert transform to compute the instantaneous attributes of oscillatory signals^{10}, the possibility for realtime computation using ecHT opens exciting opportunities in neuroscience and beyond (e.g. to monitor rotating engines and structural defects^{23}, speech analysis^{24}, and geophysics^{25}).
We then used the ecHT to demonstrate the causal role of synchronous cerebellar activity in human participants with ET. By deploying for the firsttime phaselocking stimulation to the cerebellum, we showed, in a doubleblinded, sham and active controlled experiment, that ET amplitude can be efficiently supressed within a few seconds. The range of phases that were efficacious in suppressing the tremor in our stimulation was small but varied between participants and within participants between days of experiments perhaps due to differences in the electrodeskin capacitance. Future studies may be able to adjust the target stimulation phase online using similar closedloop strategies currently deployed to adjust the target amplitude or frequency of DBS^{26}. Our results exemplify the importance of accurate phaselocking to successfully induce a reduction in tremor amplitude. The fact that the tremor amplitude continued to drop during the stimulation period suggests that a longer stimulation period may yield an even larger suppression. The sustained drop in tremor amplitude after the end of the stimulation period may hold potential for a therapeutic effect via neural plasticity. To start testing the reproducibility of the stimulation effect, we validated the effect in a subset of participants a few years after the initial experiment and share the phaselocking methodology to allow other researchers to easily reproduce the experiment.
The rational of targeting the cerebellum in ET has been motivated by the recent discoveries of cerebellar abnormalities in ET patients and its strong connectivity to the basal ganglia (via the thalamic nuclei)^{27}. Invasive phaselocked DBS of the thalamic Vim near the region receiving input from the cerebellum showed benefit in ET^{28}. Nevertheless, numerous noninvasive cerebellar stimulation studies have failed to demonstrate a clear effect on ET severity even after multiple days of stimulation (see recent reviews^{27,29,30}). For example, a prior study applying tACS to the cerebellum, but without phaselocking, found only a phase entrainment of the tremor with no effect on its amplitude^{31}. There has been an original report that showed that noninvasive phaselocked stimulation of the motor cortex can ameliorate tremor in Parkinson’s disease (PD) patients^{32}. Although both ET and PD are caused by aberrant oscillations in the motor system, their anatomical origins and degree of coupling between the central oscillators are very distinct^{33}. Of course, the effect of stimulation on the activity of a brain circuit is complex, involving mixtures of local activation and inactivation pathways and interactions with downstream and upstream brain regions^{34}, and hence cannot be extrapolated across brain locations, brain states and diseases^{27}. In fact, even a small change in stimulation parameters was shown to result in different and sometimes opposite effects^{35,36} which may be particularly true in the case of the cerebellum given its both inhibitory and excitatory effects on the motor cortex^{37,38}. There has also been a report that a periodic stimulation of the motor cortex at the tremor frequency without phaselocking, can entrain the phase of ET in patients undergoing DBS with an efficiency that was correlated to the somatosensory sensation underneath the electrodes^{39}. In our study, the changes in the circular phase distribution and amplitude of the tremor were not dependent on the subjective sensation of the patients.
Finally, we showed, using datadriven statistical learning approach, that ET severity is linked to the temporal coherence of the movement, and that stimulation that disrupts the temporal coherence can reduce its severity. Hitherto investigations of the tremor coherence have focused on the correlation between two different tremor signals, such as the bilateral hand movement^{40}, intermuscular electromyography^{41}, and corticomuscular^{42}. These studies have elucidated important differences between diseases (e.g. ET vs. PD) however have not found a relationship to the severity of the tremor. The causal relationship between the amplitude of ET and its temporal coherence provides an important insight into the dynamics of the central oscillator underlying the disease. This is particularly interesting given the distinct relationship between the instantaneous frequency of ET and its fluctuation^{43}.
With almost a third of ET patients discontinuing medications due to insufficient benefit, medical contraindications, or the emergence of adverse effects^{44}, there is a pressing need for a novel treatment strategies for ET. Invasive DBS of the Vim nucleus is an alternative treatment for drugrefractory ET patients however, it is limited by the need for a brain surgery and the development of adverse side effects such as dysarthia and dysphagia^{6,45,46}. Our results may provide the foundation for a new interventional strategy for ET. The mechanism of action of such an interventional strategy will be based on an active disruption of the cascade of coherent activities that generate the tremor oscillation in the olivocerebellar loop. Our computational modelling suggests that it may be attributed to a timely perturbation of the generation of complex spikes in the PCs. Future computational studies may be able to explain the underlying mechanisms of those features predicting the stimulation outcome. Such a mechanism of action differs from the existing Vim DBS therapy for ET that masks the tremor oscillation in the thalamocortical (TC) loop but does not mitigate its generation in the olivocerebellar loop^{15}. Future studies with larger patient cohorts and longer stimulation periods, are needed to better pinpoint the magnitude and duration of the tremor reduction and to assess the safety profile. In the future, neuromodulatory strategies that target the temporal coherence of the pathology may offer new opportunities to treat a wide range of brain disorders underpinned by aberrant synchronous oscillations.
Methods
Endpoint corrected Hilbert transform (ecHT)
A discrete analytic signal is most accurately and efficiently computed by deriving the discrete Fourier transform (DFT) of the signal, zeroing the Fourier components of the negative frequencies and doubling the ones of the positive frequencies, and constructing the analytic signal using the inverse discrete Fourier transform^{11}. However, Gibbs phenomenon distortion^{9} in the derivation of the analytic signal at the ends of finitelength signals has rendered an accurate computation of the instantaneous phase and envelope amplitude at the last data point impossible^{12}. Since the Gibbs phenomenon stems from a nonuniform convergence of the DFT at a discontinuity between the beginning and the end of the analytic signal^{47}, we hypothesised that by applying a causal bandpass filter to the DFT of the analytic signal we would establish a continuity between the two ends of the signal and remove the distortion selectively from the end part of the signal. The bandpass feature of the filter reduces extraneous DFT coefficients, limiting the oscillatory properties to the target frequencyband, while balancing the phaselag introduced by the lowpass component of the filter with the phaselead introduced by the highpass component of the filter. The causality feature of the filter restores the linear increment of the phase at the end of the analytic signal by projecting the oscillatory properties from the adjacent, nondistorted data points. Since the DFT treats finite sampled signals as if they were replicated periodically, the projection of the oscillatory properties would continue through the beginning of the signal, thus forcing a continued increment of the phase from the restored signal end to its beginning. The runtime complexity of the filtering is O(n/2), where n is the number of frequency points, is lower than \(O\left( {n \cdot log\left( n \right)} \right)\) of the fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT) that dominates the computation of the analytical signal.
Simulation of ecHT
Simulation of ecHT was done in MATLAB (MathWorks Inc). A discrete oscillatory test signal
was generated (i being the signal number) over a finite time interval T, where 0 < n < N − 1 was the time point number and N was the total number of time samples, A_{i} was the envelope amplitude of the signal, f_{i} was the frequency of the signal, and \(\emptyset _i\) was the phase delay of the signal. The analytic signal was computed by first computing the Fourier representation Y_{i}[k] of the signal using MATLAB’s fast FFT function (‘fft’), where 0 < k < K − 1 was the frequency bin number and K was the total number of frequency samples. Then, generating the Fourier representation Z_{i}[k] of the analytic signal by zeroing the Fourier components of the negative frequencies and doubling the Fourier components of the positive frequencies, i.e.
If ecHT was applied, the Fourier representation of the analytic signal Z_{i}[k] was multiplied with the response function σ[k]of a Butterworth bandpass filter that was obtained using MATLAB’s frequency response of digital filter function (‘freqz’) from the filter’s impulse response coefficients generated using MATLAB’s Butterworth filter design function (‘butter’). Finally, the analytic signal z_{i}[n] was computed from its Fourier representation Z_{i}[k] using MATLAB’s IFFT function (‘ifft’). The phase of the signal at the last data point was computed via \(atan\left( {\frac{{imag\left\{ {z_i\left[ N \right]} \right\}}}{{y_i\left[ N \right]}}} \right)\), where \(imag\left\{ {z_i\left[ N \right]} \right\}\) is the imaginary part of the analytic signal, i.e. the Hilbert transform of the original signal, and was compared to the actual phase of the signal at the last data point, i.e. \(2\pi f_iN  \emptyset _i\). The amplitude of the signal at the last data point was computed via \(\sqrt {imag\left\{ {z_i\left[ N \right]} \right\}^2 + y_i\left[ N \right]^2}\) and was compared to the actual amplitude of the signal at the last data point, i.e.A_{i}.
Feasibility study of cerebellar electrical stimulation phaselocked to ET
Ethics
The study was approved by the local research ethics committee in accordance with the declaration of Helsinki. All participants provided written informed consent prior to study participation. Specifically, the study was approved by the Heath Research Authority (HRA; REC 03/N018, principal investigator John Rothwell, UCL). The approval included the assessment of governance and legal compliance, undertaken by HRA, with the independent Research Ethics Committee (REC) opinion provided by the National Hospital for Neurology and Neurosurgery and the UCL Institute of Neurology (ION) Joint REC. The overarching aim of the research project was to use transcranial brain stimulation paradigms to discover mechanisms of cortical excitability and their impact on motor behaviour. The research project was not classified as clinical trial or interventional trial by the HRA and hence did not required registration (which is mandatory for all clinical trials in the UK).
Participants
Eleven human participants with ET (three females) were recruited from the outpatient department of the UK National Hospital of Neurology and Neurosurgery, London. All participants fulfilled the diagnostic criteria for ET according to the Tremor Investigation Group and consensus statement of the Movement Disorder Society^{48} and were on a stable treatment regime for their tremor for at least 30 days prior to the experiment. See Supplementary Table 1 for demographic and clinical information. Experiments were performed after overnight withdrawal of tremor medication during a single study visit in the dominant hand, or in case of slight asymmetry in the hand with the larger tremor amplitude. There were no dropouts or adverse events noted.
Participants (second cohort)
Seven human participants with ET (four females) were recruited as in the original to test whether their response can be predicted via the featurebased approach developed in the original study. See Supplementary Table 10 for demographic and clinical information. Experiments were performed as in the original cohort.
Experiment design
The experiment consisted of eight stimulation conditions, i.e. six sinusoidal stimulating currents that are phaselocked to the tremor movement at different phase lags (i.e. 0°, 60°, 120°, 180°, 240° and 300°), a control sinusoidal current at the tremor frequency but without phaselocking, and a sham stimulation condition. Each stimulation condition was applied in a block (i.e. trial) of 60 s during which the participants sat in an armchair and were instructed to maintain a tremor evoking posture, i.e. stretched, elevated arm with fingers parted, while their tremor movement was measured (see details below). The 60 s block included a 15 s of a baseline period, a 30 s of a stimulation period (including 5 s of rampup and 5 s of rampdown at the beginning and end of the stimulation, respectively) and a 15 s of poststimulation period. In sham stimulation blocks, the current was set to zero after the 5 s of rampup. Each 60 s block was preceded by a short (~4 s, 2048 data samples) calibration recording also in a tremor evoking posture to compute the tremor frequency and amplitude at the onset of the block (see details below). The eight stimulation conditions were applied consecutively with a 30 s rest interval between conditions. The sequence of eight stimulation conditions was repeated four times (apart from one participant in which they were applied three times due to fatigue) in a random order with 10 min rest period between sequences. The rest interval between conditions and the rest period between sequences were occasionally extended slightly if the participants requested.
Measurement and realtime computation of instantaneous tremor phase via ecHT
Tremor movements were measured using a 3axis analog microelectromechanical system accelerometer (MMA7361, Freescale Semiconductor, Inc.; operated at a sensitivity range of ±1.5 G) that was attached to the proximal phalangeal segment of the middle finger using a custommade adapter. The 3axis acceleration measurements were sampled using three analogtodigital converters of a microcontroller (Arduino Due with an Atmel AT91SAM3X8E processor and a single ARM Cortex M3 core; operated at a clock rate of 84 MHz) at a rate of ~500 Hz and an amplitude resolution of 12bit, and the vector amplitude sum of the three axes was computed and stored in a running window of 128 samples. The instantaneous phase and amplitude of the tremor movement, i.e. at the last sample of the running window, were computed in realtime and at the same rate, using ecHT that was implemented on the microcontroller. The ecHT implementation had a 2nd order Butterworth bandpass filter (2nd order low pass, 2nd order high pass) with a bandwidth that was equal to half the frequency of the tremor and was centred at the frequency of the tremor. The frequency of the tremor was computed using FFT from a short calibration measurement of 2048 samples (i.e. frequency resolution of ~0.25 Hz) before each 60 s stimulation block. The sampled tremor movement measurement was logged to a laptop, together with the ecHT setting and the tremor frequency and amplitude computed during calibration, using a Processing script that was also used to interface with the microcontroller.
Transcranial stimulation of ipsilateral cerebellum
Sinusoidal stimulating currents were generated by first producing voltage waveforms, pseudodifferentially via two digitaltoanalog converters of the microcontroller (with an amplitude range of ±1 V and an amplitude resolution of 12bit) and then feeding them to an isolated biphasic current source (DS4, Digitimer Ltd; operated at an input range of ±1 V and an output range of ±1 mA or ±10 mA). The frequency of each voltage waveform was equal to the frequency of the tremor computed before each 60 s stimulation block as mentioned above. To phaselock a stimulating current to the ongoing tremor movement, the phase of the voltage waveform was adjusted, at the same rate of 500 Hz, to maintain a fixed phase lag to the computed phase of the last acceleration sample. The amplitude of the stimulating currents was 2.7 ± 1 mA (mean ± st.d.) across the participants, (the amplitude was individually adjusted for each participant below any discomfort level due to extraneous somatosensory stimulation underneath the electrodes). To reduce risk of extraneous highfrequency stimulation due to low signaltonoise level, the amplitude of the voltage waveform was set to zero when the amplitude of the last acceleration sample was <1% of the amplitude during the short calibration measurement before each 60 s stimulation block. The generated stimulating voltage waveforms were logged to a laptop together with the tremor movement measurements using the same Processing script.
The stimulating currents were applied transcranially to the ipsilateral cerebellum via a 2 × 2 cm^{2} skin electrode (Santamedical, 2″ × 2″ carbon electrode pad with Tyco gel that was cut to the specified dimensions) that was placed 10% nasioninion distance lateral to inion (i.e. above the cerebellar lobule VIII) and was paired with a 5.08 × 5.08 cm^{2} skin electrode (the same carbon electrode pad but was not cut) that was placed over the contralateral frontal cortex between F3–F7 or F4–F8 of the international 10–20 system. Before the placement of the electrodes, the scalp skin was prepared using 80% Isopropyl alcohol and an abrasive skin gel (NuPrep, Weaver and Company Inc), and a conductive paste (Ten20, Weaver and Company Inc) and/or a conductive gel (CG04 Saline base Signa gel, Parker Laboratories Inc) was deposited at the target locations. The resistance between the electrodes was maintained below 8 kOhm.
Analysis of stimulation phase lag
Analysis of the stimulation phase lag was done in MATLAB. The tremor movement trace of each 60 s block was filtered with the same filter settings that were used in the realtime computation, i.e. a 2nd order Butterworth bandpass filter with a bandwidth that was equal to half the frequency of the tremor and centered at the frequency of the tremor computed and logged at the short calibration period preceding each block. The instantaneous phase of the stimulating waveform trace and the instantaneous phase of the filtered tremor movement trace were computed using MATLAB’s ‘hilbert’ function, and the instantaneous phase lag between the two traces was calculated and then epoched in intervals of 1 s. The stimulating trace in the sham condition was a virtual sinusoidal waveform at the tremor frequency.
The statistics and statistical tests of the phase lag values were computed, using MATLAB’s CircStat toolbox^{13}, in the following periods—the whole stimulation period (20 s since 5 s rampup time and the 5 s rampdown time at the beginning and the end were excluded, respectively), the first half of the stimulation period (10 s since 5 s rampup time was excluded), the second half of the stimulation period (10 s since 5 s rampdown time was excluded). First, the unimodality of the phase distribution of each stimulation condition was validated using Watson’s test against a von Mises distribution (set phase 0°, p < 10^{−5}; 60°, p < 10^{−5}; 120°, p < 10^{−5}; 180°, p < 10^{−5}; 240°, p < 10^{−5}; 300°, p < 10^{−5}; no phaselock, p = 0.6). The phase distribution during stimulation with phaselocking was not different from von Mises distribution but since the phase distribution during stimulation without phaselocking was different from von Mises distribution, we used nonparametric statistical tests. Next, the circular spread of the phase distribution of each stimulation condition was quantified by computing the length of the mean resultant vector R and its uniformity was assessed using the Omnibus test. Then, the difference between the mean phase of the stimulation conditions was assessed using Fisher test and the difference between the mean resultant vector length R of the stimulation conditions was assessed using ANOVA with posthoc analysis using Wilcoxon signedrank test. Finally, the effect of the tremor parameters, i.e. amplitude and frequency, on the length of the mean resultant vector R was assessed via Pearson correlation.
Analysis of change in tremor amplitude
Analysis of the tremor amplitude was done in MATLAB. The tremor trace of each 60 s block was filtered as in the ‘Analysis of stimulation phase lag’. The instantaneous amplitude was computed using MATLAB’s ‘hilbert’ function and was epoched in intervals of 1 s. To express the tremor amplitude relative to the amplitude of the baseline period, the amplitude value of each epoch was zscored by subtracting the mean value during the baseline period and then dividing by the st.d. of the value during the baseline period. The statistics and statistical tests of the tremor amplitude values were computed in the following periods—the baseline period (10 s between 3 s and 13 s from block onset), the whole stimulation period (as in ‘Analysis of the stimulation phase lag’), the first half of the stimulation period (as in ‘Analysis of the stimulation phase lag’), the second half of the stimulation period (as in ‘Analysis of the stimulation phase lag’), and the poststimulation period (10 s between 3 s and 13 s from stimulation offset). To assess the change in the tremor amplitude relative to the change in the tremor amplitude during the sham stimulation condition, the zscore amplitude values during stimulation and during poststimulation periods of each stimulation condition were subtracted by the corresponding median zscore values of the sham stimulation condition.
To assess the effect of phaselocking the stimulation to the tremor movement, the change in the tremor amplitude due to stimulation with phaselocking and without phaselocking was analysed. First, the change in the tremor amplitude due to each type of stimulation, i.e. without phaselocking and with phaselocking (data from all six phaselags of stimulation was combined) was assessed across the participants in each epoch using unpaired ttest. Next, the change in tremor amplitude of individual participant due to each stimulation condition was assessed (i.e. data including four repetition trials from each phaselag of stimulation was treated separately) during stimulation and poststimulation periods using unpaired ttest as well as using surrogate distributions (i.e. 1000 zscores values with the same st.d. but zero mean value), where the p value threshold of the stimulation conditions with phaselocking (but not without phaselocking) were Bonferroni corrected for the six phase lag conditions. Then, the number of participants that showed statistically significant increase/decrease of zscore amplitude was assessed using Fisher’s exact test against the number of participants who did not show a change in the zscore tremor amplitude (participants could have a significant increase of zscore in one phaselag and a significant decrease of zscore in another phaselag). Finally, the zscore amplitude of the subgroup of subjects that showed a statistically significant increase/decrease of zscore amplitude was assessed using unpaired ttest.
To assess the effect of the phase lag value during stimulation, the change in the tremor amplitude due to stimulation with different phase lags was analysed. First, the change in the tremor amplitude due to each phaselag of stimulation was assessed across the participants during the stimulation period using unpaired ttest. Next, the change in tremor amplitude of individual participant was assessed during stimulation again using unpaired ttest. Then, the number of participants that showed a statistically significant increase/decrease of zscore amplitude was assessed using Fisher’s exact test. Finally, to account for differences in phase response across participants, the phase lags were expressed relative to the phase lag that resulted in the largest reduction in the tremor amplitude, and the change in tremor amplitude of individual participant and the number of participants with statistically significant change were reanalysed.
Prediction of participants’ response to stimulation from features of tremor movement
Dataset
Timeseries of tremor movement during the baseline period, i.e. 10 s (5000 data points) from 5 s after the onset of tremor posture till 5 s before the onset of the phaselocked stimulation, were extracted from all the recorded trials with phaselocked stimulation, resulting in a dataset of 301 timeseries trials (28 trials per participants except participant 3 in which only 21 timeseries trials were recorded). The timeseries were assigned a ‘responder’ or a ‘nonresponder’ label if the participant responded or did not respond to the stimulation, respectively. A participant was conservatively labelled as a ‘responder’ if his/her tremor amplitude significantly decreased in at least one of the tested stimulation phases relative to sham and did not significantly increase in any of the tested stimulation phases relative to sham, and was labelled a ‘nonresponders’ if his/her tremor amplitude significantly increased in at least one of the tested stimulation phases relative to sham or did not significantly change in any of the tested stimulation phases relative to sham.
Extraction of timeseries features
For each timeseries trace, 7873 features were computed using the highly comparative timeseries analysis (hctsa)^{14}, resulting in a 301 × 7873 feature matrix. The computed features included autocorrelations, power spectra, wavelet decompositions, distributions, timeseries models (e.g. Gaussian Processes, Hidden Markov model, autoregressive models), informationtheoretic quantities (e.g. Sample Entropy, permutation entropy), nonlinear measures (e.g. fractal scaling properties, nonlinear prediction errors) etc. All features with infinity or not a number (NaN) values and features with zero variance across the dataset were removed from the feature matrix, resulting in a reduced feature matrix of 301 × 6196. The value of each feature was individually normalised to the interval [0,1].
Classification
The feature space was partitioned, i.e. classified, using a linear SVM classifier, implemented with the classify function of MATLAB’s Statistics Toolbox, which returned a threshold that optimally separated the two classes, i.e. ‘responders’ and ‘nonresponders’ timeseries. The accuracy of the classification was quantified by first computing the balanced classification accuracy \(a = \frac{\mathrm{{{precision + recall}}}}{2}\), and then computing the harmonic mean of precision and recall, i.e. F_{1} score, \(F_1 = \frac{{2 \cdot {\mathrm{precision \cdot recall}}}}{\mathrm{{{precision + recall}}}}\), where precision is the fraction of true positive classified samples over the total of positively classified samples and recall is the fraction of true positive classified samples over the total true positive and false negative classified samples. The classification was performed using a tenfold crossvalidation to reduce bias and variance.
Performancebased feature selection
The univariate classification performance of each feature was evaluated against the class labels. A subset of 40 features with the highest singlefeature classification accuracy was selected. To reduce the redundancy within the subset of features, the Pearson correlation distance, d_{ij} = 1−ρ_{ij} was computed for each pair of features, where ρ_{ij} is the Pearson correlation coefficient between feature i and feature j, and a hierarchical clustering was performed using a complete linkage threshold of 0.2, resulting in clusters of features that were intercorrelated by ρ_{ij} > 0.8. The clusters of highly correlated features were then represented by the feature that was located most centrally within the cluster (i.e. at the cluster’s centre).
Featurebased prediction of participant response
The centroid of individual participants in the feature space (including the extracted 14 most informative features) was computed by averaging the feature values across the corresponding trials. The centroid of the participant class (i.e. ‘responders’ or ‘nonresponders’) in the same feature space was computed by averaging the features values across the corresponding trial dataset. The Euclidean distance between feature centroids was computed with pdist function of MATLAB.
Visualisation using principal component analysis
To facilitate visualisation of the feature space, principal component analysis (PCA) was performed. In this case, a covariance matrix was computed for the normalised set of features from which the eigenvectors and eigenvalues were extracted. Each principal component was constructed as a linear combination of the initial features. The first two principal components were then used to display 2D scatter plots of the features.
Change in features of tremor movement due to stimulation
Dataset
Timeseries of tremor movement during stimulation (10 s; 5000 data points; from 10 s after the onset of stimulation till 10 s before the offset of stimulation) and during baseline (10 s; 5000 data points; same as in ‘Classification and prediction of participants’ response to stimulation’) from all trials with phaselocked stimulation (301 traces of stimulation and baseline each) were extracted and assigned a ‘stimulation’ class label or a ‘baseline’ class label, respectively. The ‘stimulation’ and ‘baseline’ timeseries were then divided into three datasets according to the change in the tremor amplitude during stimulation, i.e. ‘decrease’, traces in which the tremor amplitude decreased during stimulation relative to sham (58 timeseries of stimulation and baseline each, 11 subjects); ‘increase’, timeseries in which the tremor amplitude increased during stimulation relative to sham (51 timeseries of stimulation and baseline each, 10 subjects); ‘nochange’, traces in which the tremor amplitude did not change during stimulation relative to sham (192 timeseries of stimulation and baseline each, 11 subjects). In addition, in a subset of the analysis, the same ‘stimulation’ and ‘baseline’ tremor traces were extracted from all the blocks with sham stimulation (‘sham’; 43 timeseries of stimulation and baseline each, 11 subjects).
Extraction of timeseries features, classification, and performancebased feature selection
Same as in ‘Classification and prediction of participants’ response to stimulation’.
Temporal coherence analysis
The tremor temporal coherence vs. frequency of each tremor trace was quantified by computing the magnitude squared coherence across 1 s epochs during ‘stimulation’ period and ‘baseline’ period using MATLAB’s mscohere function with a frequency range of 0–31 Hz and a 1 Hz frequency resolution. The computed values during ‘stimulation’ were then zscored relative to the mean and st.d. of the values during ‘baseline’. The tremor temporal coherence at the tremor frequency band was quantified by computing the mean zscore across the 4–8 Hz frequency bins. The tremor temporal coherence vs. time of each tremor trace was quantified by computing the magnitude squared coherence between 1 s epoch and its preceding one during ‘stimulation’ period and ‘baseline’ period using the same MATLAB’s mscohere function, zscore the ‘stimulation’ values relative to ‘baseline’ in the same way, and then computing the mean zscore across the 4–8 Hz frequency bins. Statistical significance of magnitude squared coherence at a frequency bin was characterised for each dataset (i.e. decrease’, ‘increase’, and ‘sham’) using unpaired ttest with Bonferroni corrections for multiple comparisons of frequency bins and datasets.
Neurophysiological modelling
Model description
The CCTC network model under ET condition was simulated as in Zhang et al.^{15}. The model is available on ModelDB (http://modeldb.yale.edu/266842). It consisted of 425 singlecompartment, biophysicsbased neurons from the olivocerebellar and TC loops, including 40 inferior olivary nucleus (ION) neurons in the brainstem, 200 PCs and 20 granular layer clusters (GrL; 3 distinct neurons per cluster, 60 neurons altogether) in the cerebellar cortex, 5 glutamatergic deep cerebellar projection neurons (DCNs) and 5 nucleoolivary (NO) neurons in the dentate nucleus, 5 ventral intermediate thalamus (Vim) TC neurons, 100 pyramidal neurons (PYN), and 10 fastspiking interneurons (FSI). As in our previous study^{15}, the ET condition was simulated by reducing the conductivity and increasing the decay time of the PCs’ GABAergic currents to the DCN, which mimics the loss of GABA_{A} α_{1}receptor subunits and an upregulation of α_{2}/α_{3}receptor subunits in the cerebellum. Five instances of the model were considered and for each instance, simulations were repeated under normal condition, ET condition with no stimulation, and ET condition with stimulation of the cerebellum. Each simulation lasted 11,500 ms (integration step, 0.0125 ms). ET condition was initiated after 1000 ms and stimulation started after 1500 ms and lasted till the end of the simulation.
Hitherto computational studies of the effect of electrical stimulation on tremor activity have used a range of models ranging from a single cell with detailed biophysical and morphological representations^{49} to thousands of cells in which their activity is represented by a simplified pointmass function^{50}, revealing complimentary insights. Neural network modelling has an inevitable tradeoff between the scale and biological complexity of representation with both the size of the network and the biological complexity of individual cells affect the dynamics^{51}. We chose to use a middleground approach with detailed biophysical representation but reduced morphological representation —an approach proven to be successful in the past by us^{52} and others^{53,54}. This approach may be particularly suited for ET since neural mass or meanfield models cannot represent the complex change in spiking pattern (rather than mean firing rate) observed in ET patients^{55,56}. Furthermore, by maintaining a detailed biophysical representation of the cells, we could explore the effect of the stimulation on the interaction between the highfrequency simple spiking and lowfrequency complex spiking of PCs that has been causally linked to ET^{57}.
To simulate the cerebellar stimulation, a current I_{stim} was added to all the PCs in the model. I_{stim} was sinusoidal with a frequency that is equal to the frequency of the ET and amplitudes between 15pA evoking small subthreshold depolarisations expected in our experiment. Specifically, I_{stim} with an amplitude of 1pA induced a periodic depolarisation of ~0.5 mV amplitude in the singlecompartment PC model which is similar to the depolarisation that was induced by an extracellular electric field with an amplitude of 2 V/m, predicted from our FEM modelling of the experiment (Fig. 2b), in the multicompartment PC model (Supplementary Fig. 4a–b).
To validate that the direct response of the cerebellar cortex to the stimulating electric fields is dominated by the PCs, we simulated the response of the most abundant cell types in this region, i.e. PC and granule cell (GrC) to extracellular electric fields. To best capture the spatiotemporal dynamics, we used multicompartmental models with detailed 3D geometrical reconstruction of the PC^{58} and GrC^{59}. We exposed the cells to homogenous extracellular electric fields that were aligned with the dendritesomatic axes of the cells and quantified the induced depolarisation. As in the original study with the PC model^{58}, we removed the sodium and calcium channels from the axonal initial segment of this cell to reduce its spontaneous pacemaker activity (see Supplementary Fig. 4a–b).
The amplitude of I_{stim} was normalised to the average amplitude of the endogenous synaptic current to PCs, measured under ET state over 4000 ms (see also Perkel et al.)^{60}, with I_{stim} of 1pA equals 4% of the average endogenous synaptic current to PCs. To phase lock the sinusoidal current to the ET oscillation, first the spike count trace of the TC neurons of the Vim was computed with a temporal resolution of 1 ms and then filtered using a 2nd order Butterworth bandpass filter with cutoff frequencies of 6 Hz and 10 Hz. Then, the instantaneous phase of the spike count trace was computed online every 10 ms using ecHT on a running window of 1000 ms, and the phase of the stimulating current was adjusted at those timepoints to maintain the target phase lag.
Computation of PCs phaselocking value
The spike count trace of the PCs was computed with a temporal resolution of 1 ms (spikes were summed across PCs) and low pass filtered using a 2nd order Butterworth filter with a cutoff frequency of 30 Hz. Then, the instantaneous phases of the spike count trace and the stimulating current were computed offline using MATLAB’s ‘hilbert’ function, and the instantaneous phase lag between the two was calculated every 1 ms. The phaselocking value of each PC was computed as in Lachaux et al.^{61} and then averaged across the PCs.
Computation of Vim power spectrum density
First, the spike count trace of the TC neurons in the Vim was computed with a temporal resolution of 1 ms (spikes were summed across TC neurons). Then, the power spectral density (PSD) of the spike count trace was computed using Welch’s method with 2000 ms Hanning window and 1000 ms overlap, and normalised to the total power between 0 Hz and 25 Hz. Tremor PSD was estimated as the peak PSD at the tremor frequency band, i.e. between 4 and 12 Hz.
Computation of DCN and Vim temporal coherence
The spike trains of the DCN and TC neurons of the Vim were low pass filtered using a 2nd order Butterworth filter with a cutoff frequency of 30 Hz, and the magnitudes squared coherence were computed using MATLAB’s mscohere function with a frequency range of 0–30 Hz. Then the magnitude squared coherence in DCN and Vim during stimulation was expressed relative to baseline by subtracting the mean value during baseline and dividing by the st.d. value during baseline, i.e. zscore.
Sensitivity analysis to the model size
To explore the effect of the model size on the simulation outcome, we first repeated the simulation with a fivefold increase in the number of cells in the olivocerebellar circuit while keeping the other parts of the model unchanged, i.e. ‘Model expansion 1’. Model expansion 1 consisted of 1425 cells, including 200 ION neurons, 1000 PCs and 20 GrL clusters (60 neurons altogether), 25 DCNs, 25 NO neurons, 5 Vim TC neurons, 100 PYN, and 10 FSI. We randomised the synaptic connections between the TC neurons and the DCNs with adjusted weights (20% of the original value) due to model expansion. Then, we repeated the simulation with a fivefold increase in the number of all cells in the model i.e. ‘Model expansion 2’. Model expansion 2 consisted of 2125 cells, including 200 ION neurons, 1000 PCs and 100 GrL clusters (300 neurons altogether), 25 DCNs, 25 NO neurons, 25 Vim TC neurons, 500 PYN, and 50 FSI. We randomised the synaptic connections between the different neuron types along the olivocerebellar circuit, and between TC neurons and DCNs, with adjusted weights (20% of the original value) due to model expansion.
Transcranial electric field modelling
Finite element method (FEM) electromagnetic simulations were performed in Sim4Life V.4 (ZMT ZurichMedTech AG, Zurich), using a quasistatic ohmiccurrent solver. Electrodes were created within the platform using Sim4Life’s CAD functionalities and applied to the scalp of the MIDA anatomical head model^{62}. Dirichlet (voltage) boundary conditions were assigned to the electrodes, and tissues electrical conductivities were assigned according to the IT’IS LF database^{63}. A uniform rectilinear grid of 0.6 mm was used. The current between the electrodes was calculated integrating the current flux density on a closed surface surrounding one electrode and field magnitude were normalised to 2 mA input current.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Data availability
Source data are provided with this paper. The tremor recording datasets used in this paper are available on the Harvard Dataverse repository https://doi.org/10.7910/DVN/Z6EN2I. Source data are provided with this paper.
Code availability
The endpoint corrected Hilbert transform (ecHT) code implemented in Matlab is available as a supplementary file ‘Supplementary_Code_1’. The highly comparative timeseries analysis (hctsa) is available on GitHub https://github.com/benfulcher/hctsa. The Matlab code of the most informative features in Figs. 4 and 5 is also available as a supplementary file ‘Supplementary_Code_2’. The NEURON model of CCTC network under ET condition and phaselocked electrical stimulation is available on the ModelDB repository http://modeldb.yale.edu/266842. The FEM model of the transcranial cerebellar electrical stimulation is available on the Harvard Dataverse repository https://doi.org/10.7910/DVN/H7RHQF.
References
 1.
Borst, A. & Theunissen, F. E. Information theory and neural coding. Nature Neurosci. https://doi.org/10.1038/14731 (1999).
 2.
Llinás, R. R. The intrinsic electrophysiological properties of mammalian neurons: insights into central nervous system function. Science (80.). https://doi.org/10.1126/science.3059497 (1988).
 3.
Vanneste, S., Song, J. J. & De Ridder, D. Thalamocortical dysrhythmia detected by machine learning. Nat. Commun. https://doi.org/10.1038/s41467018028200 (2018).
 4.
Llinás, R. R., Ribary, U., Jeanmonod, D., Kronberg, E. & Mitra, P. P. Thalamocortical dysrhythmia: a neurological and neuropsychiatric syndrome characterized by magnetoencephalography. Proc. Natl Acad. Sci. U.S.A. https://doi.org/10.1073/pnas.96.26.15222 (1999).
 5.
Louis, E. D. & Ferreira, J. J. How common is the most common adult movement disorder? Update on the worldwide prevalence of essential tremor. Movement Disord. https://doi.org/10.1002/mds.22838 (2010).
 6.
Deuschl, G., Raethjen, J., Hellriegel, H. & Elble, R. Treatment of patients with essential tremor. Lancet Neurol. https://doi.org/10.1016/S14744422(10)703227(2011).
 7.
Raethjen, J. & Deuschl, G. The oscillating central network of Essential tremor. Clin. Neurophysiol. 123, 61–64 (2012).
 8.
Klein, J. C. et al. The tremor network targeted by successful VIM deep brain stimulation in humans. Neurology https://doi.org/10.1212/WNL.0b013e318249f702 (2012).
 9.
Gibbs, J. W. Fourier Series. Nature 59, 606 (1899).
 10.
Orfanidis, S. J. Introduction to Signal Processing. (Englewood Cliffs, NJ: PrenticeHall, 1995).
 11.
Lawrence Marple, S. Computing the discretetime analytic signal via fft. IEEE Trans. Signal Process 47, 2600–2603 (1999).
 12.
Bracewell, R. The Hilbert Transform. in The Fourier Transform and Its Applications 267–272 (McGrawHill, New York, 1999).
 13.
Berens, P. CircStat: a MATLAB Toolbox for Circular Statistics. J. Stat. Softw. https://doi.org/10.18637/jss.v031.i10 (2009).
 14.
Fulcher, B. D. & Jones, N. S. hctsa: a Computational Framework for Automated TimeSeries Phenotyping Using Massive Feature Extraction. Cell Syst. https://doi.org/10.1016/j.cels.2017.10.001 (2017).
 15.
Zhang, X. & Santaniello, S. Role of cerebellar GABAergic dysfunctions in the origins of essential tremor. Proc. Natl Acad. Sci. U.S.A. https://doi.org/10.1073/pnas.1817689116 (2019).
 16.
Vandeven, H. Family of spectral filters for discontinuous problems. J. Sci. Comput. https://doi.org/10.1007/BF01062118 (1991).
 17.
Gottlieb, D. & Shu, C.W. On the Gibbs Phenomenon and Its Resolution. SIAM Rev. https://doi.org/10.1137/s0036144596301390 (2003).
 18.
Chen, L. L., Madhavan, R., Rapoport, B. I. & Anderson, W. S. Realtime brain oscillation detection and phaselocked stimulation using autoregressive spectral estimation and timeseries forward prediction. IEEE Trans. Biomed. Eng. 60, 753–762 (2013).
 19.
Riviere, C. N., Rader, R. S. & Thakor, N. V. Adaptive canceling of physiological tremor for improved precision in microsurgery. IEEE Trans. Biomed. Eng. 45, 839–846 (1998).
 20.
Schaworonkow, N. et al. μRhythm Extracted With Personalized EEG Filters Correlates With Corticospinal Excitability in RealTime PhaseTriggered EEGTMS. Front. Neurosci. https://doi.org/10.3389/fnins.2018.00954 (2018).
 21.
Zrenner, C., Desideri, D., Belardinelli, P. & Ziemann, U. Realtime EEGdefined excitability states determine efficacy of TMSinduced plasticity in human motor cortex. Brain Stimul. https://doi.org/10.1016/j.brs.2017.11.016 (2018).
 22.
Schönhage, A. Equation solving in terms of computational complexity. Proc. Int. Congr. Math. 3, p. 40 (1986).
 23.
Feldman, M. Hilbert Transform Applications in Mechanical Vibration. Hilbert Transform Applications in Mechanical Vibration https://doi.org/10.1002/9781119991656 (2011).
 24.
Smith, Z. M., Delgutte, B. & Oxenham, A. J. Chimaeric sounds reveal dichotomies in auditory perception. Nature https://doi.org/10.1038/416087a (2002).
 25.
Huang, N. E. & Wu, Z. A review on HilbertHuang transform: method and its applications to geophysical studies. Rev. Geophys. https://doi.org/10.1029/2007RG000228 (2008).
 26.
Bouthour, W. et al. Biomarkers for closedloop deep brain stimulation in Parkinson disease and beyond. Nat. Rev. Neurol. https://doi.org/10.1038/s4158201901664 (2019).
 27.
Miterko, L. N. et al. Consensus Paper: experimental Neurostimulation of the Cerebellum. Cerebellum https://doi.org/10.1007/s12311019010415 (2019).
 28.
Cagnan, H. et al. Stimulating at the right time: phasespecific deep brain stimulation. Brain https://doi.org/10.1093/brain/aww286 (2017).
 29.
Shih, L. C. & PascualLeone, A. Noninvasive brain stimulation for essential tremor. Tremor Other Hyperkinetic Mov. https://doi.org/10.7916/D8G44W01 (2017).
 30.
Maas, R. P. P. W. M., Helmich, R. C. G. & van de Warrenburg, B. P. C. The role of the cerebellum in degenerative ataxias and essential tremor: insights from noninvasive modulation of cerebellar activity. Mov. Disord. https://doi.org/10.1002/mds.27919 (2020).
 31.
Cagnan, H. et al. The nature of tremor circuits in parkinsonian and essential tremor. Brain https://doi.org/10.1093/brain/awu250 (2014).
 32.
Brittain, J. S., ProbertSmith, P., Aziz, T. Z. & Brown, P. Tremor suppression by rhythmic transcranial current stimulation. Curr. Biol. https://doi.org/10.1016/j.cub.2013.01.068 (2013).
 33.
Nieuwhof, F., Panyakaew, P., Van De Warrenburg, B. P., Gallea, C. & Helmich, R. C. The patchy tremor landscape: recent advances in pathophysiology. Curr Opinion Neurol. https://doi.org/10.1097/WCO.0000000000000582 (2018).
 34.
Chiken, S. & Nambu, A. Mechanism of Deep Brain Stimulation: inhibition, Excitation, or Disruption? Neuroscientist https://doi.org/10.1177/1073858415581986 (2016).
 35.
Tse, N. Y. et al. The effect of stimulation interval on plasticity following repeated blocks of intermittent theta burst stimulation. Sci. Rep. https://doi.org/10.1038/s4159801826791w (2018).
 36.
Giordano, J. et al. Mechanisms and effects of transcranial direct current stimulation. DoseResponse https://doi.org/10.1177/1559325816685467 (2017).
 37.
Iwata, N. K. et al. Facilitatory effect on the motor cortex by electrical stimulation over the cerebellum in humans. Exp. Brain Res. https://doi.org/10.1007/s002210041979x (2004).
 38.
Cooper, I. S. & Upton, A. R. M. Use of chronic cerebellar stimulation for disorders of disinhibition. Lancet https://doi.org/10.1016/S01406736(78)910383 (1978).
 39.
Asamoah, B., Khatoun, A. & Mc Laughlin, M. tACS motor system effects can be caused by transcutaneous stimulation of peripheral nerves. Nat. Commun. https://doi.org/10.1038/s4146701808183w (2019).
 40.
Chakraborty, S. et al. Intermittent bilateral coherence in physiological and essential hand tremor. Clin. Neurophysiol. https://doi.org/10.1016/j.clinph.2016.12.027 (2017).
 41.
Kramer, G., Van der Stouwe, A. M. M., Maurits, N. M., Tijssen, M. A. J. & Elting, J. W. J. Wavelet coherence analysis: a new approach to distinguish organic and functional tremor types. Clin. Neurophysiol. https://doi.org/10.1016/j.clinph.2017.10.002 (2018).
 42.
Cao, C. et al. Corticosubthalamic coherence in a patient with dystonia induced by choreaacanthocytosis: a case report. Front. Hum. Neurosci. https://doi.org/10.3389/fnhum.2019.00163 (2019).
 43.
Di Biase, L. et al. Tremor stability index: a new tool for differential diagnosis in tremor syndromes. Brain https://doi.org/10.1093/brain/awx104 (2017).
 44.
Louis, E. D. Treatment of essential tremor: are there issues we are overlooking? Front. Neurol. 2, p. 91 JAN, (2012).
 45.
Fasano, A. & Deuschl, G. Therapeutic advances in tremor. Mov. Disord. https://doi.org/10.1002/mds.26383 (2015).
 46.
Deuschl, G. New hope for severe essential tremor? Lancet Neurol. https://doi.org/10.1016/S14744422(13)700620 (2013).
 47.
Hewitt, E. & Hewitt, R. E. The GibbsWilbraham phenomenon: an episode in fourier analysis. Arch. Hist. Exact Sci. https://doi.org/10.1007/BF00330404 (1979).
 48.
Deuschl, G., Bain, P. & Brin, M. Consensus statement of the Movement Disorder Society on Tremor. Ad Hoc Scientific Committee. Mov. Disord. 13(S3), pp. 2–23. (1998).
 49.
McIntyre, C. C., Grill, W. M., Sherman, D. L. & Thakor, N. V. Cellular Effects of Deep Brain Stimulation: modelBased Analysis of Activation and Inhibition. J. Neurophysiol. https://doi.org/10.1152/jn.00989.2003 (2004).
 50.
van Albada, S. J. & Robinson, P. A. Meanfield modeling of the basal gangliathalamocortical system. I. Firing rates in healthy and parkinsonian states. J. Theor. Biol. https://doi.org/10.1016/j.jtbi.2008.12.018 (2009).
 51.
Dayan, P. & Abbott, L. F. Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems. (MIT Press, Cambridge, Massachusetts 2001).
 52.
Santaniello, S. et al. Therapeutic mechanisms of highfrequency stimulation in parkinson’s disease and neural restoration via loopbased reinforcement. Proc. Natl Acad. Sci. U.S.A. https://doi.org/10.1073/pnas.1406549111 (2015).
 53.
McCarthy, M. M. et al. Striatal origin of the pathologic beta oscillations in Parkinson’s disease. Proc. Natl Acad. Sci. U.S.A. https://doi.org/10.1073/pnas.1107748108 (2011).
 54.
Ardid, S. et al. Biased competition in the absence of input bias revealed through corticostriatal computation. Proc. Natl Acad. Sci. U. S. A. https://doi.org/10.1073/pnas.1812535116 (2019).
 55.
Molnar, G. F., Pilliar, A., Lozano, A. M. & Dostrovsky, J. O. Differences in neuronal firing rates in pallidal and cerebellar receiving areas of thalamus in patients with Parkinson’s disease, essential tremor, and pain. J. Neurophysiol. https://doi.org/10.1152/jn.00881.2004 (2005).
 56.
Sedov, A. et al. Pallidal Activity in Cervical Dystonia with and Without Head Tremor. Cerebellum https://doi.org/10.1007/s12311020011195 (2020).
 57.
Gaffield, M. A., Bonnan, A. & Christie, J. M. Conversion of Graded Presynaptic Climbing Fiber Activity into Graded Postsynaptic Ca2+ Signals by Purkinje Cell Dendrites. Neuron https://doi.org/10.1016/j.neuron.2019.03.010 (2019).
 58.
Masoli, S., Solinas, S. & D’Angelo, E. Action potential processing in a detailed Purkinje cell model reveals a critical role for axonal compartmentalization. Front. Cell. Neurosci. https://doi.org/10.3389/fncel.2015.00047 (2015).
 59.
Diwakar, S., Magistretti, J., Goldfarb, M., Naldi, G. & D’Angelo, E. Axonal Na+ channels ensure fast spike activation and backpropagation in cerebellar granule cells. J. Neurophysiol. https://doi.org/10.1152/jn.90382.2008 (2009).
 60.
Perkel, D. J., Hestrin, S., Sah, P. & Nicoll, R. A. Excitatory synaptic currents in Purkinje cells. Proc. R. Soc. B Biol. Sci. https://doi.org/10.1098/rspb.1990.0074 (1990).
 61.
Lachaux, J. P., Rodriguez, E., Martinerie, J. & Varela, F. J. Measuring phase synchrony in brain signals. Hum. Brain Mapp. (1999). 10.1002/(SICI)10970193(1999)8:4<194::AIDHBM4>3.0.CO;2C
 62.
Iacono, M. I. et al. MIDA: a multimodal imagingbased detailed anatomical model of the human head and neck. PLoS ONE https://doi.org/10.1371/journal.pone.0124126 (2015).
 63.
Hasgall, P. A. et al. “IT’IS Database for thermal and electromagnetic parameters of biological tissues,” Version 4.0, May 15, 2018, https://doi.org/10.13099/VIP21000040. itis.swiss/database.
Acknowledgements
S.R.S. was supported by the Swiss National Science Foundation, Swiss Neurological Society and European Academy of Neurology and EMDO Foundation fellowship. X.Z. was supported in part by the CT Institute for the Brain and Cognitive Sciences IBRAiN Fellowship. S.S. was supported in part by the US NSF CAREER Award 1845348. K.P.B. was supported by Wellcome Trust MRC strategic neurodegenerative disease initiative award (WT089698), the Dystonia Coalition, Parkinson’s UK (G1009). M.B. was supported by EPSRC award EP/N014529/1 funding the EPSRC Centre for Mathematics of Precision Healthcare. E.S.B. acknowledges Lisa Yang, John Doerr, NIH R01MH117063, Edward and Kay Poitras, HHMI, and DARPA D20AC00004. N.G. was funded by the UK Dementia Research Institute (UK DRI)—an initiative funded by the Medical Research Council, Alzheimer’s Society and Alzheimer’s Research UK, Wellcome Trust fellowship (097443/Z/ 11/Z), Science & PINS Award for Neuromodulation, and NIHR IBRC Confident in Concept Award. The authors would like to thank Elisabeth Rounis and Tom Foltynie for helping with participant recruitment.
Author information
Affiliations
Contributions
S.R.S., designed and conducted clinical study, oversaw phaselocking, tremor amplitude and featurebased statistical learning analyses, wrote the paper. D.W., developed ecHT, design and implemented ecHTbased phaselocking brain stimulator. R.L.P., designed, developed and conducted featurebased statistical learning analysis. J.L., developed and conducted phaselocking and tremor amplitude analyses. E.R. and A.L., conducted the repeated clinical study. E.P., helped developing ecHT theory. A.M.C., conducted the FEM simulation. E.S.B., developed ecHT, oversaw experiments. M.B., designed and oversaw featurebased statistical learning analysis. X.Z. and S.S., designed, developed and conducted neurophysiological simulation of ET. K.P.B. and J.R., designed and oversaw clinical study. N.G., developed ecHT, designed and conducted clinical study, designed, developed and oversaw phaselocking and tremor amplitude analyses, oversaw featurebased statistical learning analysis and neurophysiological simulation of ET, and wrote the paper.
Corresponding authors
Ethics declarations
Competing interests
N.G., D.W. and E.S.B. have applied for a patent on the ecHT technology, assigned to MIT, and founded a company that utilises it. K.P.B., received funding for travel from GlaxoSmithKline, Orion Corporation, Ipsen, and Merz Pharmaceuticals, LLC; serves on the editorial boards of Movement Disorders and Therapeutic Advances in Neurological Disorders; receives royalties from the publication of Oxford Specialist Handbook of Parkinson’s Disease and Other Movement Disorders (Oxford University Press, 2008); received speaker honoraria from GlaxoSmithKline, Ipsen, Merz Pharmaceuticals, LLC, and Sun Pharmaceutical Industries Ltd.; personal compensation for scientific advisory board for GSK and Boehringer Ingelheim; received research support from Ipsen and from the Halley Stewart Trust through Dystonia Society UK. The rest of the authors declare no competing interests.
Additional information
Peer review information Nature Communications thanks the anonymous reviewers for their contributions to the peer review of this work. Peer review reports are available
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Source data
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Schreglmann, S.R., Wang, D., Peach, R.L. et al. Noninvasive suppression of essential tremor via phaselocked disruption of its temporal coherence. Nat Commun 12, 363 (2021). https://doi.org/10.1038/s41467020205817
Received:
Accepted:
Published:
Further reading

Realtime estimation of phase and amplitude with application to neural data
Scientific Reports (2021)

Essential tremor
Nature Reviews Disease Primers (2021)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.