Cortical network organization reflects clinical response to subthalamic nucleus deep brain stimulation in Parkinson's disease.

The degree of response to subthalamic nucleus deep brain stimulation (STN-DBS) is individual and hardly predictable. We hypothesized that DBS-related changes in cortical network organization are related to the clinical effect. Network analysis based on graph theory was used to evaluate the high-density electroencephalography (HDEEG) recorded during a visual three-stimuli paradigm in 32 Parkinson's disease (PD) patients treated by STN-DBS in stimulation "off" and "on" states. Preprocessed scalp data were reconstructed into the source space and correlated to the behavioral parameters. In the majority of patients (n = 26), STN-DBS did not lead to changes in global network organization in large-scale brain networks. In a subgroup of suboptimal responders (n = 6), identified according to reaction times (RT) and clinical parameters (lower Unified Parkinson's Disease Rating Scale [UPDRS] score improvement after DBS and worse performance in memory tests), decreased global connectivity in the 1-8 Hz frequency range and regional node strength in frontal areas were detected. The important role of the supplementary motor area for the optimal DBS response was demonstrated by the increased node strength and eigenvector centrality in good responders. This response was missing in the suboptimal responders. Cortical topologic architecture is modified by the response to STN-DBS leading to a dysfunction of the large-scale networks in suboptimal responders.

INTRODUCTION

Deep brain stimulation of the subthalamic nucleus (STN‐DBS) is an effective and well‐established treatment of motor symptoms in Parkinson's disease (PD) (Schuepbach et al., 2013). As the STN is a part of the basal ganglia (BG)‐thalamocortical circuits that are also involved in various nonmotor functions, cognitive and affective activities can also be influenced by STN‐DBS (Bostan, Dum, & Strick, 2018). Recommended indication and exclusion criteria (Benabid, Chabardes, Mitrofanis, & Pollak, 2009) have reduced the risk of side effects; however, adverse effects still may complicate the otherwise successful therapy. Clinical practice has shown that also the degree of response to STN‐DBS treatment is individual and that the level of clinical improvement and occurrence of adverse effects in each patient is hardly predictable. In the previous study (Bočková et al., 2020), we have described a simple tool for identification suboptimal responders to STN‐DBS. We have reported that the prolongation of reaction time (RT) during “on” stimulation state as compared to DBS “off” state differentiates the suboptimal responders from the optimal responders with shortened RT during “on” stimulation state as compared to DBS “off.” In this study, we focused on cortical response to DBS. We analyzed the cortical bioelectrical activity using the network analysis approach that provides a macroscopic perspective of cortical connections. Network analysis describes the cortical topologic architecture and properties at the whole‐brain level. We hypothesized that DBS‐related changes in cortical network organization are related to the clinical effect. If so, they might also be used for identification of optimal and suboptimal responders to STN‐DBS.

METHODS

Subjects

We recorded a total of 40 patients with late motor complications treated with STN‐DBS (Activa PC or Libra XP stimulators); some patients were excluded from the final analysis for technical reasons. In the end, 32 PD patients were included in our analysis (see Table S1, Supporting Information). We excluded patients with signs of dementia. There was no further patient selection. All patients were informed about the nature of this study and gave their informed consent. The study received the approval of the local ethics committee (Ethics Board of the Faculty of Medicine, Masaryk University). The Unified Parkinson's Disease Rating Scale (UPDRS) at the time of experimental recording and a neuropsychological examination were used to evaluate each patient's current clinical condition. Patients did not express signs of dementia or major depression and did not have any other important disorders according to previous detailed neuropsychological examination. The neuropsychological tests performed during the experimental session were Stroop Test, Word List (Wechsler Memory Scale‐III), and Verbal Fluency Test. The therapeutic stimulation parameters were used for recording during the STN‐DBS “on” condition (see Table S1).

Experimental protocol and recordings

Recordings were performed in a Faraday shielded room with a constant temperature. A high‐density electroencephalography (HDEEG) Electrical Geodesics, Inc. (EGI) system with 256 channels was used for the scalp EEG recording. The sampling rate was 1 kHz with Cz as a reference electrode. For motor and cognitive testing, a visual oddball three‐stimuli protocol was used (Bočková et al., 2013; Polich, 2007) (see Table 1). The frequent (nontarget, standard) stimuli, which were 70% of all the stimuli, were small blue circles. These were not to be followed by any reaction. The target stimuli, which were 15% of all the stimuli, were larger blue circles, and the patients had to press a response button at the time of the target detection. The distractors (rare nontarget stimuli), which were 15% of the stimuli, were black and white checkerboards; no response was required. The interstimulus interval was 4 s. The duration of the stimulus exposure was 200 ms. Performance of the whole task lasted 14 min, containing 200 trials in total for each patient. The visual stimuli were presented in random order on a monitor. Patients were in the “off” medication condition (12‐hr medication withdrawal) and repeated the experimental paradigm in both STN‐DBS “on” and “off” states. Fourteen patients were recorded first in the “off” state and later in the “on” state, 18 patients were recorded first in the “on” condition and then in the “off” condition to exclude the effect of learning. EEG recording during DBS “off” state was started approximately 10–15 min after turning off the stimulator.

TABLE 1

Experimental paradigm

Stimulus	Description	Image	Response	Trials	Proportion (%)
Target	Large blue circle		Press a button	30	15
Nontarget	Small blue circle		No response	140	70
Distractor	Black and white checkerboard		No response	30	15

Note: Distractor‐related responses are related to frontal focal attention and working memory, and target‐related responses are related to temporal–parietal activity and subsequent memory processing (Polich, 2007).

Data analysis

The whole analytical process is shown in Figure 1 and described in the text. This process includes three phases: preprocessing, source reconstruction, and network analysis. Recordings from both conditions (DBS “on” and “off”) in each patient underwent the same processing steps. The statistical evaluation of the findings obtained from electrophysiological data was performed by nonparametric Wilcoxon tests. Neuropsychological battery and UPDRS scores were assessed by a multiway analysis of variance (ANOVA). Due to the exploratory nature of the analysis, we did not use correction for multiple comparisons in the evaluation of the differences in network parameters between DBS “on” and DBS “off” state.

FIGURE 1

EEG data preprocessing and analysis pipeline

Reaction times

Reaction times (RT) during target stimulation in “off” and “on” states were determined for each patient. RT was measured as the delay between the target stimulus and the patient's reaction (button press). Subsequently, the statistical difference between RT‐DBS‐“on” and RT‐DBS‐“off” in each patient was tested using the Wilcoxon rank‐sum test. Based on this difference and significance (p < .05), three groups of patients were distinguished. RT faster in DBS “on” (Group 1, optimal responders), RT faster in DBS “off” (Group −1, suboptimal responders), and no significant RT difference between DBS “on” and DBS “off” (Group 0). Response accuracy was measured as the percentage of correct responses. No significant difference between groups was found. The mean accuracy was 96%.

High‐density EEG preprocessing and artifact suppression

The EEG data were processed off‐line using the EEGLAB toolbox (Delorme & Makeig, 2004) complemented by an in‐house solution running under MATLAB 2014b (The MathWorks, Inc., Natick, MA). We reduced the number of channels to 204, discarding facial and neck‐line electrodes as those are usually contaminated with muscle artifacts (Coito, Michel, Vulliemoz, & Plomp, 2019). Data were filtered to 0.1–100 Hz bandwidth with Butterworth filter of second order, 12 dB/octave roll‐off and forward and backward passes. Residues of DBS‐related artifacts under 100 Hz were visually detected in the frequency domain (several sharp peaks with substantially higher magnitude than background activity) in each patient and suppressed by a fast Fourier transform (FFT) filter. Independent component analysis (ICA) was used to eliminate common artificial signals from blinking and ECG (no more than four components were suppressed). Data were visually inspected in SignalPlant software (Plesinger, Jurco, Halamek, & Jurak, 2016). Channels with frequent artifacts (less than 5%) were interpolated using the spherical spline method with polynomials of degree 1. Continuous EEG recording was segmented into trials lasting 3 s (1 s before and 2 s after stimuli onset). Bad trials were marked to be discarded from further analysis (fewer than 5% of the trials were excluded for each patient). Because of the reconstruction into the source space as a consecutive step, data were re‐referenced to the average reference.

Source reconstruction

Electrical source imaging (ESI) was performed in Cartool software (D. Brunet, cartoolcommunity.unige.ch), complemented by an in‐house solution running under MATLAB 2014b. The Montreal Neurological Institute (MNI) template was used for forward locally spherical model with anatomical constraints (LSMAC) model construction and low‐resolution electromagnetic tomography (LORETA) was used for inversion. Reconstructed signals were parcellated based on the automatic anatomical labeling (AAL) atlas (90 areas excluding cerebellum and vermis). Electrical dipoles from each area were projected to the refined average dipole orientation based on the approach by Coito, Michel, van Mierlo, Vulliemoz, and Plomp (2016). Only the centroid signal from each area was selected for subsequent analysis.

Electrode positions

DBS electrode positions were verified in all patients using Lead‐DBS software (www.lead-dbs.org; Horn & Kühn, 2015) (see Figure 2). Postoperative CT images were coregistered to preoperative MRI using Advanced Normalization Tools (ANTs). Images were normalized into MNI ICBM 2009b NLIN ASYM space by ANTs based on preoperative MRI. Each step was manually checked in each patient. DBS electrodes were then localized within MNI space and overlaid with DISTAL atlas (Ewert et al., 2018).

FIGURE 2

Deep brain stimulation electrode localization. DBS electrode location in the STN. Green, electrodes in patients from Group 1; yellow, electrodes in Group 0; red, electrodes in patients in Group −1. Notice that the grouping in the three groups is not related to the position of the electrodes

Network analysis

Ninety AAL regions of interest (ROIs) were used as network nodes. The edges were defined by a phase‐lag index (PLI) as a measure of phase synchronization (Stam, Nolte, & Daffertshofer, 2007) separately for each stimulus type and these frequency bands: 1–8, 8–20, 20–45, and 55–80 Hz; PLI equal to zero means no phase coupling; PLI equal to one means perfect phase locking. Wider frequency bands were selected because of the relatively short temporal window for PLI calculation compared to (Geraedts et al., 2018; Hatz, Meyer, Zimmermann, Gschwandtner, & Fuhr, 2017; Stam et al., 2007). Baseline interval in each trial was 600–100 ms before stimulus onset; PLI changes after stimulation were evaluated in the 500 ms window 200 ms after stimulus onset to capture the motor‐cognitive part of the response (Polich, 2007). Followed by Fisher's Z transformation (Lowe, Mock, & Sorenson, 1998), a weighted connectivity matrix was established for each patient, frequency band, and stimulus type, and analyzed by global and regional measures.

The Brain Connectivity Toolbox (Bassett & Sporns, 2017; Fornito, Zalesky, & Bullmore, 2016; Rubinov & Sporns, 2010) was used and the average node strength, average clustering coefficient, characteristic path length, and modularity were computed on a global level. Further, normalized measures of average clustering coefficient and characteristic path length were adapted—normalization ensured by dividing the unnormalized values by values computed and averaged for 50 networks of the null model (random network with preserved strength and degree distribution). On the regional level, node strength and eigenvector centrality were computed as measures of a ROI's importance in the network.

RESULTS

The differences in RT between DBS “on” and “off” conditions identified three subgroups: Group +1, containing 12 patients, was defined by decreased RT in DBS “on” state (p < .05); Group 0, containing 14 patients, showed no statistically significant difference in RT; and Group −1, comprising six patients, displayed negative effects of DBS: the RT were longer during the DBS “on” state (p < .05) (Bočková et al., 2020). Further analysis confirmed Group +1 as the optimal responders, Group 0 as responders, and Group −1 as suboptimal responders. The neuropsychological tests showed significantly decreased (F = 5.605, p = .009) results in semantic memory test WL3 (word list recognition) for Group −1, and this trend was also detectable in WL1 and WL2 (word list 1 – immediate verbal memory, 2 – delayed verbal memory) tests (see Figure 3). The clinical effect of DBS on motor symptoms as measured by UPDRS scores was also lower in Group −1. The UPDRS “on” scores were higher in Group −1 than in the other groups, and the “on”/“off” state UPDRS difference was significantly lower than in Group +1 (see Figure 3). On the other hand, Group +1 patients were the best responders to the DBS therapy: their “off” state scores were the highest and the “on” state scores were the lowest, and their “on”/“off” state UPDRS difference was significantly higher than Group 0 (p = .04) and Group −1 (p = .02).

FIGURE 3

Clinical parameters. Upper panels: Reaction times and UPDRS differences. Each box covers the data from 25th to 75th percentiles, the red line in each box represents the median of explained variability over subjects, and whiskers represent 1.5 times the interquartile range (IQR). Red crosses show the outliers. Black stars represent significant differences using the nonparametric Wilcoxon rank‐sum test (p < .05). Lower panel: Semantic memory test results. Crosses show score means, Whiskers represent standard error. Semantic memory test results were statistically analyzed using a multiway analysis of variance (ANOVA). Difference in memory testing results in Group −1 was significant in the word semantic memory test – word list 3. The difference in the word list 1 did not reach statistical significance, but the trend to decreased performance is evident. Green, Group 1; yellow, Group 0; and red, Group −1

Global network organization

There were no significant changes in global network measures in the higher 8–20 and 20–45 Hz frequency bands in the DBS “on” condition as compared to the DBS “off” condition. However, global connectivity in Group −1 was impaired in the 1–8 Hz band in reaction to target stimulus in the DBS “on” condition compared to the DBS “off”; no changes were observed in Groups +1 and 0. The average node strength was decreased as well as the average clustering coefficient and the characteristic path length was increased in the DBS “on” state in Group −1, compared to DBS “off” (Figure 4). In the other groups, the opposite trend was noticeable, although not significant. In the 55–80 Hz bandpass, we detected different global network organization changes in “on” and “off” states as a reaction to all three‐stimuli types in Groups −1 and 0 in contrast to Group +1, where no STN‐DBS influence was observed in the gamma range connectivity; see Tables 2 and S3 summarizing all significant global network changes.

FIGURE 4

Global connectivity measures – average node strength in 1–8 Hz band. Each box covers the data from 25th to 75th percentiles, the red line in each box represents the median of explained variability over subjects, and whiskers represent 1.5 times the interquartile range (IQR). Red crosses show the outliers. Black stars represent significant differences using the nonparametric Wilcoxon rank‐sum test (p < .05). Decreased global connectivity (decreased node strength and clustering coefficient, increased characteristic path length) was observed in Group −1 patients during stimulation in contrast to the other groups where DBS produced no significant change

TABLE 2

Significant (Wilcoxon rank sum test, p < .05) changes in global network measures

	1–8 Hz		55–80 Hz
Group −1	Target		Target
	↓ w (p = .015) C (p = .015)	↑ L (p = .009)	↓ L (p = .041)	↑ w (p = .041) C (p = .041) Q (p = .026) λ (p = .015)
			Distractor
				↑ w (p = .041) C (p = .041)
			Frequent
				↑ λ (p = .041)
Group 0	–		Target
			↓ w (p = .041) C (p = .046)	↑ Q (p = .026) γ (p = .001) λ (p = .018)
			Distractor
			↓ w (p = .029) C (p = .026)	↑ L (p = .018) Q (p = .029) γ (p = .001) λ (p = .005)
			Frequent
			↓ w (p = .037) C (p = .033)	↑ L (p = .037) Q (p = .033) γ (p = .001) λ (p = .002)
Group +1	–		–

Note: The ↑ and ↓ stand for the increase/decrease in the DBS “on” state as compared to the DBS “off” state. Network measures are abbreviated as follows: w, average node strength; C, average clustering coefficient; L, characteristic path length; Q, modularity coefficient; γ, normalized average clustering coefficient; λ, normalized characteristic path length.

↓

↑

↓

↑

↑

↑

↓

↑

↓

↑

↓

↑

Local network organization

We have also found differences between “on” and “off” states in local network organization (using node strength and eigenvector centrality measures) among the three subgroups. Again, Group −1 displayed a different pattern compared to Groups 0 and +1. The most important finding was the significantly decreased node strength in the 1–8 Hz frequency range in several mainly frontal areas during DBS “on” state after target stimulus; this was not found in the other groups. In contrast, Group +1 patients displayed significantly increased node strength in the same frequency in the left and right supplementary motor area (SMA) in the “on” state after target stimulus. Group 0 was without changes in these parameters and frequency band, but we detected a significantly increased eigenvector centrality in the 8–20 Hz range in the left and right SMA and in the left gyrus parietalis superior and left paracentral lobule; see Figure 5 and Tables S2 and S4 with all regional differences. The increase in node strength in the SMA in Group +1 was equal to higher PLI coefficients on average to all ROIs (SMA is on average more phase‐locked with signals from other ROIs in the “on” state than in the “off” state). On the other hand, an increase in EC in the SMA in Group 0 suggests a closer connection (better synchronization) to other EC‐important ROIs in the “on” state than in the “off” state.

FIGURE 5

Local network measures in 1–8 and 8–20 Hz bands. The increased importance of the SMA was observed in Group 1 (increased node strength in 1–8 Hz) and Group 0 (increased eigenvector centrality in 8–20 Hz) during DBS. In contrast, this DBS‐related SMA change was absent in Group −1, where a decrease in importance was present (decreased node strength in 1–8 Hz) in several mainly frontal areas during the DBS “on” state. EC, eigenvector centrality; w, node strength

Again, there were minimal changes in local network measures observed in the higher frequency ranges in Group +1 in contrast to the other two subgroups.

DISCUSSION

The general aim of the study was to identify changes in cortical connectivity induced by STN‐DBS. We examined changes related to the performance of a cognitive‐motor task. We have also focused on the stimulation related interindividual differences in clinical outcomes and network organization. To study the influence of STN‐DBS on motor and cognitive circuits, a three‐stimuli visual experimental paradigm was used (Bočková et al., 2013; Polich, 2007). We focused on responses to target stimuli inducing a motor response that is associated with complex cognitive processes. All patients profited from STN‐DBS in motor symptom improvement and none of them manifested severe neuropsychiatric symptoms or global cognitive decline. There was variability in RT of the motor response to the target stimulus. Reaction times in PD depend on both motor and nonmotor components and have been used to study the efficacy of the drug treatment as well as of STN‐DBS and Gpi‐DBS (Jordan, Sagar, & Cooper, 1992; Kojovic et al., 2014; Kumru, Summerfield, Valldeoriola, & Valls‐Solé, 2004). RT deficits in PD were reported to be compatible with a deficit in higher‐order processes, including both cognitive and motor responses to a stimulus. Prolonged simple reaction time was proposed as a biomarker for cognitive impairment in PD (Cosgrove et al., 2016). Based on differences in RT, the patients were divided into three subgroups (Bočková et al., 2020). In one group of 12 patients, the RT as expected was shortened under the DBS “on” condition (Group +1). In the group of 14 patients, no significant difference was observed (Group 0). This subgroup may have a slightly lower response to STN‐DBS than Group +1, but the differences are mostly not significant and we consider both subgroups as good responders. Even though all the patients profited from motor improvement after DBS, a subgroup of six persons with suboptimal responses to DBS could be defined (Group −1). An unexpected prolongation of RT during “on” stimulation state as compared to DBS “off” state was identified. This reverse motor reaction to DBS appeared to be a manifestation of a complex response to DBS that differentiates the suboptimal responders from the good responders. The effect of DBS on motor PD symptoms was less expressed than in other patients. A neuropsychological examination showed a decrease in semantic memory test scores as compared to the other subgroups (Figure 3). As the UPDRS scores also improved in the suboptimal responders, albeit to a lesser degree, we presume that the prolongation of RT to target stimulus under the DBS “on” condition was caused by a dysfunction of the cognitive rather than of the motor part in this cognitive‐motor paradigm. In the previous study (Bočková et al., 2020), focused on oscillatory changes in the areas of interest linked to the experimental paradigm (temporal and parietal cortex, premotor regions, and the thalamus), the majority of patients (Groups +1 and 0) expressed and enhanced alpha and beta power decreases in the DBS “on” stimulation state as compared to the DBS “off” stimulation state. A power decrease (event‐related desynchronization [ERD]) is considered to be an expression of cortical activation (Pfurtscheller, 2001). This reactivity correlated with improved motor‐cognitive functioning during the “on” condition in Groups 0 and +1. An opposite reaction pattern was observed in Group −1, where longer RT were linked with lower ERD; that is, diminished cortical activation mainly in the alpha bandpass, known to be associated mainly with nonmotor functions.

The network analysis approach indicates a dysfunction of the large‐scale cerebral networks in suboptimal responders. The Group −1 patients displayed decreased average node strength and with it a related decrease (due to the weighted approach to the computation of network measures) in the average clustering coefficient and an increase in the characteristic path length in the 1–8 Hz band during the DBS “on” state. This means decreased connectivity and slower communication within the network. The reduction in overall functional connectivity has been described as being linked to PD‐related dementia (Ponsen, Stam, Bosboom, Berendse, & Hillebrand, 2013) and therefore we suspect that the development of cognitive deterioration is more probable in the suboptimal responders than in the optimal responders. A long‐term follow‐up study would be needed to confirm or reject this hypothesis.

The functional role of the SMA was enhanced in the DBS “on” state; in Group +1 we detected an increased node strength for the left and right SMA in the 1–8 Hz range; in Group 0, we observed an increased eigenvector centrality for the left and right SMA in 8–20 Hz range. Anatomic and functional connections between the STN, SMA, and frontal cortical structures are crucial in cognitive control over motor actions (Aron, Behrens, Smith, Frank, & Poldrack, 2007). The SMA is known to be functionally coupled to the STN mainly in the beta frequency band, and beta activity underlies the main motor symptoms in PD (Brown, 2003; Litvak et al., 2011). The time‐frequency analysis of SMA after stimulus onset computed in our previous study (Bočková et al., 2020) showed a dominant cluster of power change covering the alpha, low, and high beta bands (Figure 1). The relationship between power changes and connectivity strength (Klimes et al., 2018) indicates that beta synchronization in SMA can be linked to eigenvector centrality increase in 8–20 Hz. Frequencies in the theta‐alpha range in the STN functional couplings with frontal and temporo‐parietal areas reflect cognitive processes (Litvak et al., 2011; Zavala et al., 2014; Zavala et al., 2016). There was no DBS‐related difference in the area of the SMA in Group −1 patients. Moreover, a widespread decrease in node strength of several frontal areas and decreased global connectivity in low frequencies was observed in this group, indicating that the topological organization in particular of the frontal cortex was disturbed in patients with suboptimal responses. Reduced functional connectivity and decreased local integration in the alpha frequency band were reported in PD dementia and nonmotor symptom severity. These proposed markers might be useful in the screening process for DBS (Geraedts et al., 2018; Utianski et al., 2016). The disrupted topological organization of brain networks reflects decreased information transmission efficiency in patients with suboptimal response to DBS.

The stimulation may lead to the modification of the physiological balance in network functioning and could impair behavior and cognitive performance (Brittain, Sharott, & Brown, 2014); this was probably the case with Group −1. On the other hand, the best responders in our study, Group +1 patients, displayed no changes in overall connectivity, and minimal changes in local connectivity were also detected in lower and higher gamma bandpass. The suboptimal reactivity and clinical response in the Group −1 patients was not linked to the electrode position (see Figure 2), drug therapy, or DBS parameter setting; see Supporting Information.

RT evaluation and surface EEG can be used for identifying responsiveness to STN‐DBS. Our study confirmed the hypothesis of cortical network dysfunction in suboptimal responders to STN‐DBS. This study also confirmed that network analysis reflects the clinical outcomes of the DBS treatment and may reveal important information about brain dysfunction in neurodegenerative disorders such as PD. As PD is a heterogenous disease with different clinical courses and various phenotypes (motor and nonmotor), network analysis could help reveal different connectivity profiles in PD patients that could serve in clinical practice as predictive biomarkers. Further electrophysiological studies focused on network biomarkers are of recent high research interest (Litvak, Florin, Tamás, Groppa, & Muthuraman, 2021).

Limitations of the study

There is an ongoing discussion about the reliability of reconstructing electrophysiological activity in the deep brain areas from scalp recordings. Recently, a positive response to this issue was presented (Seeber et al., 2019), showing direct evidence that resting‐state scalp EEG data can detect subcortical activity. The electrical source imaging pipeline used in our work is very similar. The second limitation is linked to the small number of patients, mainly in Group −1 (n = 6). The groups of patients are based on RT in a visual motor cognitive paradigm. A long‐term clinical observation is necessary to evaluate the exact clinical outcomes in these patients. For these reasons, we consider our data as preliminary, needing to be confirmed and further developed in a larger prospective study.

CONCLUSION

Network organization changes after STN‐DBS in PD using HDEEG analysis correspond to the responsiveness and clinical outcomes of this advanced therapy. It remains to be clarified in future studies if such analysis could also serve as a predictive preoperative biomarker for clinical practice.

CONFLICT OF INTEREST

The authors declare no potential conflict of interest.

ETHICS STATEMENT

The study received the approval of the local ethics committee (Ethics Board of the Faculty of Medicine, Masaryk University). All subjects were informed about the nature of this study and gave their informed consent.

Appendix S1. Supporting Information

Click here for additional data file.