Imaging the neural substrate of trigeminal neuralgia pain using deep learning

1. Introduction

The 3rd edition of the International Classification of Headache Disorders (ICHD-3) classifies trigeminal neuralgia (TN) into idiopathic, classical, or secondary trigeminal neuralgia diagnoses (Arnold, 2018). Idiopathic TN, by definition, has no known etiology, while classical TN is a diagnosis given when there is verification (e.g., visualization during surgery and/or via neuroimaging) of blood vessel contact against the root of the trigeminal nerve on the ipsilateral side of the pain complaint. Secondary TN is caused by an identifiable pathological condition, including multiple sclerosis or tumor impingement on the trigeminal root trunk, most commonly at the cerebellopontine angle. For classical and idiopathic TN, recent diagnostic subdivisions also include either (1) a purely paroxysmal quality or (2) having intermittent paroxysmal bouts with a more continuous background of burning and/or aching pain (Arnold, 2018; Bendtsen et al., 2020). After diagnosis, first line treatments of TN involve anti-convulsant medications such as carbamazepine and oxcarbazepine, which often have intolerable side effects and become progressively less effective with time (Prasad and Galetta, 2009). Following failure of medication trials, patients usually escalate treatment to invasive surgical procedures. The most common surgical procedure is called microvascular decompression (MVD), which involves moving an impinging blood vessel (e.g., superior cerebellar artery) off the trigeminal nerve and placing a Teflon pad to maintain this neurovascular separation. While surgical procedures provide lasting relief for some patients, for others, the relief can be short-lived, and pain returns after a few months to a few years. According to available data, approximately 4% of patients per year experience recurrence of TN pain after MVD (Cheng et al., 2019), and over a period of 10–20 years, the recurrence rate may exceed 10% (Barker et al., 1996). Moreover, while classical TN is proposed to have a neurovascular insult as the etiology, many cases of TN do not show this impingement; conversely, many patients have impingement, but no pain. The lack of an etiology and the unpredictable response to treatment suggest that our understanding of TN remains inadequate and a range of peripheral and especially central mechanisms underlying TN pain remain to be better understood (Gambeta et al., 2020).

Functional magnetic resonance imaging (fMRI) remains the main neuroimaging technique for investigating the central mechanisms of pain. In particular, resting-state fMRI (rs-fMRI), which measures blood-oxygen-level-dependent (BOLD) signals while the patient is not engaged in any systematic thought or activity, has been applied extensively in TN (Dou et al., 2016; Yuan et al., 2018). Yuan et al. (2018) found that, compared to age- and sex-matched healthy control subjects, TN patients exhibited significantly increased regional homogeneity (ReHo) and signal amplitude in several brain regions, including the posterior lobe of the cerebellum, anterior cingulate cortex, middle temporal gyrus, temporal lobe, putamen, occipital lobe, limbic lobe, precuneus, and the medial and superior frontal gyrus, and decreased ReHo in insula. They suggested that these abnormal activities played a role in the development and maintenance of chronic TN pain. Wang et al. (2015) found increased ReHo in the inferior temporal gyrus, thalamus, inferior parietal lobule, precentral and postcentral gyri and decreased ReHo in the amygdala, parahippocampal, and cerebellum in TN patients. Yan et al. (2019) detected that TN patients had decreased ReHo in the left middle temporal gyrus, superior parietal lobule, and precentral gyrus and increased ReHo in the thalamus in the resting state studies. Moving beyond regional analysis, at the network level, Zhu et al. (2020) applied functional connectivity (FC) methods to show that TN patients exhibited significantly higher degree centrality values in the right lingual gyrus, right postcentral gyrus, left paracentral lobule, and bilateral inferior cerebellum. They proposed that these changes reflect the adaptation of the cerebral cortex to frequent pain attacks over a long period of time. In addition to resting-state fMRI, Moisset et al. (2011) evoked pain in TN patients with stimulation of the cutaneous trigger zone and found increased activity in postcentral and precentral cortex, contralateral supplementary motor area, thalamus, anterior and posterior insula, prefrontal cortex, putamen, ipsilateral midcingulate cortex, hippocampus/parahippocampal area and cerebellum. Table 1 summarizes these neuroimaging studies, highlighting the brain regions identified as having abnormal neural activities in TN subjects. Note that, with the exception of the cerebellum, there is a relative lack of consistency for identified brain structures across the studies, suggesting additional studies are needed for understanding and identifying brain structures that are functionally important in the generation and maintenance of TN pain.

Chronic pain can cause changes in many structures in the brain (Rodriguez-Raecke et al., 2013). While the neural activity in a brain area may appear to be altered in rs-fMRI studies, it may not imply that the area is directly involved in pain generation or perception. Additionally, stimulus-evoked pain may differ fundamentally from spontaneous or continuous pain experienced by the patient with TN (Wan et al., 2018). A neuroimaging technique that directly addresses the neural substrate of the naturally occurring pain is percept-related fMRI (Apkarian et al., 2001). In this technique the patients indicate their moment-to-moment pain levels while their brain activities are being recorded (Wilcox et al., 2015). Correlating brain activity and pain level fluctuations one obtains information on brain structures that directly underlie pain generation and perception. Kwan et al. (2005) applied this percept-related fMRI technology to patients with irritable bowel syndrome (IBS) and revealed abnormal urge- and pain-related forebrain activity during rectal distension. Baliki et al. (2006) applied the method to patients with chronic back pain (CBP) and found that the insular region was active when pain level increases. This knowledge is not only important for understanding the neural mechanisms of the specific pain condition it also has the potential to provide biomarkers for evaluating disease progression and treatment effectiveness. To date, however, no studies have applied the percept-related fMRI approach to investigate TN pain.

Although the correlation analysis adopted in previous percept-related fMRI studies of pain is intuitive and easy to apply (Apkarian et al., 2001), it has limitations, including (1) it can only detect linear relationships and (2) brain activations are often reported without applying multiple comparison corrections, raising robustness concerns. The emergence of AI-inspired deep learning methods such as convolution neural networks (CNN) and graph convolution neural networks (GCNN) offers an avenue to overcome these limitations. These deep learning methods differ from traditional machine learning techniques such as support vector machine (SVM) and logistic regression (LR) in that the hidden layers in these models are capable of encoding and utilizing more complex features of the data to provide more accurate predictions of experimental conditions and a deeper understanding of the data (Van Der Miesen et al., 2019). Santana et al. (2019) applied CNNs to resting-state fMRI data to distinguish chronic pain patients from pain-free controls and demonstrated that CNNs could achieve higher classification accuracy compared to other machine learning models. Li et al. (2019) applied GCNN to task-fMRI to patients suffering from Autism Spectrum Disorder (ASD) and discovered regions and networks in the brain that can serve as biomarkers to distinguish patents from controls. It is expected that by combing AI-inspired methods with the more conventional correlation analysis and by requiring additionally that the results be consistent across methods one will have enhanced ability to obtain new and robust insights into the neural basis of TN pain.

We recorded fMRI data while TN patients (1) rated their spontaneous pain levels in the pain tracking session (percept-related fMRI) and (2) rested in the resting state session (resting-state fMRI). The data were first analyzed using the conventional correlation method and then subjected to the analysis by CNN and GCNN. A number of validation analyses were also carried out. The goal was to identify a set of brain regions, called signature centers, that are robustly activated during TN pain.

2. Materials and methods

The study was approved by the WCG Institutional Review Board (IRB). Patients were recruited through the clinical care population within the University of Florida Health System and from referrals provided by the Facial Pain Research Foundation (FPRF). Screening was done either in person or via a phone call. Eligible subjects were consented and completed a study packet that included a Health History Questionnaire, Oregon Health Science University (OHSU) Trigeminal Neuralgia – Diagnostic Questionnaire, Beck Depression Inventory-II (BDI-2), Pain Anxiety Screening Scale (PASS), and Pain Catastrophizing Scale (PCS). The study coordinator read a prepared standard script explaining the study procedures. A focused medical history, a trigeminal cranial nerve exam, and a physical exam was completed by a trained clinical fellow and the PI (JN). Vital signs (blood pressure, temperature, and pulse) were also recorded prior to scanning.

2.1. Pain ratings

Subjects indicated on a 100mm visual analog scale (VAS) anchored on the left with “no pain sensation” and on the right with “most intense pain sensation imaginable” their daily experienced pain (past month). This rating is henceforth referred to as “usual pain”. Subjects also used the VAS to rate their “current pain” just prior to entering the MRI scanner. During the scanning procedure, subjects completed continuous pain tracking (see Experimental Paradigm below, Figure 1B) in the pain tracking session where they could visualize a computer screen in the scanner via mirrors and rate their pain levels in real time using a tracking ball. Besides being used in neuroimaging analysis, these pain ratings were also applied to compare different subgroups of patients (e.g., male vs female), using t-test.

2.2. Participants

2.2.1. Inclusion criteria

• Male and female subjects.

• Age 18–75 years old. This represents the age range of many chronic pain patients, including patients with TN.

• American Society of Anesthesiologists (ASA) status 1, 2, or 3, deemed in good general health.

• Only subjects with reported average usual pain of moderate to severe (VAS of 30–100 mm) at the time of screening were included.

• Subjects were diagnosed with trigeminal neuralgia, per the International Headache Society (IHS) Disorders criteria (Arnold, 2018), as having:

A. Recurrent paroxysms of unilateral facial pain in the distribution(s) of one or more divisions of the trigeminal nerve, with no radiation beyond, and fulfilling criteria B and C.

B. Pain has all of the following characteristics:

1. lasting from a fraction of a second to two minutes

2. severe intensity

3. electric shock-like, shooting, stabbing or sharp in quality

C. Precipitated by innocuous stimuli within the affected trigeminal distribution.

D. Not better accounted for by another ICHD-3 diagnosis.

• Subjects included having symptoms as purely paroxysmal pain or having intermittent paroxysmal bouts with a more continuous background of burning and/or aching pain.

Note that during the screening process, individuals reported being diagnosed with TN by their physician would have this diagnosis verified by the study staff. Additionally, subjects completed the OHSU TN – Diagnostic Questionnaire as additional verification of TN diagnosis.

2.2.2. Exclusion criteria

• Patients diagnosed with trigeminal neuralgia attributed to space-occupying lesion (ICHD-3 code: 13.1.1.2.2) or other cause (e.g., multiple sclerosis, ICHD-3 code: 13.1.1.2.3), painful trigeminal neuropathy (ICHD-3 code: 13.1.2), trigeminal post-herpetic neuralgia (ICHD-3 code: 13.1.2.2), trigeminal neuropathic pain (ICHD-3 code: 13.1.2.4), and idiopathic painful trigeminal neuropathy (ICHD-3 code: 13.1.2.5).

• ASA status 4–5 and Emergency operation.

• Presence of chronic disease (e.g., cardiovascular disease, liver disease, kidney disease, diabetes, etc.), other than trigeminal neuralgia.

• Pregnant females.

• No exclusions were made based on race, gender, or religion.

In total, 55 patients gave written informed consent and participated in the study (69% female, mean age ± standard deviation (SD) = 53.9 ± 14.9). Sixteen patients were rejected due to a combination of the following reasons: (1) not meeting diagnosis criteria (n = 1), (2) not completing the whole experiment (n = 11), (3) technical difficulties during fMRI recording (n = 2), and (4) excessive movements inside the scanner (n = 2). The data from the remaining 39 patients were analyzed and reported here. Of the 39 patients, diagnostic concordance between the PI and the OHSU Trigeminal Diagnostic Questionnaire was 38/39 (97.5%). For the 1/39 patient, the OSHU TN diagnosis was nervus intermedius neuralgia, which is characterized by an intermittent stabbing deep pain in the ear with associated tinnitus. However, during the examination, it was determined that this one subject met the inclusion criteria for idiopathic TN, having intermittent paroxysmal bouts with a more continuous background of burning and/or aching pain. The vital signs taken just prior to scanning were within normal limits for all subjects (data not shown) and no adverse reactions or events were reported during any of the procedures.

The clinical and neuropsychological data were analyzed to compare subgroups of patients (e.g., male vs female) using t-test.

2.3. Experimental paradigm

Subjects underwent functional, structural, and diffusion magnetic resonance imaging. As shown in Figure 1A, there are two types of functional scans: resting state (resting-state fMRI) and pain tracking (precept-related fMRI). During the resting state scan (7 min), the patients were instructed to fixate on the cross at the center of the monitor screen, stay still, and not think any systematic thought. During the pain tracking scan (10–20 min), fMRI data were acquired while the patients rated their momentary pain levels using a tracking ball. The tracking ball controlled the movement of a cursor along a straight line with 0 and 100 indicated at the two ends of the line on the monitor. An example of a pain tracking time course from one patient is shown in Figure 1B.

For the first 10 subjects, the pain tracking session was divided into two parts. For the first 10 min the subjects tracked their pain level fluctuations as described above. The second 10 min was a motion tracking session in which a marker moved on the monitor between 0 and 100 according to the subject-indicated pain level fluctuations from the previous 10 min of pain tracking. The subject was asked to move the tracking ball to track the movement of the marker. We had to discontinue the motion tracking session after the first 10 subjects because it was becoming apparent that 10 min of actual pain tracking scanning was not enough to produce sufficient data for the intended analyses. For the remaining patients the pain tracking session lasted the entire 20 min. It is worth noting that increasing the length of the overall experiment was not an option because of the burden it would place on the patient.

2.4. Data acquisition and preprocessing

Functional MRI images were collected on a 3T Philips Achieva scanner (Philips Medical Systems, the Netherlands) equipped with a 32-channel head coil. The echo-planar imaging (EPI) sequence parameters were as follows: repetition time (TR), 1.98 s; echo time (TE), 30 ms; flip angle, 80; field of view (FOV), 224 mm; slice number, 36; voxel size, 3.5 mm × 3.5 mm × 3.5 mm; matrix size, 64 × 64. In addition to the functional scans, a high-resolution anatomical T1-weighted MRI image was also acquired for each subject using the following parameters: FOV, 240 mm × 240 mm; TR, 8.0566 ms; TE, 3.686 ms; resolution, 1 mm × 1 mm; flip angle, 80^°.

Statistical parametric mapping (SPM) was used to preprocess the functional MRI data (Friston, 2003). The preprocessing steps include slice timing correction, realignment, co-register, spatial normalization, and spatial smoothing. Data segments with large head movements were removed from five subjects. For the fMRI data that were subjected to further analysis, six head motion variables (translations: x, y, z and three rotations: pitch, yaw, and roll) were regressed out, and a bandpass filter [0.01, 0.1 Hz] was applied to reduce low-frequency and high-frequency noise.

2.5. Correlation analysis

The correlation analysis has been applied in prior percept-related fMRI studies of chronic pain. In this analysis, the time course of pain ratings (see, for example, Figure 1B) was first convolved with the hemodynamic response function (HRF), and then correlated with the BOLD time course from every voxel in the brain (Davis and Moayedi, 2013). We note that this analysis is similar to the general linear model (GLM) analysis of fMRI data. In a typical GLM analysis, the experimental design matrix contains the onset times of various events in the experiment. Convolving the design matrix with the HRF function we generate the predicted BOLD response which is then compared to the actual BOLD data through a correlation analysis carried out in the framework of general linear model (Friston et al., 2006; Baliki et al., 2009; Poldrack et al., 2011). For the percept-related fMRI analysis, we convolved the pain rating function, which is akin to the design matrix, with the HRF function and compared the predicted pain-related BOLD responses with the actual BOLD data using correlation. For each subject, the correlation values were Fisher-transformed to achieve approximate normal distribution and subjected to the population level analysis. A statistical threshold was set at P < 0.05 uncorrected. To reduce the influence of possible spurious correlations, only voxels that are part of clusters of at least 10 voxels meeting this statistical threshold were considered in the brain map. Both positively correlated voxels and negatively correlated voxels were considered. Regions showing strong correlations were considered as potential signature centers of TN pain. The pipeline of the correlation analysis is illustrated in Figure 2.

2.6. CNN analysis

Increasingly, artificial intelligence (AI)-inspired techniques such as deep neural networks (DNNs) are being applied to analyze fMRI data, providing insights not possible with other techniques (Kalantar et al., 2021). We implemented a convolutional neural network (CNN) model, which is a type of DNN, to predict pain ratings from fMRI data. Given that CNNs are more adapted to predict discrete labels rather than continuous values, we divided continuous pain ratings into two categories: low pain (pain ratings <= 15) and high pain (pain ratings > 15). Here the threshold of 15 was chosen so that the number of data points in the high and low pain categories were approximately equal across the entire patient population. Numerically, low pain was given a value of zero whereas high pain a value of one. The CNN, as shown Figure 3A, consisted of ResNeXt101 with four cardinalized res-blocks with 32 independent paths within each block (Hara et al., 2018). Following Hara et al.’s study (Hara et al., 2018), the convolutional filters were modified and made three-dimensional so that they could be applied to the three-dimensional fMRI data using Conv3d. BatchNorm3d and Avgpool3d were used for the batch normalization and average pooling layer. We also modified the final global average pooling layer to fit our input and used two units in the final output layer to predict low pain and high pain. The detailed structure of the proposed ResNeXt101 model is given in Supplementary Table 5.

The 39 patients were divided into eight groups of four to five patients each. Seven groups were chosen as the training dataset and the remaining group was chosen as the testing dataset. This process was repeated eight times (eight-fold cross validation). The reported decoding accuracy and the receiver operating characteristic (ROC) curve were the averages from the eight repetitions. For model training, the weighted cross-entropy was used as the loss function and the stochastic gradient descent (SGD) as the optimizer, with a momentum of 0.9 and a weight decay of 5.0e-4. The initial learning rate was set as 1.0e-3. The number of training epochs was 50 with a batch size of two. All the hyperparameters for the proposed CNN are given in Supplementary Table 6. We trained all the models on a server containing 8 GeForce GTX 1080 Ti GPUs, and the training time for one CNN is around 1.5 days.

For statistical analysis of the CNN decoding results, a t-test was used to determine whether the average decoding accuracies were significantly higher than chance level (50%), with above chance level decoding signifying that the CNN model can predict pain levels using patterns of brain activity. Additionally, as part of the validation analysis, we investigated whether the average predicted pain level by the CNN model over the pain tracking session was consistent with the average of the reported pain level by the patient over the pain tracking session, using correlation analysis. The expectation was that these two averages should be correlated if the CNN adequately modeled the relationship between brain activity and reported pain level. Furthermore, also as part of the validation analysis, we conducted a correlation analysis to determine the relationship between the model-predicted pain during resting state and that during pain tracking, with the hypothesis being that if the model correctly captured the brain-behavior relationship, these two predicted pain levels should be significantly correlated.

Once the CNN model was shown to have the ability to predict pain ratings from fMRI data, we proceeded to identify the essential brain regions that contribute to the prediction performance. The occlusion method was used for this purpose. See Figure 3B. For each region of interest (ROI) in the Lausanne atlas (Daducci et al., 2012), the BOLD values in all of its voxels were replaced by zero and fed into the CNN model, and the decoding accuracy was recorded. The amount of decoding accuracy decrease compared to the full data decoding accuracy was taken as a measure of the importance of the ROI in model prediction. The more the prediction accuracy decreases from occluding a brain region, the more that brain region contributes to the prediction performance, and the more weight it gets in the resulting heatmap. The reported heatmaps were the average from the eight models described earlier.

2.7. GCNN analysis

It is increasingly recognized that pain processing involves multiple brain areas and their interactions and is a network phenomenon (Garland, 2012). The correlation analysis and the CNN analysis described above do not take into account functional relationships between different brain regions during pain processing. To address this problem, we implemented a GCNN approach to predict pain ratings from fMRI data. The GCNN, shown in Figure 4A, consisted of two Graph Isomorphism Network (GIN) layers (Xu et al., 2018), and one fully connected layer containing two output units for predicting low pain vs high pain. The Lausanne atlas was used to divide the brain into 129 ROIs (Daducci et al., 2012). BOLD signals within each ROI were averaged and ROI-to-ROI dynamic interactions were assessed by cross-correlation in moving windows of 30 s in duration. After the cross-correlation matrices were computed, a binarization was carried out by applying a threshold, where the values greater than the threshold were set to one and smaller than the threshold were set to zero. The binarized connectivity matrix along with the average BOLD signals from each ROI were taken as input features for the GCNN to predict the pain ratings in the middle of the same 30 s moving window. The patients were again divided into eight groups of four-five each (same eight-fold validation as CNN). We applied the weighted cross-entropy as our loss function and SGD as our optimizer with a momentum of 0.9 and a weight decay of 5.0e-4 during the model’s training phase. The initial learning rate was set as 1.0e-3. The number of training epochs was 50 with a batch size of one. The reported decoding accuracy and the ROC curve were averages of the eight-fold results (see above). All the hyperparameters for the proposed GCNN are given in Supplementary Table 6. The training time for one GCNN is around one day with the same 8 GeForce GTX 1080 Ti GPUs described above. We performed the same statistical analyses of the GCNN results as those conducted for the CNN analysis (see Section “2.6 CNN analysis”).

Once the models were shown to have the ability to predict pain ratings from fMRI data, we proceeded to apply the occlusion method to identify essential network nodes through which the functional interactions among different brain areas play an essential role in pain prediction. See Figure 4B. Each ROI along with the connections through the ROI were replaced by zero and the prediction performance decline was calculated. The degree of decline is taken to indicate the importance of the ROI in mediating pain related network processing. A heatmap was derived based on this principle.

3. Results

We recorded fMRI data while n = 39 TN patients (1) tracked their pain level fluctuations in the pain tracking session and (2) rested without thinking any systematic thought in the resting state session. Both conventional correlation analysis and AI-inspired analyses (CNN and GCNN) were applied to reveal the signature centers in the brain underlying the generation and perception of TN pain.

3.1. Correlation analysis

The data from the pain tracking session was analyzed here. The BOLD time course from each voxel was correlated with the HRF-convolved pain rating time course. A voxel was considered potentially part of the neural substrate underlying TN pain if this correlation is strong regardless of the sign of correlation. For the subject in Figure 1B, the HRF-convolved pain rating time course had a strong positive correlation with the BOLD signal from a voxel located in the precentral gyrus (Figure 5A, R = 0.55), a low correlation with the BOLD signal from a voxel located in the superior frontal gyrus (Figure 5B, R = -0.03), and a strong negative correlation with the BOLD signal from another voxel located in the superior frontal gyrus (Figure 5C, R = -0.61). Across all patients, using uncorrected P < 0.05 and a minimum cluster size of 10 voxels as the criteria, a set of brain regions showing both strong positive and negative correlations with pain were identified (Figures 5D, E), including the precentral gyrus, the lingual gyrus, the thalamus, the superior frontal cortex, and the superior temporal cortex. See Tables 2, 3. Some regions may have more than one cluster of correlated voxels. In addition, given that the ROIs are relatively large, positive correlation clusters and negative correlation clusters may appear in the same ROI (e.g., superior frontal cortex). Combining Tables 2, 3, there are 21 distinct ROIs that showed strong correlations with spontaneous pain fluctuations.

3.2. Deep learning-based analysis with CNNs

A CNN was trained and applied to fMRI data to predict pain levels. The batch size was two and 50 training epochs were used with early stopping (see Figure 3). Dividing the patients into eight groups of 4 or 5 each and training the model on seven groups and testing on the remaining group, the prediction accuracy averaged across all eight folds was 73% (Figure 6A), which is significantly higher than the chance level of 50% (P < 2 × 10^–3). The confusion matrix was given in Table 4. During the training phase, we recorded the value of the validation accuracy for each epoch, which, as shown in Figure 6B, stabilized after 15 training epochs. According to the ROC curve in Figure 6C, the area under the ROC curve [area under the curve (AUC)] for the low pain class was 70%, and the high pain class was 70%, demonstrating that the proposed CNN decoding approach could distinguish low pain and high pain levels via fMRI BOLD signals and detect more true positive and true negative samples than false positive and false negative samples.

To identify the contribution of each brain region to the CNN classification of pain levels, we performed a sensitivity analysis of the trained CNN models by occluding the BOLD activities from each brain region and examined the decoding accuracy change without the contribution of the voxels in the region (Zeiler and Fergus, 2014). A larger decrease in decoding accuracy is taken to indicate that the brain region being occluded plays a more important role in pain generation and perception. The important brain regions thus identified, as shown in Figure 6D, included lingual gyrus, superior frontal cortex, thalamus, and dorsal anterior cingulate cortex (dACC). A list of CNN-identified brain regions underlying TN pain is shown in Table 6.

3.3. Deep learning-based analysis with GCNNs

A GCNN was trained to take BOLD signals as well as fMRI FC patterns as input to predict pain levels. The batch size was one and there were 50 training epochs with early stopping. Dividing the patients into eight groups of 4 to 5 patients each, the prediction accuracy averaged across all eight folds was 72% (Figure 7A), which is significantly higher than the chance level of 50% (P < 2 × 10^–3). The confusion matrix of the proposed GCNN was given in Table 5. During the training phase, we recorded the value of the validation accuracy for each epoch, and as shown in Figure 7B, after 20 epochs, the validation accuracy stabilized. According to the ROC curve in Figure 7C, the AUC for the low pain class was 54%, and the high pain class was 66%, demonstrating that the GCNN approach worked well to decode low and high pain from function connection patterns.

We applied the same sensitivity analysis to the proposed GCNN model by occluding the BOLD activities and associated connections for each ROI. ROIs were ranked according to the prediction accuracy decrease from occlusion. Top-ranked brain regions, as shown in Figure 7D, included dACC, fusiform, and superior temporal. A list of GCNN-identified brain regions underlying TN pain brain is shown in Table 6.

3.4. Validation analyses

Both CNN and GCNN analyses are AI-inspired deep learning methods. This study, to the best of our knowledge, is the first applying these methods to neuroimaging data from TN patients. We performed two additional analyses to further test the validity of the two methods. First, it is reasonable to expect that CNN- and GCNN-predicted pain levels during the pain tracking session and the patient’s self-reported pain levels be related. Assigning low pain the value of zero and high pain the value of one, the pain ratings predicted by CNN and GCNN based on fMRI averaged over the pain tracking session were plotted against patients’ reported pain ratings average over the pain tracking session in Figures 8A, C, respectively; a significantly positive correlation was seen for both CNN and GCNN (R = 0.62 and 0.60 with both P < 2 × 10^–3), suggesting that the CNN- and GCNN-predicted pain levels from fMRI brain scans indeed reflected the level of pain experienced by the patients. Second, as indicated in Methods, in our experimental paradigm, pain tracking and resting state were recorded back-to-back. It is reasonable to expect that the average pain levels during the two sessions are correlated. If the model-predicted pain levels indeed reflect the actual pain levels, a strong correlation between model-predicted pain during resting state and that during pain tracking should exist. This was found to be the case in Figures 8B, D.

3.5. Pain tracking vs. motion tracking

As described in Methods, for the first 10 patients, 10-min of pain tracking was followed by 10 min of motion tracking. We tested, on the data from these 10 patients, whether the fMRI activities during pain tracking mainly reflected spontaneous pain fluctuations rather than movement executions. The BOLD time courses were extracted from each voxel in the pain tracking session (0–10 min) and the motion tracking session (10–20 min). We used the trained CNN models to compute the accuracy of pain rating predictions for the two sessions. The prediction accuracy for the pain tracking session was 68%, well above chance level of 50% (P = 0.04), and that for the motion tracking session was 42%, which is not different from chance (P = 0.33). This result was expected because the “pain level” indicated by the patient during motion tracking was not the actual pain level experienced by the patient at the time of motion tracking but the cursor position on the computer monitor.

3.6. Subgroup analysis

Besides the male and female subgroups, among the 39 TN patients, there were 20 who had prior surgery for their conditions, but the pain returned at the time of the experiment. The remaining 19 patients never had surgery. A question of interest is whether these subgroups of patients have common neural substrate of pain generation and perception.

To test whether the subgroups shared common neural substrate of pain generation and perception, we trained the CNN and GCNN models on one subgroup and evaluated the performance of the model on the other subgroup. First, for models trained on female patients and tested on male patients, the decoding accuracy achieved by CNN and GCNN were 76.4 and 66.6% respectively, both significantly above chance level of 50% (both P < 2 × 10^–3). We did not train models on male patients to test them on female patients because the number of male patients is too small for model training purposes (n = 14 for male patients compared to n = 25 for female patients). For models trained on non-surgery patients and tested on surgery patients, the decoding accuracy for CNN and GCNN were 55.7 and 53.3%, respectively, both significantly above chance level of 50% (both P < 2 × 10^–3), whereas for models trained on surgery patients and tested on non-surgery patients, he decoding accuracy achieved were 66 and 72% for CNN and GCNN models respectively, also both significantly above chance level of 50% (both P < 2 × 10^–3). See Figure 9. These results demonstrated that common neural substrate is shared among TN patients of different subgroups.

The clinical and psychological characteristics of the non-surgery vs. surgery and female vs. male subgroups were also evaluated. There were no significant differences between the surgery/non-surgery groups in age (P = 0.3848), disease duration (P = 0.3593), current pain intensity (P = 0.1027), or usual pain intensity (P = 0.9055) (Supplementary Figures 1, 2). When evaluating the different psychological outcomes from the different test inventories, there was a significant difference in the BDI-2 score, with surgical subjects having a significantly higher score (mean ± SD: 15 ± 9 vs 9 ± 8, P = 0.0487). The other psychological scores from the PCS (P = 0.5237), PASS (P = 0.3346) (including subscales: fear (P = 0.8547), cognitive (P = 0.4162), escape/avoidance (P = 0.4071), and physiological (P = 0.2511) were not significantly different between the surgery and non-surgery groups (Supplementary Tables 1, 2). When evaluating sex differences for these same clinical and psychological outcome measures, there was a significant sex difference in the age of the subjects, with males (mean ± SD: 65 ± 10 y.o.) being significantly older (P = 0.0025) than female subjects (mean ± SD: 53 ± 11). There were no significant differences between female and male subjects when comparing disease duration (P = 0.8997), current pain intensity (P = 0.6229), and usual pain intensity (P = 0.0759) distributions (Supplementary Figures 3, 4). The psychological scores from the PCS (P = 0.5592), BDI-2 (P = 0.5171), and PASS (P = 0.5777) (including subscales: fear (P = 0.7299), cognitive (P = 0.9528), escape/avoidance (P = 0.5393), and physiological (P = 0.0776) tests were also not significantly different between female and male subjects (Supplementary Tables 3, 4).

3.7. Signature centers of TN pain

Three analyses were applied in this study to identify signature centers of TN pain: correlation analysis, CNN analysis, and GCNN analysis. We listed the top 23 brain regions that appeared in at least two of the analyses in Table 6. Six regions were identified by all three methods, including superior temporal cortex, insula, fusiform gyrus, precentral gyrus, superior frontal gyrus, and supramarginal gyrus. The 17 remaining regions included dACC, thalamus, lateral occipital gyrus, inferior temporal cortex, postcentral gyrus, lingual gyrus, and inferior parietal gyrus. All these regions have been implicated in chronic pain especially TN pain in past studies (see section “4. Discussion”).

4. Discussion

In this study we sought to identify the neural substrate of TN pain by recording and analyzing fMRI data from TN patients while they rested or tracked their spontaneous fluctuations in pain levels. By applying both conventional and AI-inspired approaches, we obtained converging evidence implicating a common set of brain regions, including the insula and the precentral gyrus, as playing an important role in the generation and perception of TN pain. Additional complementary insights into TN pain-related brain regions are offered by each of the three approaches.

4.1. Signature pain centers identified by three analyses

Table 6 shows that there are six brain regions that appeared in all three analyses, including superior temporal cortex, insula, fusiform gyrus, precentral gyrus, superior frontal gyrus, and supramarginal gyrus. All these regions have been implicated in TN pain in prior studies (Tsai et al., 2018). Yuan et al. (2018) and Xiang et al. (2019) reported increased ReHo in the superior frontal gyrus in TN patients compared to healthy controls; cortical thickness was shown to decrease with increase in pain duration (Moon et al., 2018). The superior temporal gyrus (STG), functionally known for its role in memory and language processing (Burton, 2005), showed gray matter volume abnormalities in TN, and the degree of grey matter volume alterations in the left STG may reflect pain severity. The abnormal structural alterations in the temporal lobe may be due to the generation and/or maintenance of emotions and perception of pathological pain (Tang et al., 2020). The important role of precentral and postcentral gyrus in pain processing is well established. The primary motor cortex, located in the precentral gyrus, may be involved in inhibiting motion so as not to aggravate pain conditions (Apkarian et al., 2005), whereas the postcentral gyrus, where the primary somatosensory cortex is located, is the cortical gateway through which nociceptive input is processed and transmitted to other brain structures (Peyron et al., 1999; Bantick et al., 2002; Valet et al., 2004; Hu and Iannetti, 2019). In a resting-state fMRI study, Wang et al. (2015) found that ReHo in the precentral gyrus and the postcentral gyrus is positively correlated with patients’ daily experienced pain severity, suggesting a link between local synchronization of intrinsic brain activity and pain modulation. In addition, in a recent meta-analysis, the somatosensory cortex is found to be structurally changed in TN compared with healthy controls (Henssen et al., 2019). The fusiform gyrus, a visual area known for its involvement in face-related processing, is thought to be important in mediating mental imagery processes related to pain perception (Diekhof et al., 2011; Ter Minassian et al., 2013). In TN, the fusiform is shown to have decreased gray matter volume compared to controls (Parise et al., 2014; Li et al., 2017). In addition, the left fusiform gyrus is found to be involved in pain anticipation and perception, leading to a negative correlation with pain ratings (Ter Minassian et al., 2013). The insular cortex has also been implicated in mediating pain intensity as well as negative emotions (Desouza et al., 2013). A recent MRI study found altered insular morphology and FC and abnormal diffusion parameter in the white matter adjacent to the insular cortex in TN (Wang et al., 2018). Furthermore, significant local gyrification index (LGI) reductions in the left insular cortex were found in patients with TN compared with control groups (Wang et al., 2018). The authors argued that pain perception results from nociceptive representation being transformed into subjective magnitude assessment within the insula (Baliki et al., 2009). Tian et al. (2016) found that TN patients exhibited significantly increased long-range FC density in the right supramarginal gyrus. As can be seen, these past studies rely either on structural information or on resting state data. The functional roles of these regions in TN pain generation and perception remain to be better established. Our results, by utilizing a functional paradigm in which TN patients tracked pain levels, shed new light on this issue. In particular, we found that the neural activities in these six regions not only closed tracked the pain level fluctuations (correlation analysis and CNN analysis), they also mediated network-level communications among different brain regions (GCNN analysis) during TN pain.

4.2. Signature pain centers identified by two of the three analyses

Among the common set of regions identified by both AI-inspired deep learning analyses but not the correlation analysis, dACC is known for its role in pain processing (Henssen et al., 2019); it is consistently activated in human imaging studies of pain. As a central hub in the pain matrix, the dACC is highly connected to other brain areas involved in cognition, emotion, and negative affect, all of which are associated with chronic pain (Zeidan et al., 2011; Eisenberger, 2012; Roy et al., 2012). It has been further suggested that the cingulate cortex is also important for the transition from acute to chronic pain (Tsai et al., 2019). Previous TN works, including the Yuan et al.’s (2018) resting-state fMRI data and Moisset et al.’s (2011) task fMRI data, have found abnormal activities in ACC. Structurally, Mo et al. (2021) found that TN patients exhibited reductions in cortical indices in the ACC, the midcingulate cortex (MCC), and the posterior cingulate cortex (PCC) relative to healthy controls group, indicating that the ACC may play a role in pain adaptation, habituation, distraction, and the engagement of the endogenous pain control system (Qu et al., 2011). In addition to dACC, the thalamus, the postcentral, the lingual, and the inferior temporal gyrus are found to be activated in two of the three analyses. These regions have also been implicated in previous research on TN pain. In particular, the previous research has firmly established the importance of the thalamus in pain processing; it receives nociceptive sensory information from the periphery, integrates this information with arousal and attention, and sends outputs to broad regions of the cerebral cortex for further processing (Coghill et al., 1999; Peyron et al., 2000; Tracey and Mantyh, 2007). In TN, the grey matter increases in the thalamus for TN patients relative to controls, suggesting a link between thalamic structural change and activity-dependent plasticity in S1 via thalamocortical projections (Woolf and Salter, 2000; Desouza et al., 2013). In our data, whereas the GCNN analysis was not able to show that the thalamus is an important region underlying TN pain, the correlation and the CNN methods both found that the thalamus is important for predicting the spontaneous fluctuations of TN pain, consistent with the established role of the thalamus in pain processing including TN pain (Gustin et al., 2011). Regarding the inferior temporal gyrus, past work found that the gray matter volume of the left inferior temporal gyrus was negatively correlated with current pain intensity and disease duration in TN patients (Parise et al., 2014; Wang et al., 2017a). Decreases in the amplitude of low-frequency fluctuation (ALFF) in the right inferior temporal region were found in TN patients (Wang et al., 2017b). Zhu et al. (2020) also found that compared with the healthy control (HC) group, patients with TN showed the degree centrality value (calculated by counting the number of significant suprathreshold correlations (the degree of the binarized adjacency matrix) for each individual) changed in the right lingual gyrus, right postcentral gyrus, left paracentral lobule, left inferior cerebellum, and right inferior cerebellum.

To further contrast the conventional correlation analysis against the AI-inspired methods, it is informative to combine the regions identified by the two AI-inspired approaches and compared them with the ones identified by the correlation method. Although both types of methods have identified important pain-related regions, such as the thalamus, the correlation method failed to identify other significant pain-related regions, including the caudate, dACC, and postcentral gyrus, all areas of the so-called pain matrix; these regions were readily identified by the AI-inspired approaches (Ab Aziz and Ahmad, 2006; Freund et al., 2009; Iannetti and Mouraux, 2010; Lieberman and Eisenberger, 2015; Zhu et al., 2020). The power of the AI-inspired approaches may derive from their ability to model non-linear structures in the data. On the other hand, the correlation approach has identified regions that do not overlap with that identified by the AI-inspired approaches; e.g., the frontal pole. We thus suggest that integrating the conventional and AI-inspired approaches is a fruitful direction going forward for gaining comprehensive insights into the neural substrate of TN pain generation and perception.

4.3. Clinical considerations

The clinical TN population that was recruited for this study was consistent with the known disease demographics and characteristics: female > male and age > 50 years old (Jainkittivong et al., 2012). The majority of subjects (n = 16) reported moderate to severe usual pain levels, seven of which reported maximal VAS pain levels of 100 (Supplementary Figure 2B). At the time of the scanning procedure, the current pain intensity was distributed roughly equally across zero to maximum pain range, indicating the spontaneous nature of this pain disorder. There was a significant age difference between the males and females; however, with the males being on average older; the disease duration was similar between males and females. This likely indicates that males were diagnosed later in life as compared to females. The only psychological factor that was elevated was the BDI-2 score for surgical subjects, as compared to non-surgical subjects. While significant (P = 0.0487), the average male BDI-2 score is considered within the “mild-to-moderate depression” range (mean ± SD: 15 ± 9), while the female scores were within the no or minimal depression range (mean ± SD: 9 ± 8). On average, the PCS ratings for the subjects are within the mild range (21–40), and the PASS scores were not significantly different when comparing surgical status or sex differences. In the context of these clinical findings, our results that both CNN and GCNN models constructed based on patients from one subgroup (surgical vs. non-surgical, male vs. female) can decode the patients from the other subgroup can be seen as reflecting shared neural substrate rather than driven by differences in clinical conditions.

4.4. Methodological considerations

The correlation analysis has been used in numerous pain studies and remains the principal method of percept-related fMRI for identifying the neural substrate of naturally occurring pain (Apkarian et al., 2001). It is simple and intuitive (Desouza et al., 2013). The weaknesses are that it only detects linear relationships, and the statistical effects are not strong. In our data, if we applied any kind of multiple comparison correction, we would find no activations. Deep learning is an emerging area of machine learning, and it is in the early stages of being applied to analyze neuroimaging data. In this work, in addition to showing that CNN and GCNN are able to predict pain fluctuations from fMRI data, we did several additional analyses to validate the approach. First, we showed that the predicted pain levels by CNN and GCNN are correlated with the pain ratings indicated by the patients, as would be expected. Second, the CNN and GCNN predicted pain levels from the resting-state data and from pain tracking data are highly correlated, again as would be expected, demonstrating that the CNN and GCNN model predictions are robust. Combining these AI-inspired methods with conventional method to seek converging evidence may become a promising way in future studies of pain neuroimaging.

4.5. Limitations

This study has a number of limitations. First, for the correlation analysis, as mentioned above, if a whole-brain multiple comparison approach such as false discovery rate (FDR) was applied, no statistically significant regions will appear in the activation map, despite having a reasonable sample size of n = 39. The situation is similar in prior pain studies utilizing the percept-related fMRI analysis (Geha et al., 2007). Nevertheless, the thalamus, an established region in pain processing including TN pain, only appeared in the CNN and correlation analysis map (Lenz et al., 2004; Danyluk et al., 2021) but not the GCNN map. This demonstrates that, despite the statistical weakness associated with the correlation analysis, it can still provide important information which not all AI-inspired methods can provide. Second, because of the need to increase the pain tracking session length without increasing the overall scanning time, we were only able to record motion tracking from the first 10 patients. Nevertheless, by applying the AI-inspired models to the data from pain tracking and from motion tracking in these 10 patients, we were able to establish that the fMRI data recorded during pain tracking are mainly driven by spontaneous pain fluctuations, not by the movement associated with pain tracking. Third, when assigning discrete class labels to the data, the threshold of 15 was chosen to achieve the class balance between low and high pain. This high vs. low pain level demarcation was a trade-off between clinical considerations and technical requirements and is not a clinically applicable pain level designation.

5. Conclusion

We applied advanced statistical methods to patterns of brain activation related to the paroxysms of TN pain and generated a set of “signature centers” of pain generation and perception within the brain. Our approach, combining both conventional and AI-inspired methods to yield converging findings, is highly novel in the context of TN research, and our result provides a sorely needed basis for understanding the central mechanisms of TN. We hope that the insights revealed in this study can lead to a better understanding of TN and curing of this devastating condition.

Data availability statement

The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors.

Ethics statement

The studies involving human participants were reviewed and approved by the WCG Institutional Review Board (IRB). The patients/participants provided their written informed consent to participate in this study.

Author contributions

YL involved in study design, data analysis, and manuscript writing. QZ involved in study design, subject recruitment, data collection, and data analysis. ZH, KB, and SM involved in subject recruitment and data collection. JN and MD involved in study design, subject recruitment, data collection, data analysis, and manuscript writing. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the Facial Pain Research Foundation (FPRF).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnhum.2023.1144159/full#supplementary-material

Abbreviations

AI, artificial intelligence; ALFF, amplitude of low-frequency fluctuation; ASA, American Society of Anesthesiologists; ASD, autism spectrum disorder; AUC, area under the curve; BDI-2, beck depression inventory-II; BOLD, blood-oxygen-level-dependent; CBP, chronic back pain; CNN, convolutional neural networks; dACC, dorsal anterior cingulate cortex; DNN, deep neural network; EPI, echo-planar imaging; FC, functional connectivity; FDR, false discovery rate; fMRI, functional magnetic resonance imaging; FOV, field of view; FPRF, Facial Pain Research Foundation; GCNN, graph convolutional neural network; GIN, graph isomorphism network; GLM, general linear model; HC, healthy control; HRF, hemodynamic response function; IBS, irritable bowel syndrome; ICHD-3, the 3rd edition of the International Classification of Headache Disorders; IHS, International Headache Society; IRB, Institutional Review Board; LGI, local gyrification index; LR, logistic regression; MCC, midcingulate cortex; MVD, microvascular decompression; OHSU, Oregon Health Science University; PASS, pain anxiety screening scale; PCC, posterior cingulate cortex; PCS, pain catastrophizing scale; ReHo, regional homogeneity; ROC, receiver operating characteristics; ROI, region of interest; rs-fMRI, resting-state fMRI; SD, standard deviation; SGD, stochastic gradient descent; SPM, statistical parametric mapping; SVM, support vector machine; TE, time of echo; TN, trigeminal neuralgia; TR, repetition time; VAS, visual analog scale.

References