Effects of continuous theta-burst stimulation of the primary motor and secondary somatosensory areas on the central processing and the perception of trigeminal nociceptive input in healthy volunteers.

Noninvasive modulation of the activity of pain-related brain regions by means of transcranial magnetic stimulation promises an innovative approach at analgesic treatments. However, heterogeneous successes in pain modulation by setting reversible "virtual lesions" at different brain areas point at unresolved problems including the optimum stimulation site. The secondary somatosensory cortex (S2) has been previously identified to be involved in the perception of pain-intensity differences. Therefore, impeding its activity should impede the coding of the sensory component of pain intensity, resulting in a flattening of the relationship between pain intensity and physical stimulus strength. This was assessed using inactivating spaced continuous theta-burst stimulation (cTBS) in 18 healthy volunteers. In addition, cTBS was applied on the primary motor cortex (M1) shown previously to yield moderate and variable analgesic effects, whereas sham stimulation at both sites served as placebo condition. Continuous theta-burst stimulation flattened the relationship between brain activation and stimulus strength, mainly at S2, the insular cortex, and the postcentral gyrus (16 subjects analyzed). However, these effects were observed after inactivation of M1 while this effect was not observed after inactivation of S2. Nevertheless, both the M1 and the S2-spaced cTBS treatment were not reflected in the ratings of the nociceptive stimuli of different strengths (17 subjects analyzed), contrasting with the clear coding of stimulus strength by these data. Hence, while modulating the central processing of nociceptive input, cTBS failed to produce subjectively relevant changes in pain perception, indicating that the method in the present implementation is still unsuitable for clinical application.

1. Introduction

Pain serves as an evolutionary highly important warning system to avoid possible tissue damage. However, after fulfilling this purpose, pain may become a major threat to an individual's health and quality of life, in particular, when developing towards persistence. Hence, persistent pain is listed among the major health care issues defined by the World Health Organization. Its high prevalence[12,19,22,64,65] points at insufficient treatment options. The uncovering of the highly complex pathophysiological processes underlying pain[43] has become a major basis for the development of novel analgesic drugs.[56] However, clinical experience and scientific evidence suggest that nonpharmacological treatments my offer alternative or additional analgesic therapy options.

Among nonpharmacological approaches, transcranial magnetic stimulation (TMS) seems to be particularly suited because it enables the direct noninvasive modulation of the function of pain-relevant cerebral structures.[17,69,89,93] Transcranial magnetic stimulation is increasingly used in cognitive neuroscience because of its ability to impede task-specific neuronal activities by inducing a reversible so-called “virtual lesion.”[72,78] Indeed, the application of TMS stimulating mainly the primary motor cortex (M1) has been shown to yield moderate and variable analgesic effects against clinical or experimentally induced pain,[4,6,10,18,23,35,48-51,62,68,69,75,77] whereas the stimulation of the dorsolateral prefrontal cortex mainly influenced affective and attentional dimensions of pain.[68] By contrast, stimulation of the primary somatosensory cortex (S1) produced no effect on pain.[4,34,75] This heterogeneity in results stresses the need for further research on the utility of TMS as an analgesic treatment.

This study was based on the hypothesis that impeding the activity of the secondary somatosensory area (S2) should impede the coding of the sensory component of pain intensity. This hypothesis is supported by the results of a functional magnetic resonance imaging (fMRI) study showing that S2 is involved in the decoding of quantitative changes in the physical strength of nociceptive stimuli.[71] Hence, a reduction of the function of S2 should result in a flattening of the relationship between pain intensity, or pain stimulus–associated brain activation, and the physical stimulus strength. The hypothesis was investigated through an fMRI study on healthy volunteers, using a nociceptive stimulus, ie, intranasal application of gaseous carbon dioxide (CO2) that had previously been shown to particularly activate S2.[41]

2. Methods

2.1. Subjects and study design

The study followed the Declaration of Helsinki and was approved by the Ethics Committee of the Goethe-University Frankfurt am Main, Germany (protocol number 250/11). Written informed consent was obtained from each subject before participating in the study. Effects of paired continuous theta-burst stimulation (cTBS) on the processing and perception of nociceptive trigeminal input were investigated in 18 healthy right-handed volunteers (mean age ± SD, SD: 24.2 ± 2.5 years; 7 men). All subjects had been identified as responders to TMS in a test session before the actual experiments, as published previously.[33] The subjects' health status was ascertained by medical history and physical examination. Inclusion criteria were an age between 18 and 50 years and no relevant medical history, whereas exclusion criteria comprised current diseases, drug intake during the week preceding the experiment with the exception of oral contraceptives, and the presence of contraindications for TMS and fMRI. Alcohol was prohibited for a period of 24 hours before the experiments.

The study was based on a double-blind 3-way crossover design (Fig. 1), using either (1) cTBS of the primary motor cortex (M1) hand representation (“M1 verum”), or (2) cTBS of the secondary somatosensory cortex S2 (“S2 verum”), or (3) sham stimulation as the 3 main experimental conditions. The stimulation conditions were assessed at 3 different days separated by an interval of at least 1 week to avoid carryover effects, considering that for cTBS, cortical modulation lasted at least 50 minutes[91] while in chronic pain conditions, rTMS induced a significant pain decrease for up to 8 days.[50] At each day, the stimulation coils were placed successively at position M1 and S2, independent of the type of stimulus, to preserve blindness. The sham coils generated very similar clicks as the active coils, which were sensed also on the scalp; hence, the subjects were unable to point at which one was the sham or the verum stimulation.

Figure 1.

Schematic representation of the study design. The effects of deactivating cTBS transcranial stimulation on the pain-intensity ratings of intranasal nociceptive stimuli (500-ms pulses) of gaseous carbon dioxide (CO2) at 4 different concentrations (35, 50, 65, or 80% vol/vol) were assessed using 3 main experimental conditions comprising either (1) cTBS of the primary motor cortex (M1) hand representation (“M1 verum”), or (2) cTBS of the secondary somatosensory cortex S2 (“S2 verum”), or (3) sham stimulation. The experimental conditions were assessed at random order at separate days scheduled at an interval of at least 1 week. During each study day, 2 experimental sessions with nociceptive stimulation, pain-intensity ratings, and fMRI imaging were performed. After the baseline session, cTBS (or sham) was applied according to the respective study condition. During each session, subjects received 80 CO2 stimuli (20 stimuli for each of 4 different and randomly distributed strengths) at long randomly spaced intervals of 13.5 to 28.2 seconds. After a randomized interval of 3.4 to 6.6 seconds, the stimulus intensity is estimated on a visual analogue scale (enlarged section top left). cTBS, continuous theta-burst stimulation; fMRI, functional magnetic resonance imaging.

2.2. Transcranial magnetic stimulation

2.2.1. Subject preselection and individual localization of M1 and S2

Subjects were selected based on their response to TMS targeted at the right hemispheric M1 corresponding to the representation of the first dorsal interosseous (FDI) muscle of the left hand assessed in a separate session before the actual experiments and as reported in all detail previously.[33] In brief, the raw surface electromyography (EMG) signal was amplified, filtered (bandpass of 20-2000 Hz; Counterpoint Mk2, Dantec Elektronik, Skovlunde, Denmark), analog-to-digital converted at a sampling rate of 5 kHz (CED Micro 1401; Cambridge Electronic Design, Cambridge, United Kingdom), and stored in a computer for online visual inspection and offline analysis. Focal TMS of the hand area of the right M1 was performed with a figure-of-eight coil (Cool-B65; MagVenture, Farum, Denmark, diameter of each wing 65 mm) connected to a MagPro X100 magnetic stimulator (MagVenture) using a monophasic current waveform to induce a posterior-anterior current in the brain. The optimal coil position to elicit motor evokes potentials (MEPs) in the left FDI was defined as the site where TMS at suprathreshold stimulus intensity consistently produced the largest MEP. This “hotspot” was marked with a soft-tipped pen on a swimming cap on the subject's head to assure reproducible positioning of the coil throughout the experiments. The coil was held tangentially to the scalp with the handle pointing backwards and 45° away from the midline. The stimulation intensity was then adjusted to evoke MEPs in the left FDI with peak-to-peak amplitudes of 1.0 ± 0.3 mV. This intensity was used for all subsequent MEP recordings in the same subject.

2.2.2. Individual localization of M1 and S2

In addition, to obtain precise localizations of M1 and S2, a separate fMRI session was performed. During this session, the subjects received intranasal nociceptive stimuli (25 pulses of gaseous CO2 at a concentration of 65% vol/vol; for details, see below) at an interval of approximately 25 seconds. Subsequent fMRI data analysis provided the localization of S2. Then, 2 different images were alternatingly presented to the subject on a computer screen. In detail, a green image showing the written command “move” instructed the subject to extend the left index finger repeatedly at a frequency of approximately 1 Hz, whereas an image showing the written command “hold still” instructed the subject to rest the finger until the green image appeared again. Subsequent fMRI data analysis provided the localization of M1. These measurements took place at least 1 week before the main experiments started.

2.2.3. Induction of neuronal plasticity using continuous theta-burst stimulation

The individual localizations of M1 and S2 were transferred to the TMS navigation system (LOCALITE GmbH, Schloss Birlinghoven Sankt Augustin, Germany), enabling the exact alignment of the stimulating magnetic field and its precise orientation toward M1 or S2. For this purpose, the positions of the TMS coil and the subject's head were recorded by a camera system using reflectors, which had been attached and displayed on top of the previously acquired fMRI. An important maneuver was to drive a pointer around the subjects' head, which had to be detected by the camera and to be adjusted to the anatomical MRI data for surface registration. This allowed adjusting the TMS coil accurately to the brain target areas. These 2 areas were marked with a soft-tipped pen on a swimming cap on the subject's head to assure reproducible placement of the coils throughout the measurements. Neuronavigation and cTBS application took place in different rooms. To save time between cTBS and MRI scans, cTBS was applied in the MRI area without navigation. The navigated coil targets were marked on the swimming cap to ensure the same coil positions for the actual cTBS application.

Continuous theta-burst stimulation was delivered by means of a MagPro X100 magnetic stimulator connected to a 65-mm figure-of-eight coil (Cool-B65; MagVenture) using a biphasic current waveform (first phase: anterior-posterior–induced current in the brain; second phase: posterior-anterior–induced current) with the coil placed tangentially to the scalp with the handle pointing backwards and 45° away from the midline. The active motor threshold was defined as the lowest stimulus intensity, which elicited small MEPs (≥200 μV) in at least 5 of 10 consecutive trials during a slight tonic contraction (approximately 20% of the maximal strength) of the left FDI muscle, using the relative frequency method.[32] Active motor threshold values were rounded to the nearest 1% of maximum stimulator output (MSO). Stimulus intensity was set at 80% of active motor threshold as determined with the biphasic pulse configuration. The mean (±SD) across all subjects was 34.2 ± 5.9% MSO. A total of 600 pulses was delivered, using an established stimulation pattern with 3-pulse bursts at 50 Hz, and bursts were repeated every 200 milliseconds,[38] providing a total duration of cTBS of 40 seconds. A repeat of the cTBS took place after 10 minutes (spaced cTBS[29]) before the subject underwent a further fMRI scan for postbaseline. Spaced cTBS was chosen as protocol because it had demonstrated in previous work more effective long-term depression-like decrease of corticospinal excitability after M1 stimulation, and resistance to dedepression, when compared with a single cTBS train.[29] The coil was connected to a water cooling system kept at a constant temperature of 19°C. Sham cTBS was delivered using a matched, air-cooled sham TMS coil (MCF-P-B65 placebo-coil; MagVenture). The sham coil's magnetic shield provides a field reduction of approximately 80%. The placement of the sham coil was similar to the regular Cool-B65-coil during active stimulation. Using the same stimulation pattern as in the active condition ensured a replication of the look and sound of active cTBS.

2.3. Nociceptive stimulation

For nociceptive stimulation, a chemosensory pain model[47] was applied. This pain model has been successfully used for more than 30 years in pharmacological pain studies, starting from 1985,[47] and has been described in a recent computationally aided analytical review of human pain models as one of the best predictive models for clinical effects of analgesic drugs.[70] The model as applied in this study was based on the application of short pulses (500 ms) of gaseous CO2 to the subject's right nostril by means of an olfactometer (OM/2; Burghart Messtechnik GmbH, Wedel, Germany). The pulses were embedded in a continuously flowing airstream of 8 L/min at controlled temperature (36.5°C) and humidity (80% relative humidity) to avoid concomitant excitation of thermal or mechanical sensors.[47] In this study (Fig. 1), each subject received intranasally 80 CO2 stimuli (20 stimuli for each of 4 different and randomly distributed strengths: 35, 50, 65, and 80% vol/vol) at long randomly spaced intervals of 13.5 to 28.2 seconds to minimize habituation effects.[39] After each stimulus, subjects rated its pain intensity by means of a visual analog scale (VAS), displayed randomly within 3.4 to 6.6 seconds (mean: 4.9 seconds) after stimulus application and ranging from 0 (“no pain”) to 100 (“pain experienced at maximum”). Each session of nociceptive stimulation, intensity rating, and fMRT data acquisition (see below) took approximately 27 minutes.

2.4. Acquisition of functional magnetic resonance imaging data

During the session of nociceptive stimulation, an event-related design was used for fMRI data acquisition. The blood oxygenation level–dependent response to each stimulus was recorded at a field strength of 3 T with a dedicated head scanner (Magnetom Allegra; Siemens Medical Solutions, Erlangen, Germany), equipped with a 4-channel transmit–receive head coil. To reduce motion artifacts, the subject's head was immobilized using foam pads. For acquisition of fMRI data, a T2*-weighted gradient-echo (GE) echo-planar imaging (EPI) sequence with the following parameters was used: TR = 2048 milliseconds, TE = 30 milliseconds, flip angle = 90°, echo spacing = 420 μs, matrix size = 64 × 64, and in-plane resolution = 3 × 3 mm2. A total of 810 volumes were acquired per run, each of which comprised 32 slices with 3-mm thickness and an interslice gap of 1 mm, acquired in descending order; the first 5 volumes of each scanning block were discarded to ensure steady-state conditions. A total of 6 experimental sessions was performed for each subject to obtain measurements before and after each of 3 stimuli (cTBS of M1, cTBS of S2, sham, at 3 separate days with 2 sessions, ie, before and after cTBs at each day).

For subsequent off-line correction of distortions in the EPI images due to inhomogeneities of the static magnetic field B0,[3,40] magnetic field mapping was performed through GE imaging with identical geometric parameters, and 2 different TE values (4.89 and 7.35 ms) from which magnitude images and a phase difference map were calculated directly on the scanner. In addition, a T1-weighted anatomical data set with 1-mm isotropic resolution was acquired for each subject, using a 3-dimensional (3D) magnetization prepared rapid gradient echo (MP-RAGE)[66] sequence with the following parameters: TR = 2200 milliseconds, TE = 3.93 milliseconds, flip angle = 9°, TI = 900 milliseconds, FOV = 256 × 256 mm2, one slab with 160 sagittal slices of 1-mm thickness, using parallel acquisition (GeneRalized Autocalibrating Partial Parallel Acquisition, GRAPPA[31]) with an acceleration factor of 2 in phase encoding direction and a duration of 4 minutes.

2.5. Data analysis

2.5.1. Analysis of functional magnetic resonance imaging data

Functional MRI data could be analyzed from 16 subjects because of data loss in the remaining 2 volunteers. Data were analyzed using the SPM12 Matlab toolbox (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/; Wellcome Department of Imaging Neuroscience, London, United Kingdom[24,92]) on Matlab 2017b (Mathworks, Natick, MAUS) on an Intel Core i9 computer (operating system: Ubuntu Linux 17.10.1 64-bit). During each session, 810 volumes were acquired, which resulted in a total of 4860 volumes for all 6 runs acquired from each subject. The analysis of the brain activations triggered by the CO2 stimuli was performed in 2 steps comprising (1) image preprocessing followed by (2) first- and second-level statistical analysis for the identification of pain stimulus–related brain activations and their modulation by the 3 different TMS treatments (cTBS of M1, cTBS of S2, and sham).

2.5.1.1. Functional magnetic resonance image preprocessing

The data preprocessing pipeline integrated removal of the first 5 volumes in the scan series to avert T1 equilibration effects, realignment of all volumes to the first volume to correct for subject motion, and unwarping using the respective field map. Subsequently, data were corrected for acquisition time (slice timing). The high-resolution T1-weighted anatomical image was co-registered to the mean-EPI (created during the realign and unwarp process), segmented and normalized using 4th-degree B-spline interpolation to obtain image voxel sizes of 2 × 2 × 2 mm3. The resulting spatial normalization parameters were applied to the volumes of the EPI sequence that were subsequently smoothed with an isotropic 9-mm full-width half-maximum Gaussian kernel.

2.5.1.2. Functional magnetic resonance imaging data analysis

The brain activations after administration of the nociceptive trigeminal CO2 stimuli were analyzed by adding to the design matrix a first-order parametric regressor modulating the stick functions for either of the CO2 concentrations (order: 35, 50, 65, 80% vol/vol). Furthermore, the 6 rotational and translational parameters from the rigid body transformation, obtained during image realignment, were modeled as covariates of no interest. All regressors were convolved with the canonical hemodynamic response function. Low-frequency fluctuations of the MR signal were removed with a high-pass filter with a cutoff at 128 seconds. Voxelwise regression coefficients for all regressors were estimated using least squares within SPM12.

After model estimation, the effects of interest were tested by linear contrasts, generating statistical parametric maps of t-values for each subject. The design of this analysis directly addressed the main hypothesis of this study. Therefore, for brain activations associated with the different strengths of the 4 CO2 stimuli, contrast images were created separately for each experimental condition and measurement. The contrast [-3 -1 1 3] was used, which models a linear relationship between stimulus strength and associated brain activations, yielding a single-contrast image per subject. These contrast images were generated in a first-level analysis and submitted to a group-level analysis.

The group-level analysis of the TMS effects on brain activation followed a factorial 2 × 3 analysis of variance (ANOVA) design, testing the influence of the factors “measurement” (baseline, post-cTBS) and “condition” (cTBS treatment: M1 verum, S2 verum, sham) on the linear relationship between stimulus strength and associated brain activation as found in the first-level analysis. Contrasts for main effects and the interaction terms TMS by measurement were calculated. These analyses were performed using the “full-factorial” design provided by SPM12. The contrast of main interest consisted of the interaction between the ANOVA factors. In the case of significant effects, post hoc t-contrast was performed to specifically test the hypothesis that cTBS treatment flattened the linear response vs stimulus strength relationship for the M1 and S2 verum conditions against the placebo (sham TBS) condition, using the SPM t-contrasts [1 0 -1 -1 0 1] for the M1 verum and [0 1 -1 0 -1 1] for the S2 verum conditions. In each case, the entries in the design matrix correspond to measurement1/condition1, measurement1/condition2, measurement1/condition3, measurement2/condition1, measurement2/condition2, and measurement2/condition3, according to the factors defined above. The statistical parametric maps (SPMt) resulting from each analysis were interpreted with regard to the probabilistic behavior of Gaussian random fields. Results are reported at P < 0.01 (family-wise error-corrected) at cluster level with a cluster size threshold of 5 voxels. The localization of brain activations was aided by the anatomy toolbox, version 2.5.2.[21] Locations of significant peak activations are reported in Montreal Neurological Institute (MNI) coordinates (mm).

2.6. Analysis of intensity ratings of the nociceptive stimuli

Psychophysical data could be analyzed from 17 subjects because of data loss in the remaining volunteer. A total of 8160 VAS ratings was available for data analysis (17 subjects, each of whom received and rated 80 stimuli during each of 6 experiments, 17 × 80 × 6 = 8160). Data were evaluated in 3 different ways: (1) classical ANOVA, (2) unsupervised data exploration using cluster analysis and Gaussian mixture modeling, and (3) supervised machine learning using different techniques. These analyses had 2 major goals: (1) to establish the expected effect of the CO2 concentration as a positive control for the successful design and performance of the experiments. The results were used to demonstrate that the data provided a suitable basis for the second goal: (2) to assess an effect of the TMS condition. These analyses were performed using the R software package (version 3.4.3 for Linux; http://CRAN.R-project.org/).[20]

2.6.1. First evaluation: analysis of variance

The hypothesis underlying this study, ie, that cTBS treatment flattens the linear relationship between the pain stimulus strength and the pain response, was directly targeted by analyzing the changes in the slopes of this relationship. Therefore, the slopes of this curve were calculated by means of linear regression performed using the “lm” command implemented in the base package of the R software package. Subsequently, the slopes were submitted to an analysis of variance for repeated-measures ANOVA (rm-ANOVA) with the within-subject factors “condition” (cTBS treatment: M1 verum, S2 verum, and sham) and “measurement” (baseline and post-cTBS). The design of this analysis reflected the analysis of the fMRI data (see above). In addition, the data were submitted to a complete rm-ANOVA, with the within-subject factors “condition” and “measurement” (see above) and “CO2 concentration” (35, 50, 65, or 80% vol/vol). Calculations were performed using the R “aov” command implemented in the basic R software package “grDevices.” The α-level was set at 0.05.

2.6.2. Second evaluation: unsupervised data exploration for subgroup identification

Data were submitted to a cluster analysis to assess whether the pain ratings of the CO2 stimuli provided a group structure that agreed (1) with the CO2 concentrations or (2) with the TMS conditions. To study a cluster structure that agreed with the CO2 concentrations of the applied nociceptive stimuli, which served as a positive control, the data were arranged into a 68 × 6-sized matrix, ie, 17 subjects rating the pain of CO2 stimuli of 4 different strengths (17 × 4 = 68) during 3 experimental conditions always at baseline and after TMS stimulation (3 × 2 = 6). To study a cluster structure that agreed with the experimental conditions with respect to the TMS condition, the data were rearranged into a 51 × 8 matrix, ie, 17 subjects tested under 3 experimental conditions (51) rating the pain of CO2 stimuli at 4 different strengths at baseline and after TMS stimulation (4 × 2 = 8).

These data matrices were explored for a group structure using the Ward method[88] of hierarchical clustering[42] with the Euclidean distance among the subjects' VAS ratings. To identify the optimum number of clusters in the data space, a cluster stability score was computed addressing the consensus in cluster assignments across multiple runs of the clustering algorithm on data sets created by repeated random resampling from the original data set.[63] Specifically, the stability score captures the average proportion of observations not placed in the same cluster during repeated runs.[74]

Clustering, including stability assessments, was performed using the progeny algorithm,[36] which selects the optimum cluster number that renders the most stable clustering by evaluation of clustering stability starting with an initial clustering of the full data set, followed by bootstrapping and repetitive clustering. During resampling, the algorithm randomly sampled feature values with replacement to construct new samples, so-called “progenies,” rather than directly sampling existing samples as in more common algorithms, which has been shown to outperform some classical methods in discovering clinically meaningful patient groupings in biomedical data sets.[36] A number of k = [2…8] clusters were tentatively chosen, using 10 progenies in 10 randomly created data sets and 100 iterations, which corresponds to the defaults of the R software package “progenyClust” (https://cran.r-project.org/package=progenyClust).[37] To obtain SDs of the stability measures, the procedures were repeated 100 times. The final number of clusters was chosen on the basis of the stability score criterion.[36] The clustering result was visually assessed by drawing Silhouette plots that provide a graphical interpretation and validation of data clusters[76] as implemented in the R library “cluster” (https://cran.r-project.org/package=cluster).[59]

After establishment of the cluster number and cluster membership of the subjects, the agreement of the cluster structure with (1) the CO2 stimulus intensity or (2) the TMS condition was analyzed by submitting the cross-tabulated data (cluster membership vs previous classification, ie, stimulus concentration or experimental condition) to χ2 statistics.[73] The α-level was set at 0.05.

A further unsupervised analysis addressed possible group structures in the steepness of the linear relationship between pain and physical stimulus strength. Specifically, the cTBS-mediated inhibition of S2 processing had been expected to flatten the pain intensity vs the CO2 concentration curve.[71] This analysis was based on the slope of the respective curve obtained through linear regression as described above. In detail, the slope observed at baseline was subtracted from the slope observed after cTBS application, at an individual level and separately for each experimental condition. Thus, the expected effect of S2 inhibition would be reflected by negative values, ie, a numerically smaller slope after cTBS application as compared to the slope observed at baseline.

The slope differences were explored for a group structure that was expected to (1) reflect the experimental conditions by an over-representation of slope differences calculated for the S2 verum condition in groups with negative differences and in addition, (2) a possible responder vs nonresponder partition of the cohort, analogous to the observation in the preliminary experiments of the present project.[33] Specifically, the distribution of the changes in the slope of the linear relationship between pain intensity and CO2 concentration was investigated by analyzing the probability density function (PDF) as described previously.[55,84] In brief, the Pareto density estimation, ie, a kernel density estimator particularly suited for the identification of groups in the data,[83] was used. A multimodal distribution of the slope differences was assessed by fitting a Gaussian mixture model (GMM) to the Pareto density estimations as , where N(x|mi, si) denotes Gaussian probability densities (components) with mean values mi and SDs si. The wi denote the mixture weights indicating the relative contribution of each Gaussian component to the overall distribution, which add up to a value of 1. M denotes the number of components in the mixture. Gaussian mixture model fitting was performed with the R software package “AdaptGauss” (https://cran.r-project.org/package=AdaptGauss).[84] To determine the optimum number of components, model optimization was performed for M = 1 to 5 components. The final model was selected on the basis of likelihood ratio tests.[80] Subject association to the identified subgroups was obtained using the Bayes theorem[7] that provided the probability that an individual observation belongs to mode i calculated as the posterior probability. In the case of a relevant group structure, associations with the experimental conditions were statistically analyzed through χ2 tests.[73]

2.6.3. Third evaluation: supervised machine learning to predict group membership from pain-intensity ratings

After rotation and rearrangement of the data matrix, a classification problem was established by the creation of the 2 matrices described above, comprising (1) the 68 × 6-sized subject/CO2 concentration vs VAS-rating acquisition matrix and (2) the 51 × 8-sized subject/experimental condition vs CO2 stimulus strength rating matrix. Subsequently, the data were submitted to supervise machine-learning analyses. The idea behind these analyses was to train an artificial intelligence algorithm, based on different methods of machine learning, to learn the association of the CO2 intensity ratings with (1) the physical strength (CO2 stimulus intensity) of the rated nociceptive stimulus as a positive control for a successful pain experiment and with (2) the TMS condition. The goal was to find out whether this method performed better than guessing the class association (CO2 stimulus intensity or TMS condition) and, on the basis of this result, to decide whether the pain ratings provide information relevant to the class association, or if they are not related to the underlying stimulus strength or experimental condition.

In the present analysis, the mapping of the input space to the output space was performed using different methods of supervised machine learning: (1) random forests,[11] (2) k-nearest neighbors (kNNs),[16] (3) support vector machines,[15] and (4) naive Bayesian[7] classifiers. The goal was to provide an internal validation of the results, but not to compare the performances of different machine-learning methods. The choice of supervised methods covers a variety of machine-learning classifiers previously used in pain research[58] such as prototype-based (eg, kNN), collective decision-based (eg, random forests), and neuronal network-based (support vector machines) classifiers, with the addition of a classical method based on a naive Bayesian classifier.

Random forests creates sets of different, uncorrelated, and often very simple decision trees[11] with conditions on features as vertices and classes as leaves. The splits of the features are random and the classifier relates on the majority vote for class membership provided by a large number of decision trees. In the present analysis, 1000 decision trees were built containing sqrt(d) features as the standard setting implemented in the R library “randomForest” (https://cran.r-project.org/package=randomForest).[52] The number of trees was heuristically based on visual analysis of the relationship between the number of decision trees and the classification accuracy, which indicated that beyond 100 trees, the classification-balanced accuracy remained stable and a larger number merely consumed available computation time.

The kNN classification[16] is a nonparametric method that belongs to the most frequently used algorithms in data science, although it is one of the basic methods in machine learning. During kNN model building, the entire labeled training data set is stored while a test case is placed in the feature space in the vicinity of the test cases at the smallest high-dimensional distance. The test case receives the class label according to the majority vote of the class labels of the k training cases in its vicinity. In the present implementation, the size of k was established in resampling experiments with k set at 3 or 5. Even numbers of k intuitively make a majority vote on which the class assignment is based difficult when one of the nearest neighbors belongs to class 1 and the other to class 2. We tested k = 3 and 5 as these values are often used and constitute the default in various implementations of kNN. At k = 3 and Euclidean distance, the best classification accuracy was observed and more sophisticated implementations of nearest neighbor-based class assignment such as weighting or the use of kernel of different shapes were tested but did not provide any improvements regarding the basic version. These calculations were performed using the R software package “class” (https://cran.r-project.org/package=class).[87]

Support vector machines are supervised learning methods that classify data mainly based on geometrical and statistical approaches used for finding an optimum decision surface (hyperplane) that can separate the data points of one class from those belonging to another class in the high-dimensional feature space.[15] Using a kernel function, the hyperplane is frequently selected in a way to obtain a tradeoff between minimizing the misclassification rate and maximizing the distance of the plane to the nearest properly classified data point. In the present analysis, a Gaussian kernel with a radial basis was used. The analyses were performed with the R library “kernlab” (https://cran.r-project.org/package=kernlab).[46]

Bayesian classifiers were used, which provide the probability of a data point belonging to a specific class calculated by application of the Bayes theorem.[7] In naive Bayesian classifiers, the oversimplified assumption is included that all features are conditionally independent of each other, which is a widely used technique to assign class labels to the samples from the available set of features, describing a special case of the more general Bayesian network model. The calculations were performed with the R package “e1071” (Meyer D, https://cran.r-project.org/package=e1071).

The machine-learning methods were applied to the 2 original data matrices as described above. In addition, for the pain intensity vs CO2 concentration condition, positive and negative control data matrices were created by rearranging the original VAS ratings with respect to the associated CO2 stimulus class information. Specifically, a positive control data set was obtained by sorting the original VAS rating information in descending order such that the highest ratings were artificially associated with the 80% vol/vol CO2 condition while the lowest ratings became associated with the 35% vol/vol CO2 stimuli. The expectation was that with this data scenario, the association of the CO2 stimulus classes could be almost perfectly obtained by all machine-learning methods. In addition, a negative control data set was obtained by random permutation of the VAS data such that the ratings lost correlation with the CO2 stimulus concentrations. The expectation was that, during machine learning, the association with the stimulus classes was not better than guessing and should be consistently outperformed by the mapping of the true VAS ratings to the stimulus classes.

The concept of training an artificial intelligence algorithm with information to enable it to correctly associate an individual with a stimulus class or experimental condition requires measures against overfitting,[67] which are usually based on splitting the data set into a training subset that is provided to the artificial intelligence algorithm during the learning phase and a test subset that is not accessible to the artificial intelligence algorithm during learning but provided when the trained algorithm is used for classification; usually, this procedure is repeated several times in a resampling design.[67] Therefore, in all data sets, the classifiers were trained at training data subsets comprising 2/3 of the data, and subsequently, their performance was estimated on the test data subset consisting of the remaining 1/3 of the data. This was repeated in 1000 cross-validation runs using Monte-Carlo[30] resampling and random splits of the original data set into new training and test data subsets, using the R library “sampling” (https://cran.r-project.org/package=sampling).[81] Classifier performance was primarily assessed as balanced accuracy,[13,86] which is the mean of prediction sensitivity and specificity for each olfactory diagnosis and reflects the average of correctly classified cases proportional to the number of class members in class. Furthermore, secondary measures of average classification performance across olfactory diagnoses included test sensitivity and specificity, and negative and positive predictive values calculated using standard equations.[1,2]

3. Results

All 18 subjects finished the study without experiencing any noticeable side effect of cTBS.

3.1. Pain-related brain activations

From the 18 healthy right-handed volunteers, only 16 (mean age ± SD, SD: 23.8 ± 2.3 years; 6 men) were included in the analysis, due to insufficient data about stimulus timing in relation to the image acquisition for 2 subjects. In the following, the results of the 2 × 3 ANOVA design, which was applied to assess the influence of the factors “measurement” (baseline, post-cTBS) and “condition” (cTBS treatment: M1 verum, S2 verum, sham) on the linear relationship between stimulus strength and associated brain activations addressed in the first level, will be presented.

3.1.1. Brain activations associated with different stimulus strengths

Areas showing a linear relationship between brain activation and stimulus strength were identified through SPM first-level analysis with contrast [3 1 1 3] for the nociceptive stimuli at CO2 concentrations of 35, 50, 65, and 80% vol/vol. Second-level group analysis of brain activation vs stimulus strength revealed 10 clusters (Table 1). The global maximum activation was located at MNI coordinates x = 38, y = −4, z = 10 mm (Fig. 2) in the largest cluster (1668 voxels), located at the right insula, the right Rolandic operculum (corresponding to S2) and, in addition, the right supramarginal gyrus. Results were statistically significant at the FWE-corrected level (F > 55, P < 0.00001). The second largest cluster comprised 1332 voxels (FWE-corrected significant level of F > 59, P < 0.00001) and was located in the left insular cortex (MNI coordinates x = −36, y = −8, z = 12 mm).

Table 1

Clusters of brain regions showing a linear relationship between brain activation and stimulus strength, modeled at the SPM first-level analysis as contrast [-3 -1 1 3] for the nociceptive stimuli at CO2 concentrations of 35, 50, 65, and 80% vol/vol, respectively.

Figure 2.

Brain regions showing a linear relationship between activation and stimulus strength, modeled at the SPM first-level analysis as contrast [3 1 1 3] for the nociceptive stimuli at CO2 concentrations of 35, 50, 65, and 80% vol/vol, respectively (glass brain, left). The activations are superimposed on axial slices of the canonical MR template implemented in SPM12 (right). The significance at voxel level is color coded from dark red to yellow with increasing F values (Table 1). Activations are shown at a threshold of P < 0.05 (family-wise error-corrected). The figure was created using the SPM12 Matlab toolbox (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/; Wellcome Department of Imaging Neuroscience, London, United Kingdom[24,92]) and the xjView Matlab toolbox (http://www.alivelearn.net/xjview).

3.1.2. Influence of the transcranial magnetic stimulation condition on the relationship between brain activation and stimulus strength

The influence of the cTBS condition on the relationship between brain activation and stimulus strength was modeled using a 2 × 3 ANOVA design with the factors “measurement” (baseline, post-cTBS) and “condition” (cTBS condition: M1 verum, S2 verum, and sham). This analysis revealed 7 clusters (Table 2). The global maximum activation was located at MNI coordinates x = 56, y = −6, z = 14 mm (Fig. 3) in the second largest cluster (3314 voxels), located at the right Rolandic operculum (corresponding to S2) and the right insula. Results were statistically significant at the FWE-corrected level (F > 55, P < 0.00001). The largest cluster comprised 3659 voxels (MNI coordinates x = −36, y = −10, z = 14 mm) located in the left insular and S2 areas, and in addition, in the postcentral gyrus corresponding to S1. Subsequent t-contrasts indicated that significant reductions of the steepness of the brain activation vs stimulus strength (CO2 concentration) relationship modeled at level 1 were observed after the M1 verum condition, whereas cTBS targeted at S2 had no effect (empty glass brains; Fig. 4).

Table 2

Clusters of brain regions reflecting the different influences of the cTBS condition on the linear relationship between brain activation and stimulus strength, modeled at the SPM first-level analysis as contrast [-3 -1 1 3] for the nociceptive stimuli at CO2 concentrations of 35, 50, 65, and 80% vol/vol.

Figure 3.

Brain regions reflecting the different influences of the cTBS condition on the linear relationship between brain activation and stimulus strength, modeled at the SPM first-level analysis as contrast [3 1 1 3] for the nociceptive stimuli at CO2 concentrations of 35, 50, 65, and 80% vol/vol. The results of a factorial 2 × 3 ANOVA testing the influence of the factors “measurement” (baseline, post-cTBS) and “condition” (cTBS condition: M1 verum, S2 verum, and sham) on the linear relationship between stimulus strength and associated brain activations addressed in the first-level analysis are shown as a glass brain representation (left) and superimposed on axial slices of the canonical MR template implemented in SPM12 (right). The significance at voxel level is color coded from dark red to yellow with increasing F values (Table 2). Activations are shown at a threshold of P < 0.05 (family-wise error-corrected). The figure was created using the SPM12 Matlab toolbox (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/; Wellcome Department of Imaging Neuroscience, London, United Kingdom[24,92]) and the xjView Matlab toolbox (http://www.alivelearn.net/xjview). ANOVA, analysis of variance; cTBS, continuous theta-burst stimulation.

Figure 4.

Brain regions where the cTBS treatment reduced the steepness of the linear relationship between brain activation and stimulus strength, for the M1 verum and the S2 verum conditions against the sham TBS condition, modeled as SPM t-contrasts [1 0 1 1 0 1] and [0 1 1 0 1 1], respectively, in a design matrix–ordered measurement1/condition1, measurement1/condition2, measurement1/condition3, measurement2/condition1, measurement2/condition2, and measurement2/condition3, according to the factors defined as “measurement” (baseline, post-cTBS) and “condition” (cTBS condition: M1 verum, S2 verum, and sham). Although a reduction of the steepness of this relationship was observed during the M1 verum experimental condition, no effect was obtained during the S2 verum condition. The figure was created using the SPM12 Matlab toolbox (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/; Wellcome Department of Imaging Neuroscience, London, United Kingdom[24,92]). cTBS, continuous theta-burst stimulation.

3.2. Pain-intensity perceptions

Data were lost because of a technical error (failure of recording) in one subject (#2). Therefore, the analyzed cohort comprised 17 healthy right-handed volunteers (mean age ± SD, SD: 24.2 ± 2.7 years; 7 men).

3.2.1. Analysis of variance

The linear relationship between the pain VAS ratings and the strength of the CO2 stimuli did not change after cTBS treatment. This was reflected in the results of the rm-ANOVA analysis of the slopes of this relationship, where neither statistically significant main effects of the rm-ANOVA factors “measurement” or “condition” were observed, nor was the important interaction between these factors statistically significant (all P-values > 0.5; for details see Table 3).

Table 3

Results of the analyses of variance for repeated-measures (rm-ANOVA) of the slopes of the linear relationship between pain ratings and stimulus strength and of the VAS ratings of the nociceptive CO2 stimuli (Fig. 4), specifying the number of degrees of freedom (df) and the F and P-values.

The pain-intensity ratings (individual medians per CO2 concentration, experimental condition, and measurement) on the VAS (range: 0-100) were between 0 and 88.5. Pain-intensity ratings increased with increasing CO2 concentration (Fig. 5). Although this was reflected in a statistically significant effect of the factor “CO2 concentration” in the rm-ANOVA (Table 3), the ratings neither differed between the baseline and the poststimulation sessions (nonsignificant effect of the rm-ANOVA factor “measurement”), nor did the cTBS condition have any effect on the pain intensity of the CO2 stimuli (nonsignificant effect of the rm-ANOVA factor “condition” and no significant interactions “condition” by “measurement” by “CO2 concentration”; for details, see Table 3).

Figure 5.

Intensity ratings of the CO2 stimuli on a VAS ranging from 0 to 100. For each strength (35, 50, 65, or 80% vol/vol), 20 stimuli were applied during baseline sessions, ie, before cTBS application, and after cTBS application at 3 different experimental conditions, ie, M1 verum, S2 verum, and sham (placebo). The Spaghetti plots show the individual median VAS ratings per stimulus strength (dots), connected by a straight line and color-coded for each individual subject. The quartiles and medians (solid horizontal line within the box) were used to construct the overlaid “box and whisker” plots. The whiskers add 1.5 times the interquartile range (IQR) to the 75th percentile or subtract 1.5 times the IQR from the 25th percentile and are expected to include 99.3% of the data if normally distributed. A linear regression line and its confidence interval of estimate are overlaid showing the group trend of the pain intensity vs stimulus strength relationship. The figure has been created using the library “ggplot2” (https://cran.r-project.org/package=ggplot2)[90] with the R software package (version 3.4.3 for Linux; http://CRAN.R-project.org/).[20] cTBS, continuous theta-burst stimulation; VAS, visual analog scale.

3.2.2. Unsupervised data exploration for subgroup identification

Unsupervised data analysis was used to explore whether the data displayed a systematic structure that reflected the experimental condition or hinted at subgroups of subjects in the sense of responder/nonresponder.

Cluster analysis using the progeny algorithm resulted in k = 2 clusters as the optimum number to describe the pattern of pain ratings when rearranged for CO2 concentrations (Supplementary Figure 1, available at http://links.lww.com/PAIN/A660), whereas k = 4 clusters were preferred to describe the pattern of pain ratings when rearranged for cTBS conditions. These numbers provided the most stable clustering, starting with an initial clustering of the full data set, based on the maxima observed with the cluster stability score criterion proposed for the progeny clustering algorithm.[36] The clusters identified in the pain ratings when rearranged for CO2 concentrations agreed significantly with the CO2 concentrations as indicated by significant χ2 test (χ2 = 42.118, P = 3.788 × 10−9; supplementary Figure 1, available at http://links.lww.com/PAIN/A660). This persisted when testing out k = 4 clusters, for which the heat plot and a cluster dendrogram provided support, ie, the χ2 test remained significant in that 4-cluster scenario (χ2 = 67.676, P = 4.34 × 10−11). By contrast, the clusters identified in the pain ratings when rearranged for the cTBS condition did not reflect the experimental treatment (χ2 = 3.631, P = 0.7289; supplementary Figure 1, available at http://links.lww.com/PAIN/A660).

The steepness of the linear relationship between the pain rating and the physical stimulus strength was expected to decrease after cTBS-mediated inhibition of S2 processing.[71] Therefore, the slopes of this relationship were evaluated after cTBS application and at baseline, and their differences were analyzed for group structures across the whole study. The expectation was a group separation towards larger negative differences during the S2 verum condition. In addition, possible subgroups could hint at the necessity to split the cohort into responders and nonresponders in subsequent analyses. However, visual inspection of the probability density distribution (PDF) of the slope differences suggested merely a unimodal normal distribution with a few outliers (supplementary Figure 2, available at http://links.lww.com/PAIN/A660). The outliers probably had the effect that statistical evaluation of goodness of fit favored a bimodal distribution, ie, a Gaussian mixture model using M = 2 modes as indicated by a significant likelihood ratio test between M = 1 and M = 2 Gaussian modes (difference in minus 2 times log likelihood: −19.03, P = 2.74 × 10−8). No more significant improvement of the fit was obtained when a further Gaussian was added, based on likelihood ratio tests (difference in minus 2 times log likelihood: −0.996, P = 0.574). However, the slope differences were evenly distributed between cTBS experimental conditions across the whole data range (colored dots in supplementary Figure 2, available at http://links.lww.com/PAIN/A660) discouraging smaller differences observed during the S2 verum conditions or the presence of a responder subgroup.

3.2.3. Supervised machine learning to predict group membership from pain-intensity ratings

Supervised machine learning applied in cross-validation experiments using 1000 bootstrap random resamplings of 2/3 (new training) vs 1/3 (new test) of the data provided the consistent observation that, when using the true VAS ratings, the class assignment to the CO2 concentration of the rated stimuli was better than guessing (supplementary Figure 3, available at http://links.lww.com/PAIN/A660), which established the relationship between pain increasing with increasing stimulus strength on the basis of quantitate standard test performance measures (supplementary Figure 4, available at http://links.lww.com/PAIN/A660). Thus, it can be concluded that the VAS ratings yield information about the strength of the pain-generating stimulus. This finding was further strengthened by the observation that, when permuting the VAS ratings, ie, breaking their relationship to the strength of the generating stimulus, all machine-learned classifiers were not better than guessing, ie, were unable to assign the stimulus strength based on the VAS information provided. The assignment was not perfect as indicated by the almost 100% classification accuracy achieved when the VAS ratings were sorted, so that higher ratings became associated with stronger stimuli. The latter experiment served as positive control to verify that a less than perfect classification was not due to poor implementation of the machine-learning methods.

In contrast to the observation that the VAS ratings allowed conclusions about the strength of the pain-generating stimulus, associating the ratings with the cTBS-related experimental conditions failed. Specifically, all machine-learned classifiers did not perform better than guessing when training them with VAS rating information acquired during the different experimental conditions, regardless of whether the ratings before and after cTBS application were used or their differences (supplementary Table 1 and supplementary Figure 3, available at http://links.lww.com/PAIN/A660).

4. Discussion

In this study, cTBS application to the M1 was followed by a significant impairment of the relationship between the strength of nociceptive stimuli and the resulting brain activation. This effect was mainly located in the insular cortex and the secondary somatosensory area S2 (Fig. 3), which has previously been identified as being crucially involved in the processing of pain stimuli.[71] Therefore, the results of this study clearly suggest a cTBS effect on the central processing of trigeminal nociceptive inputs. This is in line with previous reports of occasionally successful pain modulations through targeting M1.[4,6,18,23,35,48-51,62,68,75,77]

However, an analogous effect of the cTBS application on S2 was not observed. Hence, although the study was successful from a technical point of view, it did not support the hypothesis that inhibition of the function of S2 by setting a “virtual lesion” should result in a flattening of the relationship between pain-intensity and physical stimulus strength. This is despite S2 was proposed as a cortical target for orofacial neuropathic pain where its stimulation produced more pronounced analgesic effects than that of S1/M1[53] and also contrasts with findings that rTMS over S2 produced hypoalgesia to heat and cold stimuli; however, the latter with a restriction to male volunteers[85] or with observations that TMS over S2-disrupted pain-intensity coding,[54] which in turn agrees with own previous observations.[71] Therefore, this study cannot discourage the inclusion of S2 as a TMS target in future studies on the modulation of pain through this method.

The present negative findings could be due to insufficient inhibition of S2 by cTBS. It should be noted that subjects were selected according to the electromyographic observations after M1 deactivation using cTBS, which led to the report of a responder subgroup[33] from which this study participants were exclusively selected. Although this establishes the principal functioning of the experimental setup and proves the influence of cTBS on M1, subject selection was based on a motor effect generated at M1. This favoured M1 over S2, and in fact, stimulation of M1 was successful. However, the effect addressed in this study was a sensory change. We are not aware of a test scenario that would allow for verifying a positive effect on the sensation of pain without performing the whole study already during the verification phase. Therefore, although having established a general success of the chosen cTBS protocol, insufficient modulation of the function of S2 remains as a possible cause for the absent changes in pain perceptions when applying cTBS on that region. There may be even a possibility that the selection method had introduced a bias for positive findings with M1, although in that case, it remains unclear how this should have led to the enrolment of subjects particularly insensitive to cTBS on S2. Another possibility as a directed bias is a variability of the cTBS responses in S2 independent of those in M1. Indeed, it cannot be excluded that the changes in sensory perceptions after cTBS at S2 showed similarly high variability as observed in the changes in motor responses after cTBS at M1; however, without that, the observed responder/nonresponder subgroups necessarily included the same subjects. Hence, although cTBS on M1 might have suppressive effects as intended, cTBS on S2 may have exerted suppressive effects only in some participants while no effects or opposite effects in others. This again emphasizes the need of a responder selection criterion for sensory cTBS effects. Exploring the relationship between changes in intensity ratings and changes in the magnitude of the blood oxygenation level–dependent response after cTBS on S2 might guide the design of further studies addressing systematically the response rates to S2 cTBS. However, in the present sample, the relationship between stimulus strength and rating seemed unaffected by cTBS on S2 except for a single subject (green line in Fig. 5). Thus, the present sample size seems insufficient to address possible subgroups at a reliable statistical power.

The finding that cTBS application on S2 did not provide the hypothesized reduction in the processing of nociceptive input might be due to a more complex architecture of stimulus intensity processing that might not be adequately addressed by applying cTBS exclusively to this brain area. In a previous experiment using identical CO2 stimuli of different strengths, activations associated with the nociceptive stimulation were observed in most regions of the so-called “pain matrix”[5,14,61]; however, correlations between brain activation and stimulus intensity were found in the posterior insula, the primary and secondary somatosensory cortex, the amygdala, and the middle cingulate cortex. This distribute network may have contributed to the failure to observe significant reductions in stimulus intensity processing when impeding only a single component of this network. Moreover, in independent assessments, alternative sites of stimulus intensity coding were proposed. For example, in experiments using fMRI[8] or intracerebral EEG[25] recordings, stimuli perceived as painful were observed to elicit a correlated response in the insular cortex rather than in S2, while the responses in S2 increased with stimulus intensity only below the pain threshold and reached a plateau at suprathreshold levels. This contrasts with the clear findings of an S2 localization of stimulus intensity coding, on which the present hypothesis was based[54,71]; however, these observations offer an explanation of the lack of cTBS effects when targeted on S2 in the present experiments.

Among brain regions where cTBS over M1 altered the linear relationship between stimulus-related activations and stimulus strength was the thalamus (Table 2). Activation of M1 is not only a pain epiphenomenon such as related to pain-evoked movement, but a specific local response to nociceptive stimuli shown by means of intracortical evoked potentials to non-noxious electrical and noxious laser heat stimuli recorded from epileptic patients.[26] Direct connections between M1 and the thalamus have been shown in owl monkeys,[79] where each of the major somatotopic divisions of M1 was observed to be connected in a somatotopically organized fashion with an anteroposteriorly elongated territory within the ventrolateral thalamic complex, and most projections were reciprocal.[79] Hence, the effects located in the thalamus observed in this study probably directly show a part of the mechanism of pain modulation by cTBS over M1. Moreover, the thalamic effects could also be generated through corticospinal circuits. Specifically, for movement-relevant circuits, it has been shown recently that corticospinal fibers originating in the mouse motor cortex directly synapse onto spinal premotor interneurons.[82] Moreover, activation in the motor cortex (M2) induced firing of layer V neurons in S1 while blocking this pathway caused deficits in sensory perceptions in mice,[60] and corticocortical connections between M1 and S1 also exist in humans.[45] Furthermore, a direct connection between M1 and the spinal level is supported by observations that greater gray-matter volume in M1 after spinal cord injury was associated with amount of neuropathic pain,[44] and corticospinal inputs have been found in cTBS experiments in humans.[27] Thus, present observations are compatible with evidence about a more complex involvement of the motor cortex in the modulation of pain including thalamic and corticospinal pathways.

Considering this possible involvement of S1 or M1, the vicinity of these brain regions to S2 arises the possibility that cTBS delivered over S2 may have also modulated the face representation within those regions. To contemplate this further, previously observed MNI coordinates available from studies with an explicit focus on the localization of brain regions activated by the CO2 stimulus were consulted. Specifically, after right-sided intranasal stimulation, S1 was observed at MNI coordinates of x, y, z = −60, −21, 33 mm, and after left-sided intranasal stimulation, S2 was observed at 57, −6, 24 mm.[9] To estimate the distance, the Signed Differential Mapping neuroimaging software library online coordinate utilities (Radua J, https://www.sdmproject.com/utilities/?show=Coordinates) was used. The distance between these two coordinates was 118.3 mm, and when switching S2 to the same side as S1, it was 17.6 mm. Alternatively, the closest distance between activations assigned to either S1 (−39, −27, 54 mm) or S2 (−45, −18, 21 mm) in 71 was 34.73 mm. In both variants, the focal effects of the cTBs impulses of approximately 5-mm diameter are exceeded. This remains true when calculating the distance between S1 and the line linking the center of the TMS coil to S2 target. Specifically, the coordinates of the S2 target were at x, y, z = 48, 2, 10 mm (see supplementary Figure 4, available at http://links.lww.com/PAIN/A660), which indicates that the coil held over the skull was placed at approximate coordinates of x, y, z = 66, 2, 20 mm, which is when switching to the same side at a distance of 27.1 mm from the S1 coordinates taken from 9 or at 52.2 mm from the S1 coordinates taken from 71.

Similar data for M1 were not available because the face representation was not experimentally addressed as in 85 where, however, the coordinates were not reported. Therefore, a costimulation of S1 with cTBS over M1 cannot be excluded. Regardless of scattering effects of cTBS across regions, coinfluences of cTBS over S2 through above-mentioned interregional connections are possible and need to be considered when interpreting present results. Indeed, this might have played a role considering the small cTBS effects on the activation vs concentration relationship observed at the stimulated sites, according to an explorative SPM analysis of the effect sizes using the SPM toolbox rfxPlot[28] (supplementary Figure 4, available at http://links.lww.com/PAIN/A660).

Although cTBS applied to M1 clearly altered the brain processing of nociceptive stimuli with different strength in a way that reduced the slope of the linear relationship between brain activation and stimulus strength, the cerebral modifications were not accompanied by a similar change in the perceptions of experimentally induced trigeminal nociceptive stimuli. By contrast, although several different methods of analyzing the pain rating data were used, comprising both the fMRI factorial analysis of the slope of the above relationship implemented as an rm-ANOVA, and several different supervised methods including machine learning, results did not support the hypothesis of a cTBS effect on the perception of trigeminal nociceptive input. The complex analysis of the VAS data was performed to prove that the psychophysical data provided sufficient resolution of different stimulus strengths to observe this effect. This was verified by the consistent success of several different bioinformatics methods to show that from the VAS data, the strength of the originating pain stimulus could be deduced. In the psychophysical data, both the M1 and the S2 cTBS treatment were not reflected, contrasting with the clear coding of stimulus strength by these data. Hence, while modulating the central processing of nociceptive input, cTBS failed to produce a subjectively relevant change in pain perception, indicating that the method in the present implementation is still unsuitable for clinical application. It should be noted that the results were obtained on healthy subjects, rather than on pain patients. However, this should not constitute a problem because it has been shown that the results obtained with the present pain model serve as satisfactory predictors of the clinical efficacy of analgesics.[57,70]

Conflict of interest statement

The authors have no conflict of interest to declare.