- 1Institute of Neurosurgery and Neurointervention, Faculty of Medicine, Semmelweis University, Budapest, Hungary
- 2Albert Szent-Györgyi Medical School, Doctoral School of Clinical Medicine, Clinical and Experimental Research for Reconstructive and Organ-Sparing Surgery, University of Szeged, Szeged, Hungary
- 3Department of Stereotactic and Functional Neurosurgery, Medical Center of Freiburg University and Medical Faculty of Freiburg University, Freiburg, Germany
- 4János Szentágothai Doctoral School of Neurosciences, Semmelweis University, Budapest, Hungary
- 5CereGate GmbH, München, Germany
- 6Department of Neurology, Faculty of Medicine, Semmelweis University, Budapest, Hungary
- 7Center for Deep Brain Stimulation, Freiburg University, Freiburg, Germany
Introduction: Although stimulation-induced sensations are typically considered undesirable side effects in clinical DBS therapy, there are emerging scenarios, such as computer-brain interface applications, where these sensations may be intentionally created. The selection of stimulation parameters, whether to avoid or induce sensations, is a challenging task due to the vast parameter space involved. This study aims to streamline DBS parameter selection by employing a machine learning model to predict the occurrence and somatic location of paresthesias in response to thalamic DBS.
Methods: We used a dataset comprising 3,359 paresthetic sensations collected from 18 thalamic DBS leads from 10 individuals in two clinical centers. For each stimulation, we modeled the Volume of Tissue Activation (VTA). We then used the stimulation parameters and the VTA information to train a machine learning model to predict the occurrence of sensations and their corresponding somatic areas.
Results: Our results show fair to substantial agreement with ground truth in predicting the presence and somatic location of DBS-evoked paresthesias, with Kappa values ranging from 0.31 to 0.72. We observed comparable performance in predicting the presence of paresthesias for both seen and unseen cases (Kappa 0.72 vs. 0.60). However, Kappa agreement for predicting specific somatic locations was significantly lower for unseen cases (0.53 vs. 0.31).
Conclusion: The results suggest that machine learning can potentially be used to optimize DBS parameter selection, leading to faster and more efficient postoperative management. Outcome predictions may be used to guide clinical DBS programming or tuning of DBS based computer-brain interfaces.
Introduction
Deep brain stimulation (DBS) has emerged as a promising treatment for a variety of neuropsychiatric diseases (1). It is established for movement disorders such as Parkinson’s disease (2–4), essential tremor (5, 6), and dystonia (7, 8), and for epilepsy (9, 10). It is an emerging therapy in neuropathic pain (11) and selected psychiatric disorders like major depression and obsessive-compulsive disorder (12, 13). Despite its clinical efficacy, DBS is frequently accompanied by undesirable side effects, which can adversely affect patient satisfaction and treatment outcomes (14–16). Conversely, in emerging technologies such as computer-brain interfaces (CBI), implanted electrodes can be used to deliberately induce sensations like paresthesias (17, 18).
With respect to clinical applications, major challenges in the post-implant management of DBS therapy remain. In particular, optimizing stimulation parameters presents a significant bottleneck in DBS therapy. One of the most important factors is searching the vast parameter space involved in DBS programming, which is often done manually (15, 19). Computer algorithms offer significant potential to improve the efficiency and accuracy of DBS parameter selection, whether the intended outcome is maximizing symptom reduction, minimizing side effects, or eliciting specific sensations for CBI applications. This work addresses the feasibility of using machine learning to improve the efficiency of DBS parameter selection, potentially overcoming some of these limitations.
Understanding the factors that contribute to paresthesias and devising strategies to mitigate side effects are essential to optimize DBS therapy. It is proposed that paresthesias result from the activation of sensory afferents near the DBS lead (20–22). Consequently, in DBS for tremor, for example, paresthesia induction can signal the activation of critical anatomical landmarks, such as the medial lemniscus and the posterior-medial border of the cerebellothalamic tract (21, 22).
One commonly used method to assess the effects of DBS is to simulate the Volume of Tissue Activation (VTA), which estimates the spatial extent of neural activation in response to stimulation (23–25). VTA models primarily rely on stimulation pulse width and current amplitude and may not fully capture the intricate relationships among all relevant stimulation parameters, individual neuroanatomy, tissue heterogeneity, and their effects on perceptual phenomena (26–28).
Here we aimed to predict the occurrence and spatial localization of paresthesias in response to thalamic DBS. To achieve this, we used a dataset comprising more than three thousand paresthesia location records, empirically obtained from 18 thalamic DBS leads in 10 individuals across two clinical centers. We modeled the VTA for each stimulation trial. Stimulation parameters and VTA metrics were used as input variables for the prediction models, which were evaluated using cross-validation.
Methods
Study design
This is a cross sectional study of a single cohort of 10 individuals (three females), from two independent deep brain stimulation centers. We recruited people who had chosen deep brain stimulation (DBS) surgery to treat either tremor or chronic neuropathic pain. Demographics and the anatomical target regions of the DBS leads (model DB2202, Boston Scientific Corporation, Marlborough, MA, United States) are summarized in Table 1. The implantations took place between July 2019 and October 2020. The experiments took place between 1 and 4 days after the surgery, during the period when externalized extension cables were connected to the DBS leads. We conducted the experiments in accordance with local guidelines and regulations and in accordance with the Declaration of Helsinki. The ethics committees of both centers approved the study (agreement number 235/19 for Freiburg; OGYEI/23818/2019 for Budapest). All participants provided written informed consent prior to the experimental procedures.
Brain imaging
Magnetic resonance imaging
Magnetic resonance imaging (MRI) was performed without a stereotactic frame. MR examinations comprised T1-weighted, T2-weighted sequences.
In Freiburg, MR imaging was performed 1–3 days before surgery, and - if necessary - with the person under mild sedation. MR images were acquired on a whole-body 3 T MR system (PRISMA, Siemens Healthcare, Erlangen, Germany) using a 64 channel phased array head coil. For T1-weighted imaging a magnetization-prepared rapid gradient-echo (MP-RAGE) scan was used (TR 2,300 ms, TE 2.26 ms, flip angle 12°, FOV 256 mm, voxel size 1mm3). The T2-weighted scan was a fast spin-echo sequence (TR 2,500 ms, TE 231 ms, FOV 256 mm, voxel size 1mm3).
In Budapest, MR imaging was performed between 1 month and 5 days before surgery, under general anesthesia. MR images were acquired on a 3 T whole-body magnetic resonance imaging system (Philips Healthcare, Best, The Netherlands) using an 8-element phased array head coil. For T1-weighted imaging a magnetization-prepared rapid gradient-echo (MP-RAGE) scan was used (TR 8.5 ms, TE 3.8 ms, flip angle 8°, FOV 256 mm, reconstructed to 1mm3 voxels). The T2-weighted scan was a fast spin-echo sequence (TR 12.65 ms, TE 100 ms, FOV 254 mm, reconstructed to 1.44mm3 voxels).
Computed tomography
Pre-operative and post-operative computed tomography scans (CT) were both performed on the day of surgery, with the exception of p10 whose post-operative CT was taken the day after surgery. Pre-operative CT was acquired to register the stereotactic frame to the person’s anatomy. Post-operative CT was acquired to determine the positions of the implanted DBS electrode arrays.
In Freiburg, CT scans were acquired with a SOMATOM Definition AS scanner (Siemens Healthcare, Erlangen, Germany). The parameters of the preoperative CT were as follows: tube voltage 120 kV, tube current 365 mAs, collimation 19·0.6 mm, tube rotation time 1.0 s, pitch 0.55, matrix 512·512, section thickness 1.5 mm, increment 1.5 mm. The parameters of the post-operative CT were identical, except for a tube current of 320 mAs. The different post-operative scanning parameters were chosen for better electrode metal artifact suppression.
In Budapest, CT scans were acquired with a Brilliance 8,000 16-row multidetector scanner (Philips Healthcare, Best, Netherlands). The parameters of the preoperative CT were as follows: tube voltage 120 kV, tube current 350 mAs, collimation 16·0.75 mm, tube rotation time 1.0 s, pitch 0.942, matrix 512·512, section thickness 1.5 mm, increment 1.5 mm. The parameters of the post-operative CT were identical, except for the tube rotation time of 0.75 s, pitch of 0.688, section thickness of 2 mm, and an increment of 1 mm. The different post-operative scanning parameters were chosen for better electrode metal artifact suppression.
Brain stimulation
We performed the experiments in dedicated laboratory rooms in each of the two clinics. The participant sat on a chair, behind a regular office desk. At a distance of about 60–80 cm from the participant, there was a PC monitor to provide feedback to the participant during the experiment. The two participants with two leads implanted on the same were only stimulated through their thalamic electrodes. We delivered the stimuli with an external programmable neurostimulator. In Freiburg, we used a Neuro Omega stimulator (Alpha Omega Engineering, Nof HaGalil, Israel). In Budapest we used a CereStim stimulator (BlackRock, Salt Lake City, UT, USA). We connected the neurostimulator and the DBS lead using a surgical operating room cable (SC-4108; Boston Scientific, Marlborough, MA, USA). We used custom software written in Matlab (MathWorks, Natick, MA) to send instructions to the external neurostimulator interface using the manufacturer’s software development kit. We asked participants to provide feedback about stimulation-induced sensations, either orally or via a button box (ResponsePixx; VPixx Technologies, Saint-Bruno, QC, Canada), depending on their preference. The range of stimulation parameters was restricted to the established safety ranges for deep brain stimulation. Stimulation parameters included the following ranges: pulse frequency 20–205 Hz, pulse width 50–250 μs, and pulse amplitude 0.25–3.0 mA. With the CereStim stimulator, electrode contact number 8 was used for the return current path. With Neuro Omega, contact 8 was used as a passive ground. To keep the applied waveforms comparable between the different stimulators, we generated a basic waveform that mimics monophasic passive discharge. This waveform was the default throughout the experiments.
Data analysis
Image processing
Image processing was performed on a “NeuroImaging Tools and Resources Collaboratory Cloud Computing Environment” (NITRC-CE) (29), running on an Amazon Web Services EC2 P3 virtual machine equipped with NVIDIA Tesla V100 GPUs. DICOM images were converted to compressed NIFTI format using MRICron’s dcm2niix converter (30). Tools available in the FMRIB Software Library (FSL 6.0.6.2) were used for linear and non-linear image registrations.
Modeling of volume of tissue activated
The Volume of Tissue Activated (VTA) was modeled using Lead-DBS 2.6 (31) software with the FastField (32) modeling method, running in Matlab (version 2023a). For each stimulation trial, we used the respective electrode configuration, stimulation currents, and cathodic pulse width as input parameters for the VTA calculation. The axon diameter parameter was fixed to 3.5 μm. DBS contact positions and orientation were determined from the post-operative CT images using the DiODe algorithm from Lead-DBS, and subsequently visually verified using Brainlab Elements (Brainlab AG, Munich, Germany). To allow pooling of cases, the DBS contact positions were transformed from individual space into MNI space. The transformation matrix was determined by non-linear registration (12 degrees of freedom) of the person’s T1 MRI to the 0.5mm3 MNI ICBM 2009b NLIN ASYM brain template (MNI hereafter) using ANTs (33). Subsequently, the VTA modeling was done in MNI space. For further analysis, VTAs were converted into binary masks.
Data augmentation and preparation
During the experiment, in the majority of trials stimulation was applied above the sensation threshold. Thus, the empirical dataset contained relatively few instances in which no paresthesias were reported. However, our aim extended beyond predicting the somatic area of evoked sensations; we also aimed to predict the presence or absence of sensations altogether. To address this, we augmented the empirical dataset by introducing pseudotrials that were synthesized from existing empirical data. Pseudotrials were generated from the subset of stimulation trials where no paresthesias were reported. This synthesis involved systematically generating novel combinations of pulse width, frequency, and stimulation currents below the empirically observed sensation thresholds.
The empirical dataset included multiple identical stimulation and response instances (e.g., due to repeating stimulation with identical parameters). We retained only those trials with distinct combinations of stimulation parameters and responses. This also ensured that during cross-validation, the test splits exclusively comprised data patterns that had not been encountered at all during the training phase. Subsequently, a random subset (sampled without replacement) of the pseudotrials were added to the set of distinct empirical trials to balance out the dataset such that for each individual DBS lead equal proportions of trials with and without sensations were available. The augmented dataset, comprising both empirical trials and synthetically generated pseudotrials, was used for further processing, including VTA generation.
Spatial data that are related to right-sided DBS electrodes were flipped so that all data is virtually projected onto the left hemisphere (34, 35). Consequently, analyses were performed in a unilateral frame of reference. For reported paresthesias, contralateral denotes the hemibody opposite to the stimulation location.
Definition and cross-validation of prediction models
For the predictions, a series of binary classifiers were trained. One binary classifier was trained to predict whether any sensation would occur or not (i.e., without predicting in which specific body parts). For this classifier we only provided stimulation-energy related parameters as features (i.e., Pw, Freq, Current, Hemisphere; see Table 2). Since the purpose of this classifier is to predict whether or not any sensation occurs (without concern for the specific location in the body), we limited the features to energy-related parameters. These features are expected to correlate with whether a sensation is triggered but should not reveal where the sensation occurs. Preserving spatial information for the specific somatic predictions also avoided redundancy. To predict the occurrence of sensations in specific body parts, a separate binary classifier was trained for each somatic area (e.g., finger, wrist, etc.). To these somatic prediction models we provided additional features about the spatial location of the stimulation energy (i.e., VTAspread, VTAx, VTAy, VTAz, Freq, Hemisphere; see Table 2).
The prediction models used the LogitBoost algorithm for classification and were trained and evaluated using a nested two-level (K by L) cross-validation approach (illustrated in Figure 1). In this procedure, the outer K loop iterates over the data, splitting it into training-and-validation and test sets. We used two different approaches for outer K-fold splitting, which is further explained in the next section “Evaluation of predictions.” The inner L loop further splits the K-th training set into multiple folds. Prediction models were trained using cross-validation on the inner L folds. The splits in the inner loop were based on group folding by lead-id, ensuring that all trials from the same lead were held out for validation, thus enhancing generalization. Within each L iteration, the training-validation subset for each somatic prediction model was rebalanced so that both classes (True/False) had an equal number of occurrences. This was achieved by identifying the least frequent class and randomly dropping trials from the more frequent class. Trials where no sensations occurred, and thus no somatic areas were reported, were excluded in the training of the somatic prediction models. Only somatic prediction models with a Kappa value above 0.05 after L-fold cross-validation were retained. These models were then evaluated using the test sets from the outer K fold.
Figure 1. Schematic of the nested cross-validation procedure. The outer K loop iterates over the data, splitting it into training-and-validation and test sets. In each iteration, a portion of the data is held out as a test set (outer fold), while the remaining data is used for training and validation. Within each K fold, the inner L loop further splits the training set into multiple folds to perform cross-validation. The prediction models are trained on the training folds from the inner loop and evaluated on the validation folds. After completing the inner loop, resulting models are then evaluated using the test set from the outer K fold.
Model training and evaluation analyses were carried out in R (version 4.4).
Evaluation of predictions
We evaluated the predictions in two scenarios, which differed in the way the training-and-validation and test splits were created:
• Intra-sample evaluation: The K test sets consisted of a random sample of trials with parameter-response combinations distinct from all trials in the training and validation sets. However, these hold-out trials were drawn from cases (i.e., leads) that were also included in the L-loop training and validation cycles. This means the test set includes trial instances the model has not explicitly been trained on, but these still come from the set of cases that were part of the training/validation process. Intra-sample performance reflects the model’s ability to generalize and predict outcomes for new stimulation parameter combinations within the same individuals used for training (i.e., like patients for whom prior responses are on record).
• Inter-sample (leave-cases-out) evaluation: The K test sets only included trials from cases (i.e., leads) that were entirely excluded from the L-loop iterations of the cross-validation process. This means the test set contains data from cases that were completely unseen during the L-loop training and validation. These trials represent truly unseen data from new cases (i.e., similar to de novo patients), providing a more challenging evaluation of the model’s generalization ability.
In both scenarios, resampling with replacement was used to ensure a uniform number of trials per case in the test set. The number of trials per case was set to the average trial count across cases. This approach ensured that each case contributed equally to the aggregated test performance metrics, while maintaining the original overall size of the test set. The model’s performance was evaluated on two aspects: (1) Prediction of the occurrence of paresthesia (“paresthesia predictions”), and (2) prediction of the somatic locations where a paresthesia would occur (“somatic predictions”). Since stimulation could elicit paresthesias in none, one, or multiple body regions simultaneously, predicting the somatic locations is a multi-output task, where multiple target variables are predicted in parallel. Thus, for each trial, the paresthesia prediction and the multiple somatic predictions were read out in parallel and combined. The somatic predictions in a given trial were retained only if the paresthesia prediction was positive. Otherwise, all somatic area predictions were overridden to False.
We used Kappa as the primary metric to compare model performance. Kappa is a robust measure of agreement that accounts for chance. In addition, we used accuracy, sensitivity, specificity, precision, F1-score. For the somatic prediction metrics the prediction of each of the multiple binary target variables was taken into account and aggregated. Our “distance” metric summarizes the average difference between observed and predicted composite somatic response (“thumb + finger + wrist”). A distance of 0 indicates an exact agreement between predictions and observations, whereas a distance of 0.5 signifies chance-level agreement, implying predictions are no better than random guessing. We report descriptive statistics as averages followed by the associated standard deviation, unless stated otherwise. The model performance metrics were aggregated across the K iterations. For statistical comparison of model metrics between scenarios, we used the non-parametric Wilcoxon rank sum test.
Control analyses
To evaluate the contribution of VTA information to somatic predictions, we conducted a control analysis. In this analysis, we repeated the nested cross-validation and evaluation procedure described earlier, but excluded the spatial features. Instead, we provided the model with Pw, Freq, Current, Hemisphere—the same features used for the paresthesia occurrence prediction model.
Additionally, we assessed the performance of two types of naive prediction models and used their metrics as baselines for comparison against actual predictions. The first naive model generated predictions by randomly permuting the reference (ground truth) response data. Consequently, its performance was influenced by the underlying class distribution. The second naive model always predicted True for both paresthesia occurrence and somatic area prediction.
Results
A three-dimensional overview of all 18 DBS leads, rendered together with anatomical landmarks in a common space is provided in Figure 2. A triple-plane orthographic rendering of all the VTAs in the dataset is shown in Figure 3.
Figure 2. Graphic rendering of the 18 DBS leads with anatomical landmark structures, in a common space. The anatomical structures are sourced from the DISTAL atlas (54). VIM is colored in yellow; VPL/VPM in blue; Medial Lemniscus in light-yellow. The background is an axial T1-MRI slice (z = 0 mm) of the MNI template. The mean coordinates (±SD) of the deepest contact points of the VIM-targeted leads were x = −12.6 (±2.1), y = −19.8 (±1.8), z = −3.9 (±1.8) in the left hemispheres, and x = 13.6 (±1.6), y = −19.6 (±1.7), z = −3.5 (±1.9) in the right hemispheres. For the VPL-targeted leads these were x = 10.1 (±0.4), y = −24.2 (±1.7), z = −7.8 (±0.9).
Figure 3. VTAs rendered and overlaid on a MNI T1 Magnetic Resonance Image (MRI) in axial (z = 141), coronal (y = 232), and sagittal (x = 170) slices. Right-sided VTAs were mirrored so that these are projected onto the left hemisphere. The crosshairs indicate where the slices intersect. The yellow rectangles delineate the target areas of each magnification inset. VTA pixels are color coded by lead number and flattened along the axis that is not visible in that specific view. VTA pixels are transparent to allow overlapping areas to blend. A, anterior; L, left; S, superior; P, posterior.
The empirical data from the experiment that were used as the basis for the training and evaluation are summarized in Table 3. Across the empirical dataset, the ratio between trials with paresthetic sensations vs. without was 87 vs. 13%. Figure 4 shows an overview of the somatic areas in which paresthesias were reported. Each lead evoked paresthesias in at least 2 different somatic areas, with the exception of lead 17 (participant 10) for which paresthesias were evoked exclusively in the fingers. On average (±SD), a DBS lead evoked paresthesias in 8.4 ± 5.1 different somatic areas. Paresthesias were most frequently reported in the fingers, palm (of the hand), and the thumb. For training and evaluation of the prediction model we retained the subset of empirical trials with distinct combinations of feature and response values. In total, 667 pseudotrials were added to balance trials with and without paresthesias. This resulted in a total of 1,451 unique trials, with a 50 − 50% balance between no-paresthesia and paresthesia trials. In total, the prediction dataset comprised 3,359 reported paresthesias across 21 different body regions. The somatic areas from the trials in the prediction dataset, and their counts, are summarized in Figure 5 and in Supplementary Figure S1.
Table 3. Counts and percentages of stimulation trials and their responses obtained in the experiment.
Figure 4. Body maps showing the locations of reported paresthesias for each of the 18 DBS electrode arrays. Left- and right-sided sensations are collapsed into a unilateral frame of reference, using a color code to indicate contra- vs. ipsilateral sensations. Percentages reflect the proportion of all reported sensations in each body area, calculated per DBS lead.
Figure 5. Total number of paresthesias recorded for each somatic area in the dataset used for training and testing the prediction models. The darker shaded bars represent somatic areas for which specific prediction models were trained and tested. For the remaining somatic areas (light shaded bars), no direct prediction models were trained. However, these areas were used as negative instances in other somatic area models, representing trials where the target sensation of the respective model did not occur (e.g., chin paresthesias could act as False instances for a model predicting paresthesias in the fingers).
The allocation of cases across K and L splits in the inter-sample (leave-case-out) cross-validation scenario is shown in Supplementary Figure S2. There were 8 K splits, each containing 5 L splits. The fold creation resulted in an average split ratio of 87.5% (training) vs. 12.5% (testing) for the K-folds, and 80% (training) vs. 20% (validation) within the L-folds.
The predictions and associated performance metrics of the model in the two different scenarios are summarized in Figure 6. Paresthesia prediction performance was comparable between the intra-sample and the inter-sample scenarios. Although the performance metrics suggested worse performance in the latter scenario, the difference in Kappa was not significant (0.72 ± 0.11 vs. 0.60 ± 0.26; p = 0.22).
Figure 6. Confusion matrices and performance metrics from the nested cross-validation of prediction models. Each cell shows the total count of true/false positives/negatives, with the corresponding percentage in parentheses. Percentages represent the proportion of total predictions classified as True (paresthesia occurring) or False (paresthesia not occurring), with model predictions on the vertical axis and the reference (ground truth) on the horizontal axis. The panels depict: (A) Intra-sample paresthesia predictions, (B) Intra-sample somatic predictions, (C) Inter-sample paresthesia predictions, and (D) Inter-sample somatic predictions. The tables beneath each matrix show the corresponding performance metrics.
For the somatic predictions, the difference between the two scenarios was more pronounced. Kappa was significantly lower in the inter-sample scenario than in the intra-sample scenario (0.53 ± 0.09 vs. 0.31 ± 0.10; p = 0.002). There were also significant differences in accuracy (83 ± 3% vs. 75 ± 7%; p = 0.005) and F1 scores (0.64 ± 0.08 vs. 0.46 ± 0.08; p = 0.002). The differences were mainly attributed to a significant reduction in sensitivity, from 67 ± 11% in the intra-sample scenario to 43 ± 12% in the inter-sample scenario (p = 0.002), whereas we did not find a significant change in specificity (88 ± 4% vs. 87 ± 10%; p = 0.64). This performance difference was further highlighted by a significant shift in the composite somatic response distance metric. The distance between predicted and reference (ground truth) response was 0.17 ± 0.03 in the intra-sample scenario, whereas it was 0.25 ± 0.07 in the inter-sample scenario (p = 0.005). The difference in precision was not significant (63 ± 10% vs. 55 ± 15%; p = 0.16).
To assess whether the VTA information had a relevant contribution to the somatic predictions we performed a control analysis. In this analysis, we trained and tested the predictions using a feature set that excluded the VTA-related features, leaving only Pw, Freq, Current, Hemisphere. For the intra-sample scenario, the removal of the VTA-related information led to a significant reduction in Kappa to 0.43 ± 0.06 (p = 0.016). Sensitivity dropped from 67 ± 11% to 57 ± 6% (p = 0.033), and the F1-score decreased from 0.64 ± 0.08 to 0.56 ± 0.05 (p = 0.016). Accuracy dropped from 83 ± 3% to 80 ± 3%, while the composite response distance increased from 0.17 ± 0.03 to 0.20 ± 0.03, with both changes just above the significance threshold (both p = 0.064). No significant decreases were observed in specificity (88 ± 4% vs. 87 ± 5%; p = 0.35) and precision (63 ± 10% vs. 56 ± 10%; p = 0.09). These results indicate that in the intra-sample scenario, the inclusion of spatial VTA information improved somatic predictions. In contrast, the control analysis showed less differences in the inter-sample scenario. Here the reduction in Kappa to 0.30 ± 0.08 was insignificant (p = 0.32). However, the reductions in precision (55 ± 16% to 47 ± 8%; p = 0.052) and specificity (87 ± 10% to 84 ± 3%; p = 0.052) were both close to the significant threshold. No significant changes were observed in the other metrics.
Finally, we compared the prediction performance metrics to the naive baseline predictions, as summarized in Supplementary Figure S3. For both scenarios, Kappa values for the actual predictions were significantly higher than those of the random and constant-baseline models (all p < 0.001, with one exception). In the inter-sample somatic predictions, the associated p-value (p = 0.026) was substantially larger than the p-value of the other control comparisons, indicating a more modest level of significance. In line with this finding, the naive predictions revealed a bias caused by a higher proportion (77%) of False instances in the ground truth for the somatic predictions. This was reflected in the accuracy metric of the naive models: Random predictions yielded an accuracy of 73%, while the constant True output (representing the opposite class) resulted in an accuracy of 25%. In light of these observations, the somatic prediction model appears to have marginal predictive power for unseen cases.
Discussion
In this study, we aimed to predict paresthesias’ occurrence and locations in response to thalamic DBS in individuals with two different clinical symptoms (tremor vs. neuropathic pain), and electrode locations (Vim-nucleus vs. VPL/VPM nuclei, respectively). Using a dataset comprising a large number of stimulation responses obtained empirically from 18 DBS lead implanted in 10 individuals, we explored how stimulation parameters and person specific features influenced perceptual outcomes. Our goal was to evaluate the predictive model’s performance to then in the future assess its real-world relevance for optimizing DBS therapy management and potential applications in computer-brain interfacing.
Our findings reveal several important insights: Firstly, the predictions demonstrated substantial agreement with ground truth for the occurrence of paresthesias, indicating the feasibility of using computational approaches to anticipate sensory outcomes in DBS. Secondly, the importance of incorporating individual-specific information from VTA simulations in predicting the somatic regions of DBS-induced paresthesias highlights the value of image guided DBS programming, underscoring the need to integrate such data in DBS programming and optimization workflow.
The distinction between intra-sample and inter-sample scenarios offers valuable perspective on the real-world application and generalizability of our results. We achieved comparable prediction performance for both previously encountered and unseen cases for the presence of paresthesias. However, the performance of predicting the somatic location of paresthesias decreased for unseen cases. Thus, the model appears to have captured individual characteristics better than generalizable patterns across different individuals. The control analyses with naive models revealed that the model’s somatic predictions for entirely unseen cases were only slightly better than random guessing, indicating that, at this point, the clinical relevance of these predictions may be limited. Somatic predictions seem to require individual-specific prior data, such as stimulation responses collected during prior programming sessions or even during surgery.
We extend earlier research that demonstrated the potential of machine learning models to predict motor response after deep brain stimulation in Parkinson’s disease (36) and Essential Tremor (37). Furthermore, a recent study on people with isolated dystonia demonstrated that machine-learning based DBS programming resulted in greater symptom reduction compared to clinical programming (38). To our knowledge, no prior studies have explored the use of machine-learning approaches for predicting the occurrence and location of paresthesias during DBS. Our findings support the role of sensory afferent fiber activation in mediating paresthesias during DBS (21, 22). However, we build further upon this understanding by demonstrating the predictive utility of computational models incorporating stimulation parameters and VTA metrics. This represents a significant advancement in the field by providing clinicians with a potential tool to anticipate and mitigate sensory side effects associated with DBS therapy.
Our findings could have implications for future clinical practice. Paresthesias can serve as an invaluable real-time feedback signal about the functional effects of DBS. Paresthesias can be a marker of efficacy in spinal cord stimulation (SCS) therapy, correlating with pain relief (39–41). Their role in DBS is less clear, but some studies suggest a similar relationship (42, 43). Furthermore, the location of the paresthesias can offer important information to determine the optimal stimulation site. Specifically, face or finger paresthesias suggests optimal placement within the Vim nucleus for tremor suppression. Conversely, deviations like intra-oral or leg paresthesias could indicate a need for stimulation adjustment, or electrode repositioning (21). Our paradigm can help clinicians predict paresthesias caused by DBS, allowing tailoring of treatment per individual. This has the potential to enhance therapy outcomes, reduce treatment-related morbidity, and improve overall satisfaction with DBS therapy. Moreover, our research supports the integration of computational modeling approaches into clinical practice to augment decision-making and improve patient care in neurostimulation therapies (44).
Our findings extend beyond clinical DBS and are relevant to the field of CBI, particularly in applications requiring precise sensory induction (17, 18). Our prediction paradigm could be leveraged in CBI systems that require modulation of neural activity to evoke specific sensations or perceptions. For instance, by adjusting stimulation parameters and targeting brain regions identified through paresthesia mapping, researchers could design CBI devices capable of generating controlled tactile sensations or proprioceptive feedback. These signals could be used to convey information from the external world, such as artificial touch perceived in a prosthetic limb or a phantom limb (17, 45, 46). Moreover, there is a rich body of work on sensory substitution showing that artificial sensory input can be used to restore, replace, or enhance sensory function (47–51). We have recently demonstrated that paresthesias evoked by spinal cord stimulation can be used to effectively convey a diverse range of information, from rhythmic cues to an artificial sense of balance (18). These advancements could have transformative implications for various CBI domains, including neuroprosthetics, virtual reality, and sensory augmentation technologies. By enabling enhanced sensory experiences and functional restoration for individuals with sensory impairments, our study paves the way for unlocking the full potential of DBS-based CBI, combining valuable insights from clinical therapy, paresthesia mapping, and sensory substitution.
Our study has several limitations. The study’s sample lacks a broad representation of individual characteristics and response variability. On the other hand, the sample’s heterogeneity may introduce confounding factors that could affect the generalizability of the findings. While the study showed promising results for this specific sample of individuals and experimental setting, it’s unclear how well these findings would translate to everyday clinical practice or different subpopulations. Factors such as variations in DBS electrode type and location, stimulation protocols, demographics, and medical condition could impact the model’s performance in real-world applications. Stimulation parameters and anatomical VTA information only partially account for the variation in responses to DBS (1, 52, 53). Unaccounted variables, such as medication or comorbidities, may have influenced the perception of paresthesias. However, we believe these factors are unlikely to introduce significant confounds given our primary objective: predicting stimulation effects using energy-related parameters and anatomical locations. Furthermore, with respect to a potential CBI approach our data does not cover long term paresthesia effects, since it was not the purpose of the paresthesia screening done here. Lastly, while this study provides insights into the predictive modeling of DBS-induced sensations, practical challenges related to model implementation and integration into clinical practice remain.
Looking ahead, several avenues for future research emerge. For example, further refinement and validation of predictive models incorporating larger cohorts, additional features, and biomarkers could enhance the accuracy and clinical utility of paresthesia prediction in DBS therapy. Also, investigations into the broader applicability of predictive modeling approaches across different DBS targets, outcomes, and populations would be valuable for advancing the field of neuromodulation.
Conclusion
We have demonstrated the feasibility of using machine learning to predict paresthesias in response to thalamic DBS. This predictive modeling provides clinicians with a powerful tool to optimize DBS programming, allowing them to either avoid or intentionally evoke paresthesias. Additionally, paresthesia mapping could enhance computer-brain interfaces, enabling precise sensory induction and offering transformative applications in neuroprosthetics, virtual reality, and sensory augmentation technologies. Our findings are intended to serve as a stepping stone for future, more extensive studies, such as those with larger cohorts. These advancements may pave the way for neuromodulation interventions that are better tailored to individual needs.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Ethics statement
The studies involving humans were approved by Freiburg University ethics committee (agreement 235/19); OGYÉI Budapest (agreement 23818/2019). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
LH: Conceptualization, Data curation, Investigation, Resources, Writing – original draft, Writing – review & editing. BS: Writing – review & editing, Data curation, Investigation, Resources. GM: Data curation, Writing – review & editing, Formal analysis, Software, Visualization. GE: Formal analysis, Software, Visualization, Writing – review & editing, Writing – original draft. SH: Data curation, Investigation, Software, Writing – review & editing. BV: Conceptualization, Data curation, Investigation, Supervision, Writing – review & editing. GT: Data curation, Resources, Writing – review & editing. VC: Conceptualization, Resources, Supervision, Writing – review & editing. LE: Conceptualization, Resources, Supervision, Writing – review & editing.
Funding
The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.
Acknowledgments
We would like to thank Till Dembek from the University Clinic of Cologne for sharing his scripts for batch VTA generation with Lead-DBS.
Conflict of interest
LH, VAC, and LE are compensated members of the CereGate GmbH Advisory Board. GM, GVE, SHG, and BV are employees of CereGate GmbH.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fneur.2024.1467307/full#supplementary-material
Abbreviations
CBI, Computer-brain interface; CT, Computed tomography; DBS, Deep brain stimulation; MRI, Magnetic resonance imaging; PAG, Periaqueductal gray; PVG, Periventricular gray; SCS, Spinal cord stimulation; VIM, Ventral intermediate thalamic nucleus; VPL, Ventroposterolateral thalamic nucleus; VPM, Ventroposteromedial thalamic nucleus; VTA, Volume of tissue activated.
References
1. Lozano, AM, Lipsman, N, Bergman, H, Brown, P, Chabardes, S, Chang, JW, et al. Deep brain stimulation: current challenges and future directions. Nat Rev Neurol. (2019) 15:148–60. doi: 10.1038/s41582-018-0128-2
2. Benabid, AL, Pollak, P, Louveau, A, Henry, S, and de Rougemont, J. Combined (thalamotomy and stimulation) stereotactic surgery of the VIM thalamic nucleus for bilateral parkinson disease. Appl Neurophysiol. (1987) 50:344–6. doi: 10.1159/000100803
3. Okun, MS. Deep-brain stimulation for parkinson’s disease. N Engl J Med. (2012) 367:1529–38. doi: 10.1056/NEJMct1208070
4. Limousin, P, and Foltynie, T. Long-term outcomes of deep brain stimulation in parkinson disease. Nat Rev Neurol. (2019) 15:234–42. doi: 10.1038/s41582-019-0145-9
5. Benabid, AL, Pollak, P, Gervason, C, Hoffmann, D, Gao, DM, Hommel, M, et al. Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus. Lancet. (1991) 337:403–6. doi: 10.1016/0140-6736(91)91175-T
6. Hopfner, F, and Deuschl, G. Managing essential tremor. Neurotherapeutics. (2020) 17:1603–21. doi: 10.1007/s13311-020-00899-2
7. Coubes, P, Roubertie, A, Vayssiere, N, Hemm, S, and Echenne, B. Treatment of DYT1-generalised dystonia by stimulation of the internal globus pallidus. Lancet. (2000) 355:2220–1. doi: 10.1016/S0140-6736(00)02410-7
8. Vidailhet, M, Vercueil, L, Houeto, J-L, Krystkowiak, P, Benabid, A-L, Cornu, P, et al. Bilateral deep-brain stimulation of the Globus Pallidus in primary generalized dystonia. N Engl J Med. (2005) 352:459–67. doi: 10.1056/NEJMoa042187
9. Fisher, R, Salanova, V, Witt, T, Worth, R, Henry, T, Gross, R, et al. SANTE study group, electrical stimulation of the anterior nucleus of thalamus for treatment of refractory epilepsy. Epilepsia. (2010) 51:899–908. doi: 10.1111/j.1528-1167.2010.02536.x
10. Peltola, J, Colon, AJ, Pimentel, J, Coenen, VA, Gil-Nagel, A, Ferreira, AGC, et al. Deep brain stimulation of the anterior nucleus of the thalamus in drug-resistant epilepsy in the MORE multicenter patient registry. Neurology. (2023) 100:e1852–65. doi: 10.1212/WNL.0000000000206887
11. Farrell, SM, Green, A, and Aziz, T. The current state of deep brain stimulation for chronic pain and its context in other forms of neuromodulation. Brain Sci. (2018) 8:158. doi: 10.3390/brainsci8080158
12. Sheth, SA, and Mayberg, HS. Deep brain stimulation for obsessivecompulsive disorder and depression. Annu Rev Neurosci. (2023) 46:341–58. doi: 10.1146/annurev-neuro-110122-110434
13. Wu, H, Hariz, M, Visser-Vandewalle, V, Zrinzo, L, Coenen, VA, Sheth, SA, et al. Nuttin, deep brain stimulation for refractory obsessive-compulsive disorder (OCD): emerging or established therapy? Mol Psychiatry. (2021) 26:60–5. doi: 10.1038/s41380-020-00933-x
14. Dayal, V, Limousin, P, and Foltynie, T. Subthalamic nucleus deep brain stimulation in parkinson’s disease: the effect of varying stimulation parameters. J Parkinsons Dis. (2017) 7:235–45. doi: 10.3233/JPD-171077
15. Wagle Shukla, A, Zeilman, P, Fernandez, H, Bajwa, JA, and Mehanna, R. DBS programming: an evolving approach for patients with parkinson’s disease. Parkinsons Dis. (2017) 2017:8492619. doi: 10.1155/2017/8492619
16. Hariz, M, and Blomstedt, P. Deep brain stimulation for parkinson’s disease. J Intern Med. (2022) 292:764–78. doi: 10.1111/joim.13541
17. Nanivadekar, A. C., Bose, R., Petersen, B. A., Okorokova, E. V., Sarma, D., Madonna, T. J., et al. Restoration of sensory feedback from the foot and reduction of phantom limb pain via closed-loop spinal cord stimulation, Nat Biomed Eng (2023). doi: 10.1038/s41551-023-01175-2 [Epub ahead of print].
18. Várkuti, B, Halász, L, Hagh Gooie, S, Miklós, G, Smits Serena, R, van Elswijk, G, et al. Conversion of a medical implant into a versatile computer-brain interface. Brain Stimul. (2024) 17:39–48. doi: 10.1016/j.brs.2023.12.011
19. Kuncel, AM, and Grill, WM. Selection of stimulus parameters for deep brain stimulation. Clin Neurophysiol. (2004) 115:2431–41. doi: 10.1016/j.clinph.2004.05.031
20. Kuncel, AM, Cooper, SE, and Grill, WM. A method to estimate the spatial extent of activation in thalamic deep brain stimulation. Clin Neurophysiol. (2008) 119:2148–58. doi: 10.1016/j.clinph.2008.02.025
21. Keane, M, Deyo, S, Abosch, A, Bajwa, JA, and Johnson, MD. Improved spatial targeting with directionally segmented deep brain stimulation leads for treating essential tremor. J Neural Eng. (2012) 9:046005. doi: 10.1088/1741-2560/9/4/046005
22. Brinda, A, Slopsema, JP, Butler, RD, Ikramuddin, S, Beall, T, Guo, W, et al. Lateral cerebellothalamic tract activation underlies DBS therapy for essential tremor. Brain Stimul. (2023) 16:445–55. doi: 10.1016/j.brs.2023.02.002
23. McIntyre, CC, Mori, S, Sherman, DL, Thakor, NV, and Vitek, JL. Electric field and stimulating influence generated by deep brain stimulation of the subthalamic nucleus. Clin Neurophysiol. (2004) 115:589–95. doi: 10.1016/j.clinph.2003.10.033
24. Butson, CR, and McIntyre, CC. Role of electrode design on the volume of tissue activated during deep brain stimulation. J Neural Eng. (2005) 3:1–8. doi: 10.1088/1741-2560/3/1/001
25. Chaturvedi, A, Butson, CR, Lempka, SF, Cooper, SE, and McIntyre, CC. Patient-specific models of deep brain stimulation: influence of field model complexity on neural activation predictions. Brain Stimul. (2010) 3:65–77. doi: 10.1016/j.brs.2010.01.003
26. Aström, M, Lemaire, J-J, and Wårdell, K. Influence of heterogeneous and anisotropic tissue conductivity on electric field distribution in deep brain stimulation. Med Biol Eng Comput. (2012) 50:23–32. doi: 10.1007/s11517-011-0842-z
27. Butson, CR, Cooper, SE, Henderson, JM, and McIntyre, CC. Patientspecific analysis of the volume of tissue activated during deep brain stimulation. NeuroImage. (2007) 34:661–70. doi: 10.1016/j.neuroimage.2006.09.034
28. Howell, B, and McIntyre, CC. Analyzing the tradeoff between electrical complexity and accuracy in patient-specific computational models of deep brain stimulation. J Neural Eng. (2016) 13:036023. doi: 10.1088/1741-2560/13/3/036023
29. Kennedy, DN, Haselgrove, C, Riehl, J, Preuss, N, and Buccigrossi, R. The three nitrcs: a guide to neuroimaging neuroinformatics resources. Neuroinformatics. (2015) 13:383–6. doi: 10.1007/s12021-015-9263-8
30. Li, X, Morgan, PS, Ashburner, J, Smith, J, and Rorden, C. The first step for neuroimaging data analysis: DICOM to NIfTI conversion. J Neurosci Methods. (2016) 264:47–56. doi: 10.1016/j.jneumeth.2016.03.001
31. Horn, A, Li, N, Dembek, TA, Kappel, A, Boulay, C, Ewert, S, et al. Lead-DBS v2: towards a comprehensive pipeline for deep brain stimulation imaging. NeuroImage. (2019) 184:293–316. doi: 10.1016/j.neuroimage.2018.08.068
32. Baniasadi, M, Proverbio, D, Gonçalves, J, Hertel, F, and Husch, A. FastField: an open-source toolbox for efficient approximation of deep brain stimulation electric fields. NeuroImage. (2020) 223:117330. doi: 10.1016/j.neuroimage.2020.117330
33. Avants, BB, Tustison, NJ, Song, G, Cook, PA, Klein, A, and Gee, JC. A reproducible evaluation of ANTs similarity metric performance in brain image registration. NeuroImage. (2011) 54:2033–44. doi: 10.1016/j.neuroimage.2010.09.025
34. Reich, MM, Horn, A, Lange, F, Roothans, J, Paschen, S, Runge, J, et al. Probabilistic mapping of the antidystonic effect of pallidal neurostimulation: a multicentre imaging study. Brain. (2019) 142:1386–98. doi: 10.1093/brain/awz046
35. Treu, S, Strange, B, Oxenford, S, Neumann, W-J, Kühn, A, Li, N, et al. Deep brain stimulation: imaging on a group level. NeuroImage. (2020) 219:117018. doi: 10.1016/j.neuroimage.2020.117018
36. Habets, JGV, Janssen, MLF, Duits, AA, Sijben, LCJ, Mulders, AEP, De Greef, B, et al. Machine learning prediction of motor response after deep brain stimulation in parkinson’s disease-proof of principle in a retrospective cohort. PeerJ. (2020) 8:e10317. doi: 10.7717/peerj.10317
37. Middlebrooks, EH, Okromelidze, L, Wong, JK, Eisinger, RS, Burns, MR, Jain, A, et al. Connectivity correlates to predict essential tremor deep brain stimulation outcome: evidence for a common treatment pathway. NeuroImage Clin. (2021) 32:102846. doi: 10.1016/j.nicl.2021.102846
38. Lange, F, Soares, C, Roothans, J, Raimundo, R, Eldebakey, H, Weigl, B, et al. Machine versus physician-based programming of deep brain stimulation in isolated dystonia: a feasibility study. Brain Stimul. (2023) 16:1105–11. doi: 10.1016/j.brs.2023.06.018
39. North, RB, Ewend, MG, Lawton, MT, and Piantadosi, S. Spinal cord stimulation for chronic, intractable pain: superiority of “multi-channel” devices. Pain. (1991) 44:119–30. doi: 10.1016/0304-3959(91)90125-H
40. Kumar, K, Nath, R, and Wyant, GM. Treatment of chronic pain by epidural spinal cord stimulation: a 10-year experience. J Neurosurg. (1991) 75:402–7. doi: 10.3171/jns.1991.75.3.0402
41. Deer, TR, Mekhail, N, Provenzano, D, Pope, J, Krames, E, Leong, M, et al. North, Neuromodulation Appropriateness Consensus Committee, the appropriate use of neurostimulation of the spinal cord and peripheral nervous system for the treatment of chronic pain and ischemic diseases: the neuromodulation appropriateness consensus committee. Neuromodulation. (2014) 17:515–50. doi: 10.1111/ner.12208
42. Tasker, RR, and Vilela Filho, O. Deep brain stimulation for neuropathic pain. Stereotact Funct Neurosurg. (1995) 65:122–4. doi: 10.1159/000098682
43. Mandat, V, Zdunek, PR, Krolicki, B, Szalecki, K, Koziara, HM, Ciecierski, K, et al. Periaqueductal/periventricular gray deep brain stimulation for the treatment of neuropathic facial pain. Front Neurol. (2023) 14:1239092. doi: 10.3389/fneur.2023.1239092
44. Johnson, KA, Dosenbach, NUF, Gordon, EM, Welle, CG, Wilkins, KB, Bronte-Stewart, HM, et al. Proceedings of the 11th annual deep brain stimulation think tank: pushing the forefront of neuromodulation with functional network mapping, biomarkers for adaptive dbs, bioethical dilemmas, ai-guided neuromodulation, and translational advancements. Front Hum Neurosci. (2024) 18:1320806. doi: 10.3389/fnhum.2024.1320806
45. Chandrasekaran, S, Nanivadekar, A, McKernan, G, Helm, E, Boninger, M, Collinger, J, et al. Sensory restoration by epidural stimulation of the lateral spinal cord in upper-limb amputees. eLife. (2020) 9:e54349. doi: 10.7554/eLife.54349
46. Nanivadekar, A, Chandrasekaran, S, Helm, E, Boninger, M, Collinger, J, Gaunt, R, et al. Closed-loop stimulation of lateral cervical spinal cord in upper-limb amputees to enable sensory discrimination: a case study. Sci Rep. (2022) 12:17002. doi: 10.1038/s41598-022-21264-7
47. Bach-Y-Rita, P, Kaczmarek, K, Tyler, M, and Garcia-Lara, J. Form perception with a 49-point electrotactile stimulus array on the tongue: a technical note. J Rehabil Res Dev. (1998) 35:427–30.
48. Bach-Y-Rita, P, and Kercel, S. Sensory substitution and the humanmachine interface. Trends Cogn Sci. (2003) 7:541–6. doi: 10.1016/j.tics.2003.10.013
49. Bach-Y-Rita, P. Tactile sensory substitution studies. Ann N Y Acad Sci. (2006) 1013:83–91. doi: 10.1196/annals.1305.006
50. Plaisier, MA, Sap, LI, and Kappers, AM. Perception of vibrotactile distance on the back. Sci Rep. (2020) 10:17876. doi: 10.1038/s41598-020-74835-x
51. Longin, L, and Deroy, O. Augmenting perception: how artificial intelligence transforms sensory substitution. Conscious Cogn. (2022) 99:103280. doi: 10.1016/j.concog.2022.103280
52. Bari, AA, Fasano, A, Munhoz, RP, and Lozano, AM. Improving outcomes of subthalamic nucleus deep brain stimulation in parkinson’s disease. Expert Rev Neurother. (2015) 15:1151–60. doi: 10.1586/14737175.2015.1081815
53. Rizzone, MG, Martone, T, Balestrino, R, and Lopiano, L. Genetic background and outcome of deep brain stimulation in parkinson’s disease. Parkinsonism Relat Disord. (2019) 64:8–19. doi: 10.1016/j.parkreldis.2018.08.006
Keywords: DBS programming, paresthesia, machine learning, computer-brain interfaces, prediction
Citation: Halász L, Sajonz BEA, Miklós G, van Elswijk G, Hagh Gooie S, Várkuti B, Tamás G, Coenen VA and Erōss L (2024) Predictive modeling of sensory responses in deep brain stimulation. Front. Neurol. 15:1467307. doi: 10.3389/fneur.2024.1467307
Edited by:
Thomas Schultz, University of Bonn, GermanyReviewed by:
Bowen Chang, Anhui Provincial Hospital, ChinaAndreas Husch, University of Luxembourg, Luxembourg
Copyright © 2024 Halász, Sajonz, Miklós, van Elswijk, Hagh Gooie, Várkuti, Tamás, Coenene and Erōss. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: László Halász, l.halasz@outlook.com