Topology of pain networks in patients with temporomandibular disorder and pain-free controls with and without concurrent experimental pain: A pilot study

Introduction

Chronic temporomandibular disorders (TMD) involve persistent pain in the masticatory muscles and jaw joints and in many cases dysfunction of the jaw. Peripheral nociceptive pain elicited by inflammation or injury in the temporomandibular joint and surrounding muscles and tendons can contribute to clinical pain in the area, but central nervous system changes in at least a subset of individuals with TMD can also lead to centralized, or nociplastic, pain (1). Decades of research in a variety of regional pain conditions including TMD, chronic pelvic pain, irritable bowel syndrome, migraine, vulvodynia, and others, as well as more widespread pain as is observed in fibromyalgia, have shown that chronic pain is associated with differences in brain structure, function, and neurochemistry, compared to chronic pain-free individuals (2). In addition, many patients with chronic pain exhibit alterations in experimentally evoked noxious stimulus processing, including hyperalgesia, allodynia, and impaired ability to reduce ascending nociceptive signals through endogenous activation of the descending pain modulatory system (3, 4). Thus, brain neuroimaging of TMD patients who are undergoing experimental pain could reveal underlying functional differences in nociceptive processing.

Central nociceptive functions are subserved by cortical-subcortical pathways, including the trigemino-thalamocortical (CN V-TC) and periaqueductal-rostral ventromedial medulla (PAG-RVM) pathways, as well as distributed cerebral networks, such as the sensory/discriminatory lateral and affective/cognitive medial pain systems (LPS and MPS, respectively) (5–8). These nociceptive systems, however, do not function in isolation. They are closely associated with other cortical and subcortical networks which primarily subserve affective, cognitive, motivational, and integrative functions, but are jointly activated with the above nociceptive systems and pathways during pain sensation. For example, a corpus of brain areas, originally described as a “pain neuromatrix,” are involved in the nociceptive (posterior insular, middle cingulate, and medial parietal opercular cortices) or perceptual aspects of pain (middle and anterior insular, prefrontal, posterior parietal, supplementary motor, and anterior cingulate cortices, hippocampus, striatum, and cerebellum) (5, 9–16); however, it should be noted that most of these areas are not specific to pain processing. Additionally, alterations in grey matter volume, neurochemistry, and connectivity within the default-mode network (DMN), which functions in monitoring and processing internal states (17–22), are correlated with pain sensitivity, severity, duration, catastrophizing, and rumination (16, 23–26). Finally, several regions which mediate body perception, learning and motivational functions, and pain anticipation, perception, and empathy, are often omitted from even the broadly-defined pain neuromatrix. This is due, in part, to the absence of a consensus definition of what areas constitute this neuromatrix and a lack of specificity of these regions for processing nociception. These regions include the fusiform cortex, a component of the dopaminergic reward system which is not only implicated in body and face recognition but also in pain anticipation and perception (27–32); other components of the dopaminergic reward system, such as the ventral pallidum and planum polare (33–35); the supramarginal and angular gyri, which are engaged in proprioception, body perception, and pain perception and empathy (36–38); the putamen, which subserves motor, learning and motivational functions (39–41); and frontopolar cortex, which inhibits the DMN and mediates among other intrinsic connectivity networks (42). Some areas have been implicated in regional pain syndromes, such as maxillofacial pain (e.g., angular and transverse temporal gyri) and may respond in a graded fashion as a function of pain intensity (e.g., frontopolar cortex) (43–45).

Previous attempts to characterize the connectivity of regions subserving pain perception and its sequelae have generally been limited to the well-defined CN V-TC and PAG-RVM pathways, LPS and MPS, DMN, and several of the miscellaneous or ancillary regions (5, 16, 26, 46–51). Adhering to the definition of the pain neuromatrix suggested by Melzack (9, 52), however, which includes all systems jointly activated by pain sensation and associated with perceptual-attentional, cognitive, affective, motivational, and nociceptive functions, little is known about its global topology. Specifically, it is unclear which connections have primacy (serve as network hubs) or how this topology changes in the presence of chronic pain syndromes. This exploratory study leverages resting-state functional MRI (rsfMRI) to begin to address these issues in a cohort of pain-free controls and individuals with painful TMD under basal and evoked-pain conditions, resulting in a preliminary description of the nociceptive network topology and its alteration in chronic pain conditions. TMD is one of several chronic overlapping pain conditions believed to be driven, at least in part, by central nervous system alterations. Fibromyalgia is considered by most to be a prototypical nociplastic pain condition, and higher “fibromyalgianess” scores have been found to be associated with more complex pain mechanisms and poorer outcomes in TMD and other chronic pain conditions (53–56). Therefore, in addition to the comparisons between pain and pain-free groups, we also assessed whether TMD patients degree of likely nociplastic pain (assessed with the fibromyalgia survey criteria) is correlated with the degree of central nervous system changes. The present study builds upon prior work utilizing rsfMRI and task-based fMRI which demonstrated changes in the CN V-TC, LPS and MPS, motor system, and DMN in individuals with TMD (see, for example, the recent review by 5).

Methods

Study population and testing

This study was approved by the University of Michigan Institutional Review Board and all subjects provided informed consent. 16 patients diagnosed with painful TMD (15 female, 1 male, aged 18–49 years) and 12 pain-free controls (11 female, 1 male, aged 19–57 years) participated. Participants with TMD were recruited from the Oral and Maxillofacial Surgery Center and Hospital Dentistry at the University of Michigan and controls were recruited through the University's online research match system. Two experienced physicians determined if participants with TMD met the diagnostic criteria and patients were subsequently included if they had nonspecific pain in the jaw/cheek area for a least 6 months. Pain-free control participants were excluded if they had a history of chronic pain lasting more than 6 months. All participants were excluded if they had physical impairments (e.g., amputation, blindness or deafness), severe psychiatric illnesses (e.g., major depression with suicidal ideation, schizophrenia, or substance abuse within the last 2 years), pregnancy, current tobacco use (>5 cigarettes/day, equivalent tobacco containing product), history of alcoholism, use of opioids including tramadol, sedatives including benzodiazepines or hypnotics, NMDA receptor antagonists (e.g., memantine), and eugeroics (e.g., modafinal) in the absence of appropriate washout periods, and the unstable use of the following medications: monoamine oxidase inhibitors (MAOIs), tricyclic antidepressants, selective serotonin reuptake inhibitors (SSRIs), serotonin norepinephrine reuptake inhibitors (SNRIs), antiepileptics, muscle relaxants. TMD and control groups did not differ with respect to age ($p=0.399$, Mann–Whitney U test). Participants were administered the Pennebaker Inventory of Limbic Languidness (PILL) (57), which yielded a hypervigilance score, as well as the 2010/2011 American College of Rheumatology Fibromyalgia Survey (ACRFS) (58, 59), which yielded widespread pain index (WPI), symptom severity (SS), and total ACRFS scores. Hypervigilance, WPI, and SS scores are subjective measures of patients' attention to aversive sensations, distribution of pain among 19 body areas, and general and domain-specific effects of pain symptoms, respectively (57, 59–61). Scores, as well as demographic and pain rating data, are summarized in Table 1.

MRI acquisition and preprocessing

Two fMRI datasets were acquired during the MRI session. Before the start of the scan, a 13.5 cm-wide velcro-adjusted pressure cuff, connected to a rapid cuff inflator (Hokanson E20 AG101, Hokanson Inc, Bellevue, WA, USA), was attached to the leg over the left gastrocnemius muscle, with the top of the cuff ∼4 cm distal to the bottom of the kneecap. The first rsfMRI scan was conducted with the cuff deflated (i.e., 0 mmHg). The second resting state scan was conducted with the cuff inflated to 200 mmHg; except, in cases where that pressure was intolerable, a lower pressure was used. Tolerability was determined by a 1-min inflation of the cuff outside of the scanner and again more briefly (5-s inflation) just before the start of the scan. If the pressure was intolerable, the pressure was adjusted in 20 mmHg increments until the participant found it tolerable. The cuff was inflated simultaneously with the start of the scan and deflated at the end. Perceived pain of the cuff inflation was obtained by having participants rate painfulness of the cuff on a scale from 0 to 100, where 0 means “no pain” and 100 means “the most intense pain imaginable.” These ratings were collected immediately following the cuff inflated and deflated runs.

Cuff-deflated and cuff-inflated (induced-pain) rsfMRI data were acquired as separate but consecutive scans (conditions) within a single imaging session on a 3.0 T GE Discovery MR750 (GE Medical Systems, Waukesha, Wisconsin) with a 32-channel Nova head coil. Echo-planar (EPI) rsfMRI data was acquired over a continuous 5-min period (resolution 2.386 mm × 2.386 mm × 2.500 mm, interleaved slices, TR = 1,200 ms, TE = 30 ms, FA = 70°, FOV = 210 mm, NEX = 1, multiband factor = 3). Two gradient-echo field maps (6.5 and 8.5 ms TE) and one T1-weighted 3D-spoiled gradient recalled echo (SPGR) structural scan (650 ms TR, 14.653 ms TE, 90° flip angle, echo train length (ETL) = 28) were also acquired for each subject. rsfMRI data was subjected to the default minimal preprocessing and denoising pipelines in the CONN Toolbox v19c (62, 63). These pipelines are described in detail in the Supplementary Material. All structural and denoised functional data, grey matter, white matter, and CSF masks, as well as regions of interest (see below), were manually inspected for quality assurance and to confirm registration validity. Selected quality assurance results are also presented in the Supplementary Material.

Surface-based morphometry (SBM) and SBM regions of interest

Volumetric segmentation of the T1-weighted MPRAGE was performed with Freesurfer (v6.0: surfer.nmr.mgh.harvard.edu). The Freesurfer segmentation processing stream, which is described in detail elsewhere (64–67), consists of segmentation of grey/white matter and grey/cerebrospinal fluid tissue classes as well as the subcortical white matter and deep gray matter volumetric structures, and is not restricted to the voxel resolution of the original data, but is capable of detecting submillimeter differences between groups (64, 67–70). Freesurfer's “cross-sectional” workflow consists of automated brain extraction, segmentation, intensity normalization, and topological correction using intensity and surface continuity information from two channels represented by participants' T1 and surface continuity information from participants' T1-weighted images (71). Statistical evaluations of grey and white matter differences between the control and TMD groups were performed at the cluster level by multivariate general linear model (GLM) using Freesurfer's Qdec utility as well as via ROI-based GLM analyses in SPSS. A smoothing kernel of 10 mm full width at half maximum (FWHM) was applied to all segmentations prior to statistical analysis.

Cluster-based analyses

For cluster-based analyses, two sets of GLMs were conducted, the first as a group-wise analysis using a binary encoding of “group” as a fixed factor and “age” and total estimated intracranial volume (“eTIV”) as covariates, and the second as a regression-type analysis using ACRFS total scores, WPI, symptom severity, or Pennebaker scores as predictor variables and “age” and “eTIV” as covariates. In both cases, correction for multiple comparisons was performed by Monte Carlo simulations at $Z\ge 1.30\sim \mathrm{family}-\mathrm{wise}\mathrm{error}(\mathrm{PFWE})P_{FWE}\le 0.05$. Six suprathreshold clusters, including the left pre- and postcentral gyri (L PreCG, L PostCG), left supramarginal gyrus (L SMG), and left superior temporal sulcus and temporo-occipital and posterior middle temporal gyri (L toMTG, L pMiTG, L bSTS) were indicated as “of interest” by these cluster-based SBM analyses (see “Results” section); of these, L SMG was already included as an atlas-defined ROI. L PreCG, L PostCG, L toMTG, and L pMiTG were back-projected to subject space as volumetric regions of interest (ROIs) and leveraged in connectivity analyses, as described below, and included in Supplementary Table S1.

A priori regions of interest

Additional ROIs for the present study were selected for their association with the trigemino-thalamocortical, or CN V-TC, pathway (5, 6), periaqueductal grey-rostral ventromedial medulla, or PAG-RVM, pathway (5, 72–75), lateral and/or medial, or LPS and MPS, pain systems (5, 7, 8, 76), the default-mode network, or DMN (17–22), or other components of the “pain neuromatrix” (5, 7, 9–11, 13–16, 77). These ROIs are listed in Supplementary Table S1. Although most ROIs were predefined in the CONN Toolbox as per the AAL and FSL Harvard-Oxford Maximum Likelihood Cortical and Subcortical atlases (78–81), several required definition from MNI coordinates reported in prior literature or were more precisely delineated from direct segmentation of participants' T1 volumes using Freesurfer 6.0 or SUIT v3.4 (82–86). Coordinate-defined ROIs included the periaqueductal grey (PAG), which was comprised of a medial, 6 mm-radius spherical ROI centered at MNI coordinate [$0,-32,-12$] per Coulombe, Lawrence et al. (87), and lateral and medial hypothalamus ROIs, comprised of 2 mm-radius spherical ROIs centered at [±6, −9, −10] and [±4, −2, −12], respectively, per Li et al. (88) Rostral and caudal anterior cingulate (rACC, cACC), brainstem structures [medulla, pons, midbrain, and superior cerebellar peduncle (SCP)], hippocampus (HC), amygdala (AMG), and entorhinal (EnRC) and anterior and posterior parahippocampal cortices (aPaHC, pPaHC) were delineated from T1 segmentations in Freesurfer and SUIT. In all cases, registration validity of stereotactically-defined ROIs was confirmed by visual inspection against individual structural data following back-projection to subject space. The total corpus from all sources was comprised of 86 ROIs.

ROI-ROI connectivity matrix generation and inference testing

Four connectivity matrices were generated (control-deflated, control-inflated, TMD-deflated, and TMD-inflated) by computing the matrix $\mathbf{Z}$ of Fisher-transformed, bivariate correlation coefficients among the 86 ROIs (see Supplementary Table S1) using the CONN Toolbox. These coefficients were computed from the blood oxygenation level dependent (BOLD) BOLD timeseries within each pair of ROIs $R_i$ and $R_j$:

where

Assessments of functional connectivity as a function of group, condition, and Pennebaker and ACRFS measures were performed via GLM. For first-level GLM analyses, design matrices comprising group (control and TMD), condition (cuff-deflated and cuff-inflated), and hypervigilance, WPI, and symptom severity were defined while controlling for physiological and motion artefacts using the motion parameters generated during preprocessing. Second-level, within-subjects analyses on condition were conducted by leveraging individual ROI-ROI connectivity matrices within either the control (control inflated versus deflated) or TMD group (TMD inflated versus deflated) in paired-samples t-tests. Second-level, between-subjects analyses on group leveraged individual contrast images within either the deflated (TMD versus control while deflated) or inflated (TMD versus control while inflated) in independent-samples t-tests. Finally, second-level analyses against hypervigilance, WPI, and symptom severity scores were conducted as regression-type GLMs across all subjects and conditions. One analysis was performed for each of these scores by specifying hypervigilance, WPI, or symptom severity score as the predictor variable and ROI-ROI connectivity as the dependent variable. All regression models included an intercept term.

Inference testing for all GLM analyses was performed using network-based statistics (NBS), which treats ROI-ROI connectivity matrices as network graphs and subgraphs and extend the principles of cluster-based thresholding (89–91) to the graph model. The NBS approach controls for the familywise error rate for each subgraph, thereby increasing the contrast-to-noise ratio (92). These attributes are particularly attractive for the present study, wherein the correlation between an effect (group, condition, or Pennebaker or ACRFS scores) and topology of specific networks, rather than pairwise connection strengths between ROIs, is of interest, and where the contrast-to-noise ratio is of critical importance due to small sample size. Computation of the NBS proceeds by first applying an a priori “height” threshold ($p<0.001$, uncorrected) to the ROI-ROI matrix of t- or F-statistics to yield a suprathreshold network among all ROIs. Second, this network is decomposed, by breadth-first search, into component subgraphs which are each characterized by a “network mass,” defined as the sum of $t^2$- or F-statistics over all connections within the subgraph, and the mass of each subgraph compared with a null distribution estimated from multiple permutations of the initial data. Third, subgraph-level FWE-corrected p-values (the likelihood of finding a subgraph with the same or larger mass across all ROI-ROI connections) and false discovery rate (FDR) FDR-corrected p-values (the expected proportion of false discoveries among networks with the same or larger mass across all ROI-ROI connections) are computed.

Adjacency matrix generation and graph-theoretic analyses

Group × condition adjacency matrices (control deflated, control inflated, TMD deflated, and TMD inflated) were generated from one-sample permutation tests against zero and corrected for multiple comparisons using NBS in the CONN Toolbox, as described previously (92–94). Suprathreshold and subthreshold connections were binary-encoded as 1 and 0 respectively, and a “core” adjacency matrix, consisting of connections which were suprathreshold across groups and conditions, generated by elementwise multiplication of the control-deflated, control-inflated, TMD-deflated, and TMD-inflated adjacency matrices. Nodal degree and betweenness centrality were also computed from the four group × condition adjacency matrices (but reported only for “core” connections in the “Results” section) using the Brain Connectivity Toolbox (95). These graph-theoretic metrics have been described in detail previously, and the reader is referred to those publications for a full discussion of the theoretical basis of the metrics (95–99). They are also described in detail in the Supplementary Material.

To assess differences in betweenness centrality and nodal degree among well-defined systems and pathways, ROIs were categorized as CN V-TC, PAG-RVM, LPS, MPS, or DMN, and a linear mixed-effects model conducted in Python statsmodels version 0.12.2 (100). Other ROIs (“pain neuromatrix” or “not otherwise specified”) were excluded from this analysis due to the ambiguity of their role in nociceptive processing. Separate linear mixed-effects models were conducted for the dependent variables centrality and degree. Each model was specified as a variance component model (101) using group (control or TMD) as the between-subjects factor and condition (deflated or inflated) and ROI “family” (CN V-TC, PAG-RVM, LPS, MPS, or DMN) as nested within-subjects factors, i.e.,

where a colon represents nested factors and DV the dependent variable (centrality or degree), and the variance component specified as

where 0 represents a random intercept and C denotes the categorical nature of the variable. Contrasts for “ROI family” were specified as “family” versus “non-family,” e.g., “LPS ROIs” versus “non-LPS ROIs” (all ROIs not specified as having membership within the LPS). In this example, the mixed-effects model would be specified with group (control = 0; TMD = 1), condition (deflated = 0; inflated = 1), and LPS (true = 0; false = 1) as nested variables and betweenness centrality or nodal degree as continuous dependent variables.

Finally, resolution of communities among nodes in the four adjacency matrices and “core” adjacency matrix was estimated in the Brain Connectivity Toolbox using the Louvain community detection algorithm.

Results

ACRFS and Pennebaker scores

TMD participants demonstrated more widespread pain [WPI: control median ± intraquartile range (IQR), 0 ± 1, TMD 3 ± 3, $p<0.001$, Mann–Whitney U test], greater sensitivity to aversive sensations (hypervigilance: control median ± IQR, 6.00 ± 8.00, TMD 14.00 ± 13.00, $p=0.001$, Mann–Whitney U), and heightened effects of pain on fatigue, cognitive symptoms, and other domains (symptom severity: control median ± IQR, 1.00 ± 2.00, TMD 4.00 ± 5.00, $p=0.001$, Mann–Whitney U), relative to controls. ACRFS total scores also differed between control (median ± IQR, 2.00 ± 3.00) and TMD participants (7.00 ± 9.00, $p<0.001$, Mann–Whitney U). Cuff pressure was statistically equivalent in both groups (control median ± IQR, 180.00 ± 60.00 mmHg, TMD 180.00 ± 40.00 mmHg, $p=0.755$, Mann–Whitney U) as was perceived pain intensity (control median ± IQR, 40.00 ± 23.00, TMD 40.00 ± 40.00, $p=0.516$, Mann–Whitney U; see Table 1).

Surface-based morphometry analysis

The results of the SBM cluster-based analyses are shown in Figure 1. Control and TMD groups did not differ significantly with respect to total eTIV (F = 1.476, p = 0.237). A GLM on group, controlling for eTIV and corrected for FWE by Monte Carlo simulation, indicated that TMD subjects exhibited decreased left-hemisphere superior temporal surface area relative to controls. The center of mass of this cluster, at MNI [−54, −25, −6], was localized to the ventral-posterior aspect of the superior temporal sulcus. TMD subjects also exhibited reduced left superior temporal volume with a center of mass coordinate of [−64, −49, 15], which was localized to the supramarginal gyrus. Additional GLMs, conducted across all subjects as regressions, indicated that higher ACRFS WPI and ACRFS total scores were associated with decreased white matter volumes deep to the left parietal lobe, with centers of mass located at [−27, −19, 30] (corticospinal tract/superior corona radiata) for both effects, and that higher Pennebaker hypervigilance scores were associated with decreased grey matter volume across the left precentral and postcentral gyri, or the dorsal aspect of Brodmann area 6, with a center of mass located at [−39, −19, 67]. The analyses failed to resolve any effects of symptom severity on either grey or white matter volume, nor were any effects, among group, ACRFS total scores, Pennebaker scores, or age, associated with an increase or decrease in cortical thickness or area.

ROI-ROI connectivity and graph theory analysis

Group×condition connectivity matrices

Group adjacency matrices are presented in Figure 2 (deflated condition) and Figure 3 (inflated condition). Connections were considered “present” in a group/condition if they were suprathreshold following NBS, and “absent” otherwise. Across groups and conditions, adjacencies were noted within functionally and anatomically homologous ROIs. These ROIs were grouped into ten major classes: (1) brainstem and rhombencephalic structures (medulla, pons, brainstem, superior cerebellar peduncle); (2) basal ganglia; (3) cingulate structures (subcallosal cortex, posterior cingulate and rostral and caudal anterior cingulate); (4) medial temporal and associated structures (amygdala, accumbens, entorhinal and parahippocampal cortices, hippocampus, and anterior and dorsal posterior insula); (5) frontal lobe (dorsolateral prefrontal, caudal middle frontal, inferior frontal, orbitofrontal, and frontopolar cortices); (6) fusiform areas; (7) occipital lobe; (8) pre-, postcentral, and supplementary motor areas; (9) supramarginal and angular gyri, planum polare, and parietal operculum; and (10) temporal areas (transverse temporal and right middle temporal gyri and the bank of the superior temporal sulcus). Additionally, thalamic, hypothalamic, and cerebellar ROIs comprised minor classes.

Connections between the basal ganglia and accumbens nuclei, basal ganglia and posterior temporal fusiform cortex, and posterior temporal fusiform cortex and superior and inferior lateral occipital cortices were present in the deflated condition in both groups. TMD participants exhibited more extensive connections overall, but particularly with respect to thalamic, medial temporal, rhombencephalic, frontal, temporal, and occipital ROIs, and more amygdalar connections with caudate, pallidum, entorhinal, and right anterior insula were present than in the control participants. In the inflated condition, control participants exhibited bilateral connections between pallidum, caudate, and right accumbens, left and right hippocampus and pons, inferior frontal and dorsolateral prefrontal, thalamus and frontopolar, and frontopolar and dorsolateral prefrontal ROIs, which were not present in the deflated condition, whereas connections between left anterior parahippocampal areas, the pons, and the brainstem, between the angular and supramarginal gyri, and between the right accumbens, left entorhinal, and right entorhinal ROIs, were absent. Conversely, TMD participants exhibited connections between fusiform cortex, caudate, and pallidum, left anterior insula and parahippocampal cortices, and between the amygdala, entorhinal, and left thalamus ROIs which were not present in the deflated condition, whereas connections between left anterior parahippocampal and rhombencephalic ROIs were absent, and more connections were present in TMD participants than control participants in the inflated condition. Connection strength did not differ significantly in a between-groups, within-condition analysis against either condition, or in a within-groups, between-condition analysis against the control participants. Within TMD participants, however, thalamic-right anterior insular and thalamic-right orbitofrontal connections were stronger (L thalamus-right aINS $t=+4.47$, $P_{unc}=0.004$, R thalamus-R aINS, $t=+4.14$, $P_{unc}=0.0008$), and right anterior insular-right temporo-occipital fusiform connections weaker ($t=-4.41$, $P_{unc}=0.0004$), in the inflated condition relative to the deflated condition (see Figure 4A and Table 2B; control participants were not included in Table 2B as there were no differences as a function of condition). These associations were significant at the connection level ($P_{unc}$ values listed above) and for the NBS “mass” of the pallidum-putamen-rostral anterior cingulate (graph $P_{FWE}$ = 0.0350) and orbitofrontal-thalamus-fusiform-anterior insula networks (graph $P_{FWE}$ = 0.0170), but not for their discrete subgraphs (subgraph $P_{FWE}$ = 0.1464–0.4919).

Functional connectivity versus ACRFS and Pennebaker scores

GLM analyses of ACRFS and Pennebaker scores failed to resolve statistically-significant differences in functional connectivity as a function of hypervigilance, widespread pain index, symptom severity, or ACRFS total score following FWE correction with NBS. However, connection strength within one network motif, comprised of basal ganglia and cingulate ROIs, was observed to increase with symptom severity score when assessed across all participants. In particular, participants with higher symptom severity scores tended to exhibit stronger connections between rostral anterior cingulate and the pallidum (R rostral ACC [RH rACC]-L pallidum, $t=+4.45,P_{unc}=0.0002$, L rostral ACC [LH rACC]-L pallidum, $t=+4.12,P_{unc}=0.0004$) and putamen (RH rACC-L putamen, $t=+4.02,P_{unc}=0.0005$, LH rACC-R putamen, $t=+4.01,P_{unc}=0.0006$, LH rACC-L putamen, $t=+3.86,P_{unc}=0.0008$, RH rACC-R putamen, $t=+3.79,P_{unc}=0.0009$). These findings are presented as Figure 4B and Table 2A.

Group×condition adjacency matrices, “core” network, and graph theory analysis

Adjacency matrices for each group and condition were generated from suprathreshold connections following one-sample t-tests and NBS thresholding (Supplementary Figure S2). Connections present across groups and conditions were considered to constitute a “core” network, which was generated by a binary operation (elementwise multiplication) between the control-deflated, control-inflated, TMD-deflated, and TMD-inflated adjacency matrices (Supplementary Figure S2I) and rendered as Supplementary Figure S3 and the inset of Figure 5. Nodal communities within the core network, identified using the Louvain community detection algorithm, consisted of eleven classes: (I) hippocampus, anterior and posterior parahippocampal cortex, anterior and left posterior temporal fusiform cortex, and left intracalcarine cortex; (II) amygdala, entorhinal cortex, orbitofrontal cortex, and left superior lateral occipital cortex; (III) medulla, brainstem, periaqueductal grey, basal ganglia, subcallosal cortex, nucleus accumbens, anterior and left dorsal posterior insula, thalamus and hypothalamus, medial prefrontal cortex, the bank of the right superior temporal sulcus, and right middle temporal gyrus; (IV) pons, midbrain, superior cerebellar peduncle; (V) dorsolateral prefrontal cortex, right inferior frontal gyrus, frontopolar cortex, right occipital fusiform gyrus, and cerebellum; (VI) right caudal middle frontal gyrus, pre- and postcentral gyri, and supplementary motor cortex; (VII) right posterior temporal fusiform and ipsilateral intracalcarine cortex; (VIII) posterior cingulate and right superior lateral occipital cortex; (IX) right dorsal posterior insula, anterior and posterior supramarginal gyri, angular gyrus, parietal operculum, planum polare, and transverse temporal gyrus; (X) rostral and caudal anterior cingulate; and (XI) left occipital fusiform and left and right temporo-occipital fusiform cortex and the inferior lateral occipital cortex (Figure 5). Nodal communities were also detected for each group × condition adjacency matrix and are presented as Supplementary Figures S4, S5.

Finally, nodal degree and betweenness centrality were computed for each group × condition adjacency matrix and are presented as Supplementary Figure S6. Network “hubs,” identified by betweenness centrality rankings, included posterior cerebellum (betweenness centrality = 0.1774), right parietal operculum (0.1479), left anterior parahippocampal cortex (0.1372), and left putamen (0.1262) for controls under the deflated condition; left anterior insula (0.2013), right frontopolar cortex (0.1591), left entorhinal cortex (0.1360), left putamen (0.1349), right planum polare (0.1332), and left posterior supramarginal gyrus (0.1286) for controls under the inflated condition; left posterior parahippocampal cortex (0.0732), right frontopolar cortex (0.0721), right parietal operculum (0.0587), posterior cerebellum (0.0486), and left nucleus accumbens (0.0418) for TMD participants under the deflated condition; and right anterior insula (0.0871), right caudate (0.0776), right frontopolar cortex (0.0615), left posterior parahippocampal cortex (0.0557), and right parietal operculum (0.0507) for TMD participants under the inflated condition. These results suggest that the importance or “hub-ness” of the pons, right anterior insula, left parahippocampal cortex, left angular gyrus, and orbitofrontal cortex, among others, is amplified in TMD participants across conditions, despite lower centralities within the TMD group overall (Supplementary Figure S7).

Differences in either degree or centrality were separately assessed for each “ROI family” (CN V-TC, PAG-RVM, LPS, MPS, or DMN) by linear mixed-effects models. Across ROI families and conditions, TMD participants exhibited lower betweenness centralities [mean ± standard deviation in the deflated condition 0.0173 ± 0.0166 versus 0.0346 ± 0.0410 in controls; inflated, 0.0175 ± 0.018 versus 0.0345 ± 0.042 in controls; $F\left(1,306\right)=22.558,$ mean squared error (MSE) $=0.023,t\left(306\right)=-4.75,p<0.001,$ $R_{adjusted}^2=0.066,{\eta }_{partial}^2=0.069$, Cohen's $d=0.541$] and higher nodal degree [deflated, 9.674 ± 3.981 versus 4.558 ± 2.655 in controls; inflated, 10.116 ± 4.143 versus 4.047 ± 2.415 in controls; $F\left(1,306\right)=215.276,MSE=10.794,t\left(306\right)=14.672,$ $p<0.001,R_{adjusted}^2=0.413,{\eta }_{partial}^2=0.413$, Cohen's $d=1.672$] relative to control participants (Figure 6). Neither degree nor centrality differed within groups as a function of condition (deflated versus inflated) when taken across ROIs [centrality $F\left(2,308\right)=0.0010$, effect of condition within group at inflated: $\beta =-0.017,z=-0.267$, standard error (SE) 0.064, $p=0.789$; degree $F\left(2,308\right)=10.816$, effect of condition within group at inflated: $\beta =5.935,z=0.899,SE=6.599,p=0.368$], and consequently, additional nested effects and interactions were not explored. These findings therefore indicate strong effects of group on both node centrality and degree but failed to discern any meaningful differences within groups with respect to nodal centrality or degree in the cuff-deflated or cuff-inflated conditions.

Discussion

This exploratory study represents an investigation of the network topology of systems and pathways subserving pain sensation, perception, discrimination, anticipation, learning and motivational functions, and affective function in control and TMD participants at rest and during evoked-pain conditions, with a principal objective of documenting the overall network structure as a function of group and condition, including identification of “hub” nodes and topological changes during induced pain in each group. An attempt was made to ensure that the major nociceptive pathways, including the CN V-TC, pPAG-RVM, LPS and MPS, and DMN, as well as components of the “pain neuromatrix,” such as hippocampus, fusiform cortex, striatum, frontopolar cortex, cerebellum, and ventral pallidum, were well-represented.

In the current cohort, TMD participants differed substantively from controls with respect to overall pain (ACRFS total score or “fibromyalgianess”) widespread pain index, sensitivity to aversive sensations (hypervigilance), and influence of pain on fatigue, cognition, and other quality-of-life aspects (symptom severity). By design, cuff pressure was equilibrated across individual participants, such that perceived pain intensity during the evoked-pain (cuff-inflated) condition was comparable between the two groups (Table 1).

Morphometric findings

Surface- and volume-based assessments of grey and white matter structure were conducted primarily to generate additional ROIs which differed as a function of group or ACRFS and Pennebaker scores. Whole-brain/cluster-based and ROI-based morphometric analyses were conducted. These analyses evinced decreased left superior temporal surface area and volume (largely localized to the supramarginal gyrus), decreased surface areas in the right caudal middle frontal, right supramarginal, and right bank of the superior temporal sulcus, and increased right inferior frontal and right precuneus surface areas in TMD participants, as well as negative correlations between ACRFS total and widespread pain scores and corticospinal/corona radiata white matter volume, and a negative correlation between hypervigilance scores and grey matter volume in the pre- and postcentral gyri, across all subjects (Figure 1).

The significance and direction of structural findings differs across studies depending on covariates and other factors, and in general, such findings are inconsistent across studies, possibly due to individual differences in pain duration and history (5). Several studies have reported significant volumetric differences between TMD and controls [e.g., (102), and as reviewed in (47)]; for example, chronic pain is associated with grey matter alterations in thalamus and primary somatosensory cortex (47), and maxillofacial pain tolerance, in particular, has been associated with decreased grey matter volume in the right inferior frontal, right anterior insular, right putamen, left ventroposterior and right ventrolateral thalamus, and right globus pallidus, and increased grey matter volume in the ventral pons and middle cerebellar peduncle (5, 25, 47, 103). Additionally, the current results replicate several previous findings, including a reduction in superior (24) and middle temporal (104), medial prefrontal (105), and inferior frontal (24) grey matter volume. The reduction in inferior frontal surface area, as well as the negative associations currently observed between ACRFS total and WPI scores and corticospinal tract volume, and between hypervigilance scores and pre- and postcentral grey matter volume, may also comport with previously published observations that orbitofrontal grey matter thickness and primary motor cortical thickness are negatively associated with pain unpleasantness and duration, respectively (25, 106). Higher ACRFS scores were recently shown to be associated with higher levels of clinical pain, functional limitation, and disability and also with reduced responsiveness to arthroscopy of the temporomandibular joint in TMD patients (55, 56). These central findings could help explain why the clinical pain picture is often more complex in patients with higher fibromyalgianess. Interpretation of these findings is complicated by the deficiency in theoretical and clinical understanding of the causal (or acausal) relationship between pain and neural microanatomy. As noted by Younger, it is possible that these structural differences represent pre-existing vulnerability to pain or its chronic effects, or, alternatively, are secondary, adaptive responses to chronic pain; however, neither surface-based nor volume-based morphometry facilitates evaluation of critical aspects of neurite structure or glial density which would assist in such interpretation (25). Moreover, diffusion-based analyses and subcortical morphometry (including cerebellar morphometry) were omitted from the present study, thereby precluding examination of group-wise differences in these metrics in the present context. This omission should be addressed in future studies.

Association between connectivity strength and pain scores

Although TMD and control participants in the present study differed substantively with respect to symptom severity, hypervigilance, WPI, and ACRFS total scores (Table 1), only symptom severity was associated with differences in the strength of functional connectivity when taken across all participants. Symptom severity is defined in the ACRFS as the impact of pain on fatigue, cognitive symptoms, the restorative nature of sleep, and other somatic symptoms over the previous 7 days and, in the current study, was positively associated with connectivity strength between the right rACC and contralateral pallidum as well as between the left rACC and left and right putamen (Figure 4B). In addition to their conventional role in the regulation of voluntary movement (39, 107, 108), both the putamen and the pallidum (particularly the ventral pallidum, a component of the dopaminergic reward system) subserve learning and motivational functions (33, 34, 40, 41). These functions may be effected, in part, via putamenal connections with association, sensorimotor, and limbic cortices, including the insula (39). The rACC is considered part of the MPS (5, 7, 8, 76) and a second-order perceptual component of the nociceptive-processing system (9), and is believed to subserve the affective and cognitive aspects of pain. The association between symptom severity and rACC-putamenal-pallidal connectivity in this study may therefore reflect the affective and motivational aspects of pain. Depression, which is a component of symptom severity, was recently found to be associated with connectivity between the prefrontal cortex, the insula, and other cortical regions in patients with chronic low back pain (109, 110).

In contrast to symptom severity, the ACRFS WPI and Pennebaker hypervigilance are related to the distribution of pain over the body, and the tendency to prioritize pain signals over other somatic or environmental cues (111, 112), respectively. Although these perceptual and attentional correlates should arguably be reflected in inter-regional connectivity, no effects of the WPI or hypervigilance measures on connectivity strength were detected, despite substantive differences in WPI and hypervigilance between groups. Kutch et al. (2017) found the WPI to be positively correlated with functional connectivity between the salience network and the somatosensory and motor cortex in patients with chronic pelvic pain (113). The lack of significant findings in the present study may be due to the small sample size of, or exclusion of key regions of interest from, the current study. Alternatively, increased hypervigilance and WPI may partially explain the groupwise differences in the resting (deflated) and evoked-pain (inflated) condition adjacency matrices (described below); however, this possibility was not explored.

Finally, it must be noted that symptom severity, WPI, hypervigilance and ACRFS total score reflect participant phenotypes and correspond to general aspects of pain severity and quality-of-life, rather than severity or perceived pain intensity during the MRI sessions. As cuff pressure and pain rating were equilibrated across participants in both the control and TMD groups, it is unlikely that differences in the strength of functional connectivity reflect differences in perceived pain intensity during the evoked-pain condition.

Differences in connectivity strength as a function of group and condition

ROI-ROI connectivity strength was not observed to differ in the current cohort either between the control or TMD groups (across conditions) or between conditions in the control group. However, TMD participants demonstrated stronger connectivities between the left and right thalamus, right anterior insula, and right orbitofrontal cortex, and weaker connectivities between right anterior insula and the right temporo-occipital fusiform cortex, in the inflated condition compared to the deflated condition (Figure 4A). The anterior insula is a component of the MPS which functions in pain anticipation and salience, and may serve to modulate aversive, pain-related sensations via anticipatory downregulation of nociceptive regions (5, 7, 76, 114–117), and, in the present study, was correlated with both right and left thalamus ROIs but anticorrelated with the right fusiform cortex in TMD participants during the evoked condition. The fusiform cortex is known to function in body and face recognition and differentiation (30), functions which may be of particular interest in maxillofacial pain syndromes such as TMD. More generally, however, it has been shown to subserve pain anticipation and perception (31) and contributes to the affective and motivational aspects of pain. Like the striatum, accumbens, and orbitofrontal and prefrontal cortices, the fusiform cortex is part of the dopaminergic reward system (27–29) and is activated by stimulation of Aδ fibers (32) and during postoperative pain (118). Finally, the fourth component of this thalamic-anterior insular-fusiform-orbitofrontal network, the orbitofrontal cortex (OFC), also subserves pain inhibition and, with the accumbens nucleus, ventral striatum, dorsolateral prefrontal cortex, and thalamus (among others), is associated with motivational functions via the dopaminergic “reward circuit” (27–29, 119). The association of the components of this system with pain anticipation and modulation may thus indicate top-down effects on face-specific representations in the fusiform gyrus in anticipation of the evoked pain condition. However, no causal analyses were performed in this exploratory study, and such analyses could serve as a subject for future research.

Graph-theoretic analyses

In the current study, group × condition adjacency matrices were based on one-sample tests to determine whether a given ROI-ROI connection was “present” (or suprathreshold) or “absent” in a given group and condition. The assessments of nodal degree and centrality, which are derived from these adjacency matrices, therefore differ from the correlation-based analyses of ROI-ROI connectivity strength described previously, not only with respect to approach but also to interpretation, despite the application of NBS thresholding to both one-sample and two-sample tests. An example is presented as Supplementary Figure S8, in which the normalized connectivity values of the L thalamus:R OFC connection, which differed in strength between conditions within TMD participants, are compared to those of the L PO (parietal operculum):L dpIns (dorsal-posterior insula) connection, which was “present” in the TMD-inflated adjacency matrix, but “absent” in the controls-inflated adjacency matrix, based on one-sample tests. In the paired-sample, L thalamus:R OFC case, TMD-deflated connectivities were negatively skewed, with a greater number of negatively-valued connections, whereas TMD-inflated connectivities were positively skewed, with a number of positively-valued, moderately negatively-valued, and highly negatively-valued connections, despite the Fisher normalization performed during data preprocessing. Consequently, the mean connectivity values in both conditions diverged from each other and from zero, yielding a significant result in the paired-samples test (see Supplementary Figure S8A and Figure 4A). Conversely, the L PO:L dpIns connectivity values in the inflated condition were slightly positively skewed for both control and TMD participants, indicating the presence of a slightly greater number of positively-valued connections in both groups, but only TMD participants' connectivities diverged (as a group) so substantively as to yield a suprathreshold one-sample test result (see Supplementary Figure S8B and Figure 3). One interpretation of this result is that, taken as a group, there was insufficient evidence to conclude that the correlation between the left parietal operculum and ipsilateral dorsal posterior insula ROIs was greater than the mean correlation among all ROIs ($z=0$) in the control participants during the inflated condition, whereas the statistical evidence did indicate a divergence from the mean ROI-ROI correlation in TMD participants. Conceptually, then, pairwise tests, which account for differences in two connectivity distributions, are useful for determining differences in ROI-ROI connectivity strengths between groups and conditions, whereas adjacency matrices, and the graph-theoretic metrics derived from them, are useful for describing and visualizing the network topology within each group and condition. Adjacency matrices are, however, subject to small variations in the distribution of connectivities due to the presence of a greater number of positive or negative ROI-ROI correlations in a group and condition, meaning that sample size is a strong determinant of the stability of the adjacency matrices; consequently, these results will need to be replicated with much larger samples than those used in the present, exploratory study.

For the current cohort of participants, network topology was generally conserved across groups and, to some degree, across conditions. Visual inspection of the adjacency matrices and renderings of this network in the deflated (Figure 2) and inflated conditions (Figure 3), as well as visual plots of the nodal community structure derived for connections which were present across groups and conditions (denoted in this paper as the “core” network: see Figure 5 and Supplementary Figure S3) provide an overall sense of its structure at rest and in the presence of evoked pain. Eight major and three minor functional subdivisions across subcortical areas and the prefrontal, temporal, cingulate, occipital, and parietal lobes, comprised the core network. The largest of these subdivisions was comprised of striatal, subcallosal (subgenual) cingulate, accumbens nuclei, bilateral anterior and left dorsal posterior insula, thalamus, and hypothalamus. Other major subdivisions included temporal and medial temporal structures and calcarine cortex; amygdala, orbitofrontal, and entorhinal cortex; the pons, midbrain, and the superior cerebellar peduncle; frontal/prefrontal and fusiform cortex and the cerebellum; middle frontal and somatomotor cortex; right dorsal posterior insula and parietal structures; rostral and caudal anterior cingulate; and temporo-occipital fusiform and inferior lateral occipital cortex. Notably, however, the density of connections in TMD participants across both conditions contributed to variations in community structure compared with controls (Supplementary Figures S4, S5).

Indeed, the connection density of TMD participants, across both the cuff-delated and cuff-inflated, or evoked pain, conditions is apparent throughout the current findings, including the number of connections “present” in the group × condition adjacency matrices (Figures 2, 3) and the marked decrease in nodal betweenness centrality (Figure 6 and Supplementary Figures S6, S7) and increase in nodal degree (Figure 6 and Supplementary Figure S6), in addition to the integration of nodes into functional communities. TMD nodal centralities and nodal degree were approximately 50% and 230%, respectively, of their observed values in control participants, rendering discrimination of “hub” nodes more difficult in the TMD group. For example, whereas control participants exhibited marked differences in node centrality between the two conditions, such that the primacy of the posterior cerebellum, right parietal operculum, and right nucleus accumbens were apparent in the cuff-deflated condition, and that of the left anterior insula, right prefrontal cortex, left orbitofrontal cortex, and other nodes in the cuff-inflated condition, fewer differences in node primacy were evident in TMD participants. Whereas the centrality of right anterior insula, right caudate, and right planum polare were markedly increased in TMD participants in the presence of evoked pain, for example, differences between conditions with respect to the centrality of the right frontopolar cortex, parietal operculum, and left entorhinal cortex were less substantive in this group (see Supplementary Figure S7; note the difference in y-axis scale). This loss of “hub-ness” is called “hub disruption” and has been documented in chronic pain syndromes such as migraine (120), back pain (121–124) and in rodent models of neuropathic pain (125, 126). Additionally, increased connectivity with respect to hubs such as the planum polare, occipital cortex, frontopolar cortex, and supramarginal gyri in episodic migraine (51, 127), and among default-mode components (medial prefrontal cortex, thalamus, posterior cingulate) (26) across various pain syndromes. Kaplan et al. (2019) identified the anterior insulae as being important hubs for processing in fibromyalgia patients, and they also discovered that hub membership of some pain-processing regions varied based on clinical pain intensity. The regional direction and specificity of hub disruption differs, however, depending on the nature of the pain syndrome and whether the data being analyzed is structural (128, 129) or functional, and it is as yet unclear whether such disruption is a contributor to, or a consequence of, chronic pain syndromes, and its interpretation may be region-specific. For example, Yin and colleagues have proposed that alterations in primary and supplementary motor cortices may be a function of “maladaptive neuroplasticity” (5), whereas other instances, such as aberrant connectivity between the prefrontal cortex, pregenual anterior cingulate, and amygdala [reviewed in (5)] and connectivity among default-mode components, which has been associated with pain rumination (26), are at least partially due to alterations in the interaction between nociceptive and attentional systems. Moreover, many of these studies also describe decreases in functional connectivity, such as decreased frontopolar activation (50) or reduced medial prefrontal-PAG connectivity (48), which were not reproduced in the present paper.

It should be noted that there are numerous methodological approaches to assessing functional connectivity, and such approaches may be more appropriate for some purposes than the NBS analysis of bivariate correlations used in the present paper. For example, recent research by He et al. (130) and Zhang et al. (49) have investigated the spectral power of low-frequency BOLD signals (fractional amplitude of low-frequency fluctuations, fALFF) and synchrony of spontaneous signals within brain areas (regional homogeneity, ReHo), respectively, to demonstrate improvements in local or regional function with gnathological treatment. These improvements were also correlated with subjective pain relief. Notably, these approaches are appropriate for evaluation of local brain activity, rather than distributed network activity, and should be viewed as complementary approaches with resting state.

There are several limitations to the current study which render it exploratory rather than conclusive. The most obvious of these is the small sample size, which is not only subject to sampling variability and affects statistical power, but also excludes the possibility of more interesting and powerful approaches to the characterization of participant phenotypes and network structure, such as clustering and graph similarity, and limits the present analyses to groupwise comparisons. Additionally, the FWER and FDR corrections for multiple comparisons in this small sample may be overly conservative and select for the largest effect sizes, limiting replicability (131). We attempted to address some of these concerns through the NBS approach, which focuses hypothesis tests on network topology over connection strength and enhances the contrast-to-noise ratio in smaller samples; however, the sample size also admits the possibility of model overspecification, particularly in the linear mixed effects modeling. Also, TMD duration was not included as a covariate due to the difficulty in accurately establishing the duration from self-report as well as the lack of a prior hypothesis regarding how this duration would affect the results. Nevertheless, this is an important consideration which will need to be addressed in a larger cohort. With most neuroimaging techniques in cross sectional studies, it is difficult to ascertain causation concerning relationships between symptoms and brain attributes. There is some evidence that perceived pain levels while at rest in the scanner are associated with functional connectivity (132, 133), meaning that comparisons made between individuals experiencing chronic pain while in the scanner and those experiencing no pain could be driving some group differences in resting state connectivity. Inclusion of a task (like inflation of the cuff) helps to determine causal changes in connectivity, but it does not solve the between-group comparison issues. Finally, it should be noted that network topology cannot be reliably ascertained using resting-state connectivity alone, but should be premised upon a combination of modalities, including structural analyses. Despite these limitations, however, it is our expectation that our thorough documentation of source ROIs, which were carefully selected to ensure adequate representation of major nociceptive systems as well as components of the ill-defined “pain neuromatrix,” will help facilitate future studies, including those which serve to replicate or refute our present findings.

Conclusion

This pilot study represents a first attempt to leverage graph-theoretic analyses to evaluate the topology of the various nociceptive networks and systems, including those subserving pain sensation, perception, and discrimination, internal state and body perception, affective/motivational, and anticipatory functions, as a unified whole. Its results, which must be interpreted in the context of the small sample size, suggest that first-order nociceptive and second-order cognitive, affectual, and attentional systems comprise a distributed network located primarily in prefrontal, temporal, parietal, and cerebellar areas. Furthermore, the findings suggest that, although the substrates of this network do not generally differ between control and TMD participants in resting or evoked-pain conditions, the TMD syndrome may be associated with an overall increase in functional connectivity even at rest. Future studies will be required to replicate these results, incorporate structural connectivity, and establish causal associations among the network nodes delineated.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving human participants were reviewed and approved by the University of Michigan IRB. The patients/participants provided their written informed consent to participate in this study.

Author contributions

JS, JA, CF, and DH contributed to the conception and design of the study. DH collected the data. JS analyzed the data with some contributions from other authors. JS wrote the first draft of the manuscript. CF and DH contributed sections. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National Institute of Dental and Craniofacial Research grants K99 and R00 (DE026810) awarded to DH. The data were collected using the K99 funds and analyzed using R00 funds.

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/fpain.2022.966398/full#supplementary-material.

References