Transcriptional and cellular signatures of cortical morphometric remodelling in chronic pain

Supplemental Digital Content is Available in the Text. Imaging transcriptomics bridge levels to connect genes, cell classes, and biological pathways to in vivo imaging correlates of cortical morphometric remodelling chronic pain.


Introduction
Chronic pain is a prevalent and highly debilitating condition 11,23,61 with a moderately strong genetic basis (heritability, h 2 16%-50%). 34,67 Treatment response to current standard treatments is complex and overall poor. 24 Identifying new potential targets for drug development has become over the years a priority for researchers and clinicians. However, therapeutic advances have been marred by our still poor understanding of the mechanisms underlying chronic pain. 85 Chronic pain has been increasingly recognized as a disorder of the brain. 12,22,39,66,78,89 Using different neuroimaging techniques, a rich body of evidence has shown that most chronic pain syndromes is associated with a number of spatially patterned structural and functional alterations in different cortical and subcortical regions and brainstem. 40,58,75,89 The aberrant functional 19,26,41,48 and structural 10,21,45,92 configuration of the brain network detected across patients with different chronic pain syndromes has given rise to the idea that this distributed pattern of brain changes might reflect disruption of large-scale brain networks comprising anatomically connected brain areas. 25 However, exploring this hypothesis further has been constrained by challenges in measuring anatomical connectivity in humans. 5 For instance, the use of diffusion-weighted imaging to estimate the connectivity of long-distance projections, such as those between hemispheres, remains challenging. 28 On the other hand, structural covariance analysis is not applicable at the single-subject level, and because it relies on group-based covariance, it requires large sample sizes to be reliably estimated. 5 Morphometric similarity (MS) mapping has recently emerged as a new approach to construct whole-brain anatomical networks for individual subjects, overcoming some of the methodological limitations highlighted above. 62,76,77 It quantifies the similarity between cortical regions for multiple magnetic resonance imaging (MRI) parameters measured in each area. 77 This metric has close associations with the cytoarchitectonic properties of the cortex and axonal connectivity between regions. 77 Morphometric similarity mapping is a reliable method that can capture interindividual differences in cognition 77 and clinical abnormalities in brain disorders. 52,62,76 However, its use for uncovering morphometric differences in the brain of patients with chronic pain syndromes remains untried.
Spatially diffuse correlates of chronic pain across cortical anatomy could arise from a host of biological changes in patients, such as altered neurotransmission and maladaptive synaptic plasticity 44,60,70,72 and neuroinflammation, 37,38,90,94 among others. However, most neuroimaging modalities are not sensitive to underlying molecular or transcriptional properties of brain tissue. Therefore, understanding how disease-related alterations at the microscopic transcriptional architecture might explain regional vulnerability to macroscopic brain abnormalities in chronic pain syndromes remains rather challenging.
Here, we sought to bridge these gaps by examining alterations in MS in 3 independent case-control studies of patients with chronic pain syndromes: knee osteoarthritis (OA), chronic low back pain (CLBP), and fibromyalgia (FM). Moreover, given the tight relationship between regional MS and gene expression, 77 we leveraged data from the Allen Human Brain Atlas (AHBA) (Fig. 1) to explore among potential molecular and cellular pathways that might explain regional vulnerability to MS changes during chronic pain.

Samples
We used structural MRI data from 3 prior case-control studies on knee OA, 86 CLBP, 56 and FM. 69 The full details of each sample (including inclusion and exclusion criteria) have been described in the original publications. The OA data set included 56 OA chronic pain patients (54% women, 58 6 6.96 years) and 20 agematched healthy control subjects (50% women, 58 6 6.65 years). The CLBP data set included 29 CLBP patients (59% women, 30.79 6 11.59 years) and 33 healthy controls (14 females, 31.18 6 9.65 years). The FM data set included 20 women with FM (46.4 6 12.4 years) and 20 age-matched healthy control women (42.1 6 12.5 years). The original studies were approved by the competent ethics assessment board of each respective host institution. All participants provided informed consent before enrolling the respective studies.

MRI data acquisition
We used high-resolution T1-weighted anatomical images acquired in 3T scanners. We provide a summary of the parameters used for T1-weighted anatomical images acquisition in each data set (as described in the original publications).

Chronic low back pain data set
Images were acquired in a Siemens 3.0 T Trio B whole-body scanner equipped with a 32-channel head coil with the following parameters: TR/TE 5 1900/2.52 milliseconds, flip angle 5 9˚, matrix 256 3 256, slice thickness 5 1 mm, and number of slices 5 176. 56

Morphometric similarity mapping
The T1-weighted MRI data from all participants were preprocessed using the recon-all command from FreeSurfer (version 6.0). 15 The surfaces were then parcellated using the 68 cortical regions of the Desikan-Killiany Atlas. 20 For each cortical region, we estimated 5 parameters: gray matter volume, surface area, cortical thickness, Gaussian curvature, and mean curvature. Each parameter was normalized for sample mean and standard deviation to account for variation in value distributions between the features. After normalization, MS networks were generated by computing the regional pairwise Pearson correlations in morphometric feature sets, yielding a 68 3 68 MS matrix ℳi for each participant, i 5 1, …, N, which represents the strength of MS between each pair of cortical areas. For all individuals, regional MS (ie, nodal similarity) estimates were calculated as the average MS between a given cortical region and all others. Although MS calculation allows for the inclusion of data from different imaging modalities, here, we used only features extracted from T1weighted MRI data. It has been previously demonstrated that there is high spatial concordance (r 5 0.91) between the topography of regional MS derived from T1-weighted MRI data alone and regional MS from more modalities (eg, a combination of T1-weighted and diffusion-weighted imaging). 43,77

Case-control analysis of morphometric similarity networks
The global MS for each participant is the average of ℳ i . The regional MS i,j for the ith participant at each region j 5 1,…,68 is the average of the jth row (or column) of ℳ i . Thus, regional MS is equivalent to the weighted degree or "hubness" of each regional node, connected by signed and weighted edges of pair-wise similarity to all other nodes in the whole-brain connectome represented by the MS matrix. For global and regional MS statistics alike, we fit linear models to estimate case-control differences, with age, sex, and total intracranial volume as covariates. The resulting P value for each region was thresholded for significance using false discovery rate (FDR) , 0.05, to control type 1 error over multiple (n 5 68) tests. To test for the differential contribution of single cortical features to the observed regional MS changes in each chronic pain condition, we recomputed condition-specific MS change maps with exclusion of each individual cortical feature before MS calculation and then determined which of these single-feature exclusions most change the topography of the observed MS changes. We did so by calculating pairwise Spearman correlations between the original and all other leave-one-out maps of MS changes, in each condition separately.
2.5. Cortical morphometric similarity remodelling in major depressive disorder and its association with chronic pain Chronic pain is often comorbid with major depressive disorder (MDD) 8 and these 2 entities share brain mechanisms of neuroplasticity. 80 We used another openly available data set (see the original article for further details 50 ) including high resolution structural brain data from 19 unmedicated patients with MDD (11 females, 29.45 6 11.26 years) and 20 age-matched healthy e760 D. Martins et al. · 163 (2022) e759-e773 PAIN ® controls (11 females, 33.52 6 14.07 years). We characterized changes in regional MS associated with MDD and investigated whether changes in regional MS in MDD can predict those we observed in chronic pain patients. We calculated MS using the procedures described above and tested for case-control differences between MDD patients and healthy controls using linear models, where we accounted for age, sex, and intracranial volume. Finally, we calculated pairwise Spearman correlations between changes in regional MS in MDD and those observed in the different chronic pain conditions. 2.6. Mapping case-control differences in regional morphometric similarity to established patterns of cytoarchitectonic cortical organization To help us to contextualize the case-control differences in regional MS we observed for the different chronic pain conditions, we mapped them in relation to well-established patterns of cytoarchitectonic organization of the cortex. To that end, we used the von Economo Atlas of the cortex classified by the cytoarchitectonic criteria. 95 For each subject, we quantified MS within each parcel of these atlases and then performed casecontrol comparisons using linear models, with age, gender, and intracranial volume as covariates.

Morphometric similarity hubs susceptibility
We investigated relationships between case-control differences in regional MS and the typical pattern of regional MS distribution in healthy controls with Spearman correlations. In keeping with histological results indicating that cytoarchitectonically similar areas of the cortex are more likely to be anatomically connected and that MS in the macaque cortex was correlated with tracttracing measurements of axonal connectivity, 77 we followed the approach suggested by Seidlizt et al. 76 to map each region to one of 4 patterns of changes in MS: (1) regions of low MS in healthy controls (highly differentiated from the rest of the cortex) that increase their MS during with the rest of the cortex during chronic pain (dedifferentiation), (2) regions of high MS in healthy controls (highly connected with the rest of the cortex) that increase their MS during chronic pain (hypercoupling), (3) regions of low MS in healthy controls that decrease their MS during chronic pain (hyperdifferentiation), and (4) regions of high MS in healthy controls that decrease their MS during chronic pain (decoupling). We subdivided each axis of the scatter plot in 2 sections, one above and another below 0, which resulted in 4 quadrants, each representing one of the 4 scenarios presented above. We then quantified the percentage of regions falling within each of these 4 quadrants to identify dominant patterns of change.
2.8. Defining a cross-condition pattern of changes in regional morphometric similarity during chronic pain We investigated the similarity in case-control changes in regional MS across chronic pain conditions by calculating pairwise Spearman correlations of regional case-control statistics (Zscores) from each condition. 96 We found significant correlations for all possible pairs of conditions, which indicated that remodelling of regional MS during chronic pain shares a common pattern across conditions. We then ran principal component analysis on the 3 vectors of case-control changes in regional MS to identify this shared profile of cross-condition changes. The first component alone explained 64.45% of the shared variance, the second 25.38%, and the third 10.16%. Only the first component showed an eigenvalue . 1 (1.93). Hence, we kept only the first PC because case-control changes in regional MS across our 3 chronic pain conditions seem to be well captured by one single dominant cross-condition pattern.
2.9. Microarray expression data: Allen Human Brain Atlas Regional microarray expression data were obtained from 6 postmortem brains provided by the AHBA (AHBA; http://human.brainmap.org/) (aged 24-57 years). 33 We used the abagen toolbox (https://github.com/netneurolab/abagen) to process and map the transcriptomic data to 84 parcellated brain regions from the Desikan-Killiany Atlas. 20 In brief, genetic probes were reannotated using information provided by Arnatkeviciute et al., 7 instead of the default probe information from the AHBA data set, hence discarding probes that cannot be reliably matched to genes. Following previously published guidelines for probe-to-gene mappings and intensity-based filtering, 7 the reannotated probes were filtered based on their intensity relative to the background noise level; probes with intensity less than background in $50% of samples were discarded. A single probe with the highest differential stability, DS(p), was selected to represent each gene, where differential stability was calculated as 32 : Here, r is Spearman's rank correlation of the expression of a single probe p across regions in 2 donor brains, B i and B j , and N is the total number of donor brains. This procedure retained 15,633 probes, each representing a unique gene.
Next, tissue samples were assigned to brain regions using their corrected MNI coordinates (https://github.com/chrisfilo/alleninf) by finding the nearest region within a radius of 2 mm. To reduce the potential for misassignment, sample-to-region matching was constrained by hemisphere and cortical or subcortical divisions. If a brain region was not assigned to any sample based on the above procedure, the sample closest to the centroid of that region was selected to ensure that all brain regions were assigned a value. Samples assigned to the same brain region were averaged separately for each donor. Gene expression values were then normalized separately for each donor across regions using a robust sigmoid function and rescaled to the unit interval. Scaled expression profiles were finally averaged across donors, resulting in a single matrix with rows corresponding to brain regions and columns corresponding to the retained 15,633 genes. As a further sanity check, we conducted leave-one-donor out sensitivity analyses to generate 6 expression maps containing gene expression data from all donors, one at a time. The principal components of these 6 expression maps were highly correlated (average Pearson correlation of 0.993), supporting the idea that our final gene expression maps where we averaged gene expression for each region across the 6 donors is unlikely to be biased by data from a specific donor. Because the AHBA only includes data for the right hemisphere for 2 subjects, in our transcriptomic-imaging association analyses, we only considered the left hemisphere cortical regions (34 regions).

Identifying transcriptomic correlates of cortical morphometric similarity remodelling in chronic pain
To be able to investigate associations between cross-condition changes in MS during chronic pain and brain gene expression, we used partial least square regression (PLS). 62 Partial least square regression uses the gene expression measurements (the predictor variables) to predict the regional MS changes (the response variables). This approach allows us to rank all genes by their multivariate spatial alignment with cross-condition regional MS changes during chronic pain. The first PLS component (PLS1) is the linear combination of the weighted gene expression scores that have a brain expression map that covaries the most with the map of MS changes. As the components are calculated to explain the maximum covariance between the dependent and independent variables, the first component does not necessarily need to explain the maximum variance in the dependent variable. However, as the number of components calculated increases, they progressively tend to explain less variance in the dependent variable. Here, we tested across a range of components (between 1 and 15) and quantified the relative variance explained by each component. The statistical significance of the variance explained by each component was tested by permuting the response variables 1000 times.
In our analysis, a solution with a single component explained variance in regional MS changes above chance (P boot 5 0.003). The first PLS component (PLS1) alone explained the highest amount of variance alone (24.42%). Hence, we focused our further gene set enrichment analyses on PLS1. The error in estimating each gene's PLS1 weight was assessed by bootstrapping (resampling with replacement of the 34 brain regions), and the ratio of the weight of each gene to its bootstrap standard error was used to calculate the Z scores and, hence, rank the genes according to their contribution to PLS1. 95 Genes with large positive PLS1 weights correspond to genes that have higher than average expression in regions where MS increases and lower than average expression in regions where MS decreases. Midrank PLS weights showed expression gradients that are weakly related to the pattern of MS changes. On the other side, genes with large negative PLS1 weights correspond to genes that have higher than average expression in regions where MS decreases the most and lower than average expression in regions where MS increases. Hence, from the ranked PLS1 list of genes, we then selected all genes with positive and negative weights Z . 3 and Z , 23, respectively (all P FDR , 0.05, FDR corrected for the total number of genes tested). PLS1 genes with Z . 3 are for simplicity termed PLS11 and genes with Z , 23 PLS12. Although our choice of Z 5 3 as a threshold to identify the most positively and negatively genes associated with MS changes in patients chronic pain is somehow arbitrary, we note that Z 5 3 in our case corresponds to a stringent threshold of P FDR 5 0.036 (more stringent than P FDR , 0.05). We used these 2 sets of genes for further enrichment analyses, as described below. We confirmed our enrichment analyses were not driven by the choice of this specific threshold by repeating all analyses described below considering all genes that passed a more liberal threshold of P FDR , 0.05. The overall pattern of results did not change.

Protein-protein networks and gene set enrichment analysis
We then used all genes in PLS11 and PLS12 to conduct further bioinformatics analyses investigating whether these genes map to common and relevant biological pathways. First, we used STRING (version 10.5) 84 to construct protein-protein functional interactions networks. We excluded text mining as an active interaction source and used the default medium required For each region, we averaged across all the edges involving that area to obtain a singular representation of the mean MS score for that region. We then computed case-control differences for each region, while accounting for age, gender, and total intracranial volume. (B) Gene expression analysis. We used abagen to obtain gene expression profiles from the Allen Human Brain Atlas (AHBA) in 68 regions (left hemisphere only) across the 6 postmortem brains sampled in this atlas. We excluded all genes with normalized expression values below the background (15,633 genes met this criterion). When more than one probe was available for a certain gene, we selected the probe with higher consistency in expression across the 6 donors. We used partial least squares regression (PLS) to rank all genes according to their association with the case-control changes in MS. Finally, we performed a set of enrichment analyses on the top genes positively or negatively associated with case-control differences in MS. interaction score of 0.400 to identify all possible links within our list of target genes. Second, we used the GENE2FUNC function from the Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA) 93 platform to investigate functional enrichments using rank-based Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathways enrichment analysis. We used as background all AHBA genes that passed our preprocessing criteria and hence were used in our PLS analyses (n 5 15,633). We corrected for multiple gene-set enrichment testing by applying FDR correction.

Brain cell-type enrichment analysis
We also investigated whether our PLS11 and PLS12 subsets of genes were particularly enriched for genes of specific brain cell types. We compiled data from 5 different single-cell studies using postmortem cortical samples in human postnatal subjects to avoid any bias based on the acquisition methodology or analysis or thresholding. 16,31,47,51,99 To obtain gene sets for each cell type, categorical determinations were based on each individual study, as per the respective methods and analysis choices in the original article. All cell-type gene sets were available as part of the respective articles. We generated a single omnibus gene list for each cell type by merging the study-specific gene lists and then filtered it to retain only genes sampled in the AHBA. Two studies did not subset neurons into excitatory and inhibitory, 16,99 and thus, those gene sets were excluded from the cell-class assignment. In addition, only one study included the annotation of the "Per" (pericyte) type, and thus, we did not consider that cell type. 47 This approach has already been validated elsewhere. 76 We then conducted cell-class enrichment analyses using the GeneOverlap package from R (version 1.26.0). We used as background all AHBA genes that passed our preprocessing criteria and hence were used in our PLS analysis (n 5 15,633). GeneOverlap calculates the overlap between 2 sets of genes (in our case, the set of PLS11 or PLS12, and each of the brain cell-type omnibus gene sets derived as explained above) and uses a Fisher exact test to find whether the overlap between these 2 sets is higher than one would expect by randomly selecting a subset of genes from the background with the same number of elements. Here, enrichment is quantified as an odds ratio, where values lower or equal to 1 indicate minimal overlap between sets and hence absence of enrichment. Therefore, the null hypothesis is that the odds ratio is no larger than 1. Significant odds ratio larger than 1 indicate enrichment for genes of a specific cell type. We applied FDR correction for the number of cell types tested.

Pain-related and other brain disorder-related genes enrichment
We also investigated whether our PLS11 and PLS12 subsets of genes were particularly enriched for pain-related and other brain disorder-related genes. The list of previously identified painrelated genes was defined using the following public resources: (1) the pain-related genes identified in mice gene knockout studies collected in the Pain Gene Database (http://paingeneticslab.ca/4105/06_02_pain_genetics_database.asp), 46 (2) the pain-related genes identified in humans collected in the Human Pain Genes Database (https://humanpaingenetics.org/hpgdb), 59 (3) the Pain Research Forum (https://www.painresearchforum. org/resources/pain-gene-resource), 17 (4) the genes involved in human pain diseases collected in the DisGeNET (http://www. disgenet.org), 71 and (5) the genes considered to be functioning in pain perception summarized in GO (GO:0019233). In total, we identified 2111 pain-related genes included in the above public resources. Eight hundred seven of these genes were not part of our initial list of 15,633 AHBA genes and were excluded from further analyses (please see Supplementary data S3 for the list of 1304 pain-related genes used in these analyses, available at http://links.lww.com/PAIN/B511). The lists of genes associated with other brain disorders (Alzheimer disease, Parkinson disease, Huntington disease, epilepsy, autism spectrum disorder, MDD, anxiety, bipolar disorder, and schizophrenia) were collected from the DisGeNET. We tested for gene enrichment in both PLS11 and PLS12 subsets using GeneOverlap, as explained above.

Spatial permutation test (spin test)
In several analyses in the current study, we investigated the spatial correspondence between different imaging-derived measures. Although several studies have reported significance based on the assumption that the number of samples is equal to the number of regions, this is technically inaccurate, as the number of regions is both arbitrary (due to the resolution of the chosen parcellation) and non-independent (because of spatial autocorrelation among neighbouring parcels). To overcome this issue, we used spatial permutation tests (spin test) as implemented in previous studies. 5,6,91 This approach consists in comparing the empirical correlation among 2 spatial maps to a set of null correlations, generated by randomly rotating the spherical projection of one of the 2 spatial maps before projecting it back on the brain surface. Importantly, the rotated projection preserves spatial contiguity of the empirical maps, as well as the hemispheric symmetry. Therefore, each analysis correlating values from 2 cortical maps is reported with a P-value derived from the spherical permutation (P spin ), obtained by comparing the empirical Spearman Rho to a null distribution of 10,000 Spearman correlations, between one empirical map and the randomly rotated projections of the other map. The Matlab code to implement this test can be found in https://github.com/ frantisekvasa/rotate_parcellation.

Code availability
The code for MS and gene expression association analyses is available at https://github.com/SarahMorgan/Morphometric_ Similarity_SZ.

Morphometric similarity is a reproducible measure in healthy controls across independent studies
The cortical maps of regional MS in Supplementary Figure S1 summarizes the anatomical distribution of areas of positive and negative MS in healthy controls from each of the 3 data sets (available at http://links.lww.com/PAIN/B508). The patterns of regional distribution of MS were highly correlated across healthy controls from the 3 data sets (Supplementary Figure S2, available at http://links.lww.com/PAIN/B508). The results are similar to those reported in other independent samples, 62,76,77 with positive MS in frontal and temporal cortical areas (indicating high levels of similarity with the rest of the cortex) and negative MS in occipital and cingulate cortices (indicating low levels of similarity with the rest of the cortex; hence, high levels of differentiation from the rest of the cortex). This confirms the replicability of this pattern of regional MS in healthy individuals.

Patients with chronic pain do not differ from healthy controls in global morphometric similarity
Regional MS had an approximately normal distribution over all 68 cortical regions (after regressing age, sex, and intracranial volume) in both patients and healthy controls from all 3 data sets ( Fig. 2-upper panel). We did not find any significant casecontrol differences in global MS in any of the 3 data sets (OA: T(75)  (Fig. 2-lower panel).
3.3. Differences in regional morphometric similarity in patients with chronic pain syndromes as compared with healthy controls The cortical maps in Figure 3-upper panel summarize the distribution of case-control changes in cortical MS for each chronic pain condition. In Figure 3-lower panel, we highlight only regions with case-control differences significant at P , 0.05, uncorrected (none of these regions survived FDR correction). In the OA data set, we found decreases in MS in the left superior frontal gyrus, right pericalcarine cortex and in the left posterior cingulate, and increases in the left insula and inferior temporal gyrus, and in the right bank of the superior temporal sulcus and right inferior temporal gyrus. In the CLBP data set, we found decreases in MS in the left and right superior parietal gyri and left lateral occipital cortex; and increases in the left entorhinal cortex and caudal anterior cingulate, and in the right insula. In the FM data set, we found decreases in the left superior parietal, medial, and inferior temporal and fusiform gyrus and increases in the left and right isthmus of the cingulate, left posterior cingulate, and entorhinal and parahippocampal cortices. Changes in regional MS correlated positively between conditions (Supplementary Figure S3, available at http://links.lww.com/PAIN/B508), suggesting the existence of a shared pattern of regional MS changes across the 3 conditions. To further support the existence of this pattern, we performed a principal component analysis on the regional MS changes of the 3 conditions, finding that the first PC explained most variance (64.45%) in case-control changes across the 3 conditions (PC1).
Because we used 5 different cortical features to estimate MS (gray matter volume, surface area, cortical thickness, Gaussian curvature and mean curvature), we tested for differential contributions of single cortical features to the observed regional MS changes by recomputing the condition-specific MS change maps excluding each individual cortical feature at a time before MS calculation. We then determined which of these singlefeature exclusions most changed the topography of the observed MS changes by identifying the leave-one-feature-out map that correlated the least with the total map. This leave-one-out procedure showed that cortical thickness was the feature that most contributed to the topography of observed regional MS changes across the 3 conditions ( Supplementary Figures S4 and  S5, available at http://links.lww.com/PAIN/B508). Yet, all leaveone-out maps were positively correlated in all the 3 conditions (Supplementary Figure S4, available at http://links.lww.com/ PAIN/B508).
Chronic pain is often comorbid with MDD. 8 Several preclinical and clinical studies have found considerable overlaps between chronic pain-induced and depression-induced neuroplasticity. 80 Figure 2. Case-control differences in global morphometric similarity. In the upper panel, we present case and control distributions of the residual regional morphometric similarity (MS) strength (ie, the average similarity of each region to all other regions) after regressing out the effects of age, gender, and intracranial volume, for each data set separately. In the lower panel, we present case-control comparisons of the global MS. To calculate global MS, we averaged the residual regional MS strength across all regions for each subject. In all 3 data sets, there were no differences between groups in global MS. A, osteoarthritis dataset; B, chronic low back pain dataset; C, fibromyalgia dataset. CLBP, chronic low back pain; HC, healthy controls; OA, osteoarthritis. To investigate whether the regional MS changes we report here might predominantly reflect neuroplasticity associated with mood alterations other than chronic pain per se, we used another openly available data set including high resolution structural brain data from unmedicated patients with MDD and healthy controls. 50 We used these data to define the pattern of changes in regional MS associated with MDD and to investigate whether changes in regional MS in the MDD data set can predict those we observed during chronic pain. We found a considerably different pattern of MS changes in unmedicated MDD patients, as compared with never depressed healthy controls (Supplementary Figure S6A, available at http://links.lww.com/PAIN/B508). Morphometric similarity changes in the MDD data set did not correlate with MS changes in OA and correlated negatively with MS changes in CLBP and FM. We also found a negative correlation between MS changes in MDD and PC1 capturing the crosscondition pattern of changes (Supplementary Figure S6B, available at http://links.lww.com/PAIN/B508). Hence, the pattern of MS changes we report here for patients with chronic pain is unlikely to simply reflect neuroplastic changes associated with comorbid MDD.

Mapping case-control differences in regional morphometric similarity to established patterns of cytoarchitectonic cortical organization
To help us contextualize the case-control differences in regional MS we observed for the different chronic pain conditions, we mapped them onto well-established patterns of cytoarchitectonic organization of the cortex as defined by the von Economo Atlas of the cortex. 95 In OA, we found an increase in MS in the insular cortex (P , 0.05, uncorrected). In both CLBP and FM, we found increases in MS in the limbic cortex (P , 0.05, uncorrected). We also found decreases in association cortex A in FM (P , 0.05, uncorrected) (Fig. 4B and Supplementary Table S1, available at http://links.lww.com/PAIN/B508). None of these changes survived FDR correction.
3.5. "Hub susceptibility": associations between case-control differences in regional morphometric similarity and regional morphometric similarity in healthy controls Previous studies have shown that highly connected "hub" regions are the most likely to show reduced connectivity in disease as measured in functional magnetic resonance imaging and diffusion tensor imaging (DTI) brain networks. 14 Given the tight positive association between axonal connectivity and MS, these highly connected "hub" regions map to regions with constitutive high MS. Hence, here we also tested this "hub susceptibility" model by investigating relationships between case-control changes in regional MS and the typical pattern of regional MS distribution in healthy controls. Across the 3 data sets, we found significant negative correlations between these 2 variables (OA-Spearman Rho 5 20.397, P spin 5 0.001; CLBP-Spearman Rho 5 20.379, P spin 5 0.003; FM-Spearman Rho 5 20.522, P spin 5 1.7 3 10 204 ) (Fig. 5-upper panel). Moreover, we found that most regional increases map to regions typically showing low MS in controls, whereas most regional decreases map to regions of high MS in controls (Fig. 5-lower panel). Altogether, these findings support the idea that chronic pain is associated with decoupling of MS "hubs" and dedifferentiation of highly differentiated regions.

Transcriptomic correlates of cortical morphometric similarity remodelling during chronic pain
We investigated associations between cross-condition changes in MS during chronic pain (PC1) and brain gene expression using PLS. The first PLS component (PLS1) explained the highest proportion of MS changes (24.4%) and did so above chance Figure 3. Case-control differences in regional morphometric similarity. In this figure, we present cortical maps of the distribution of Cohen d effect size quantifying case-control differences in regional morphometric similarity (D morphometric similarity) for each condition. In the upper panel, we present unthresholded maps. In the lower panel, we present only regions where we found case-control differences for a threshold of P , 0.05, uncorrected. Note that none of these regions survived FDR correction for the total number of regions tested within each condition. A, osteoarthritis dataset; B, chronic low back pain dataset; C, fibromyalgia dataset. FDR, false discovery rate.
June 2022 · Volume 163 · Number 6 www.painjournalonline.com e765 (P boot 5 0.003). PLS1 gene expression weights were positively correlated with cross-condition changes in regional MS (r 5 0.494, P spin 5 0.013) (Fig. 6, panel A). This positive correlation means that genes positively weighted on PLS1 are highly expressed in regions where MS was increased, whereas negatively weighted genes are highly expressed in regions where MS was decreased in patients. We found 338 genes with PLS1 weights Z . 3 (which we denoted PLS11) and 236 genes with Z , 23 (PLS12) (all P FDR , 0.05). The gene with the highest positive weight was the "Chemokine-Like Factor Superfamily Member 3" (CMTM3), a microglia gene related to immune cytokine activity. The gene with the lowest negative weight was the "Family with Sequence Similarity 126, Member B" (FAM126B), which is a part of a complex required to localize phosphatidylinositol 4-kinase to the plasma membrane (Fig. 6B).

Protein-protein networks and gene set enrichment
We mapped the network of known interactions between proteins coded by the PLS11 and PLS12 gene sets (Supplementary Figure S7, available at http://links.lww.com/PAIN/B508). For PLS11, the resulting network had 303 nodes and 615 edges, more than the 239 edges expected by chance (PPI enrichment Pvalue , 1 3 10 216 ). Using Gene Set Enrichment Analysis, we found enrichment for a number of GO terms-biological pathways broadly mapping to the neuroimmune response axis (Fig. 7A). For PLS12, the resulting network had 225 nodes and 195 edges, more than the 152 edges expected by chance (PPI enrichment P-value 0.0005). Using Gene Set Enrichment Analysis, we did not find enrichment for any GO terms but found 4 terms from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that reached significance. Those were "Calcium signalling," "Long-term potentiation (LTP)," "Taste transduction" and "Type II diabetes mellitus" (Fig. 7A).

Enrichment for transcriptional signatures of canonical brain cell types
We also performed cell-type enrichment analysis using omnibus lists of gene expression in different brain cells of the postmortem human brain, as characterized across 5 different studies. For PLS11, we found significant enrichment in genes typically expressed in microglia, astrocytes, and oligodendrocytes precursor cells (Fig. 7B and Supplementary Table S4, available at http://links.lww.com/PAIN/B508). For PLS12, we found significant enrichment in genes typically expressed in excitatory and inhibitory neurons (Fig. 7B and Supplementary Table S2, available at http://links.lww.com/PAIN/B508).

Enrichment for genes related to pain and other brain disorders
Finally, we investigated whether PLS11 and PLS12 are enriched for pain-related and other brain disorder-related genes as identified in previous studies. For PLS11, we found enrichment for pain-related genes(OR 5 1.40, P 5 0.04), but not genes related to any of the other brain disorders we tested (Supplementary Table S3, available at http://links.lww.com/PAIN/B508). For PLS12, we did not find enrichment for pain-related genes (OR 5 1.07, P 5 0.42) but found enrichment for genes associated with epilepsy (OR 5 1.54, P 5 0.02, P FDR 5 0.09) and MDD (OR 5 1.57, P 5 0.02, P FDR 5 0.09) (Supplementary Table S3, available at http://links.lww.com/PAIN/B508).

Discussion
In this article, we uncovered a new pattern of cortical MS remodelling across 3 chronic pain syndromes, which was different from that observed in patients with MDD and points towards the existence of shared disease mechanisms driving cortical remodelling that cut across the boundaries of specific . Mapping case-control differences in regional morphometric similarity to established patterns of cytoarchitectonic cortical organization. To help us to contextualize the case-control differences in regional morphometric similarity we observed for the different chronic pain conditions, we mapped them in relation to the von Economo Atlas of the cortex classified by the cytoarchitectonic criteria. For each subject, we quantified MS within each parcel of these atlases (after regressing out age, sex, and intracranial volume) and then performed case-control comparisons using independent sample t tests. The colours in the tile plots represent the Cohen d effect size of case-control differences for each condition. *Highlights significant case-control differences, for P , 0.05, uncorrected. MS, morphometric similarity. pain syndromes. Furthermore, we demonstrate that cortical MS remodelling in chronic pain spatially correlates with the brain-wide expression of pain-related genes and genes involved in glial immune response and neuronal plasticity, which links putative underlying molecular perturbations with regional vulnerability to brain structural changes. These findings bridge levels to connect genes, cell classes, and biological pathways to in vivo imaging correlates of chronic pain and provide food for beleieved in how future treatment development might be pursued.
Morphometric similarity mapping disclosed a pattern of cortical MS changes across chronic pain syndromes, which generally involved small-to-medium-sized increases in the insula and limbic cortices, and decreases in the occipital, sensorimotor, and frontal cortices. Although most of these changes did not survive correction for multiple comparisons and should be interpreted cautiously, these findings are interesting for several reasons. First, of the various brain regions that have been implicated in the perception of pain, the insula and limbic system are among the ones most consistently reported across studies. 35 Second, functional and structural alterations in these regions have often been reported in neuroimaging studies of patients with different chronic pain syndromes, 40,58,75,89 including increases in the connectivity of the insula with nodes of the default-mode network. 9,55,65 Third, previous studies have further demonstrated that the limbic system plays a key role in the transition from acute-to-chronic pain. 56,57 Altogether, these aspects reinforce the neuroanatomical plausibility of the pattern in regional MS changes we report in this article.
Morphometric similarity quantifies the correspondence or kinship of 2 cortical areas for multiple macrostructural features that are measurable by MRI. 77 Hence, high MS between a pair of cortical regions indicates that there is a high degree of correspondence between them for cytoarchitectonic features. This assumption has received empirical support in prior work showing that morphometrically similar cortical regions share patterns of gene coexpression and are more likely to be axonally connected to each other. 77 Therefore, here, we interpret reduced MS as indicating that there is reduced cytoarchitectonic similarity, or greater cytoarchitectonic differentiation, between these areas and the rest of the cortex, which is probably indicative of reduced anatomical connectivity to and from the less similar, more differentiated cortical areas. On the other hand, increased MS implies increased cytoarchitectonic similarity and, perhaps, axonal connectivity with the rest of the cortex. In line with previousfunctional magnetic resonance imaging and DTI brain networks studies have demonstrated that "hub" regions are more likely to be disturbed and reduce their connectivity in the presence of brain disease, 14 we found negative associations between regional MS in healthy controls and MS changes in patients with chronic pain. Therefore, it is not implausible that the decreases in MS during chronic pain we describe here might reflect an overall pattern of decreases in axonal connectivity of "hub" regions with the rest of the cortex, as observed in other brain disorders. The reverse, ie, increased connectivity, might drive increases in MS during chronic pain. Nevertheless, we Figure 5. Chronic pain is associated with decoupling of morphometric similarity "hubs" and dedifferentiation of highly differentiated regions. In the upper panel, we present scatterplots depicting significant negative associations between case-control differences in regional morphometric similarity and the mean morphometric similarity of each region in healthy controls. In the lower panel, we subdivided these scatter plots in 4 quadrants according to the relationship between the distribution of case-control differences and mean regional morphometric similarity in healthy controls. The upper left quadrant includes regions of low MS in healthy controls (highly differentiated from the rest of the cortex) that increase their MS with the rest of the cortex during chronic pain (dedifferentiation). The upper right quadrant includes regions of high MS in healthy controls (highly connected with the rest of the cortex) that increase their MS during chronic pain (hypercoupling). The lower left quadrant includes regions of low MS in healthy controls that decrease their MS during chronic pain (hyperdifferentiation). The lower right quadrant includes region of high MS in healthy controls that decrease their MS during chronic pain (ecoupling). Across chronic pain conditions, we observed a predominant pattern of increases in regional MS in regions of low MS (dedifferentiation) and decreases in regions of high MS (decoupling) in healthy controls. Significance was assessed with spatial permutation testing. A, osteoarthritis dataset; B, chronic low back pain dataset; C, fibromyalgia dataset. MS, morphometric similarity.
June 2022 · Volume 163 · Number 6 www.painjournalonline.com e767 cannot exclude that either increases or decreases in MS might simply reflect local changes in cytoarchitectonics or even a combination of local tissue changes and connectivity.
In an attempt of connecting these MS changes during chronic pain to the gene expression and cellular pathways potentially explaining regional vulnerability to those changes, we used PLS to identify the weighted combination of genes in the whole transcriptome that has a cortical expression map most similar to the cortical map of cross-condition case-control MS differences we derived for patients with chronic pain. Further reinforcing the relationship of this subset of genes with pain, we found enrichment in PLS11 for pain-related genes but not genes related to other brain disorders. In PLS12, we did not find enrichment for pain-related genes but found enrichment for genes related to epilepsy and MDD. This last finding is in keep with the idea that chronic pain might share neurobiological pathways with epilepsy 18,68,83 and depression 80 and matches well with the clinical observation that antiepileptic drugs 81 and antidepressants 27 also improve pain in patients with chronic pain syndromes. We also characterized further the top genes positively and negatively associated MS changes in chronic pain by conducting agnostic gene set enrichment analyses and enrichment for genes associated with different classes of brain cell types. In PLS11, we found predominance of genes related to the glial immune response and highly expressed in microglia, astrocytes, and oligodendrocyte precursor cells; although in PLS12, we found predominance of genes related to calcium signalling and LTP, which are highly expressed in excitatory and inhibitory neurons. Altogether, these findings suggest that the constitutive distribution of genes involved in glial immune response and neuronal plasticity, both key elements of the current pathophysiological models of chronic pain, 37,101 can explain variance in the regional vulnerability to MS cortical remodelling during chronic pain.
How could engagement of these biological and cellular pathways lead to the patterned MS changes we report here? Our imaging transcriptomics findings for PLS11 suggest that neuroinflammation, with neuronal loss, synapse removal, and glial proliferation, 1 might drive loss of cortex differentiation during chronic pain. 36,44 This hypothesis matches well with recent preclinical models highlighting the role of neuroinflammation for neuronal sensitization of pain pathways, at both spinal and brain levels. 29,37,38 In humans, the neuroinflammation hypothesis has received direct support from positron emission tomography studies showing increased binding of ligands for the 18-kDa translocator protein, currently used as a marker of neuroinflammation and glial activation in the brain, [2][3][4]54,87 in the spinal cord, and nerve roots 2 of patients with chronic pain as compared with healthy controls. On the other hand, our transcriptomic-imaging association findings for PLS12 identify synaptic plasticity of excitatory and inhibitory neuronal populations as a potential driver of the increase in the MS of regions with constitutive low differentiation we report here. 44 Maladaptive neuronal plasticity, with dysfunctional regulation of the cortical E/I balance has been suggested to contribute to chronic pain. 44,60,70,72 Shifting of the E/I towards hyperexcitability, 73,88 as a consequence of either enhanced excitation or reduced inhibition, 49,53 is believed to augment central Figure 6. Gene expression profiles related to cross-condition changes in regional morphometric similarity during chronic pain. (A) We summarized case-control changes in regional morphometric similarity (MS) across chronic pain conditions using principal component analysis (PCA). The first component (PC1) explained the large majority of variance 64.45% and was the only component with an eigenvalue . 1. The cortical distribution of the scores of PC1 is depicted in the upper part of panel A. In the lower part, we show the cortical distribution of PLS1 scores summarizing the regional weighted expression of genes associated with crosscondition changes in regional MS during chronic pain. In the scatter plot on the right, we depict a significant positive correlation between PLS1 gene expression weights and the PC1 scores summarizing case-control regional MS differences across conditions. In panel B, we present the top 3 genes positively and negatively associated with PCA1, ranked by the respective loading into PLS1. Loading (Z-score) refers to the weight of each gene in PLS1. Genes with positive weights are highly expressed in regions where MS increases in patients, whereas genes with negative weights are highly expressed in regions where MS decreases. In the right part of panel B, we provide cortical maps summarizing the regional distribution of the top genes with the highest (CMTM3) or lowest (FAM126B) weights in PLS1. pain processing. 44,60 One key idea around this model is that peripheral injury triggers plastic changes or LTP in the cortical synapses. 13 Long-term potentiation promotes the formation of synapses and remodelling of dendritic spine substructures, 98 which compartmentalize calcium. 42,98 Hence, it is possible that complexification of the synaptic structure, reflecting long-lasting plastic changes in synaptic plasticity, might manifest as decreased MS. However, we should also acknowledge that neuroinflammatory and synaptic processes in the central nervous system tend to complement each other. 63,100 Therefore, assuming that increases and decreases in MS during chronic pain might involve mutually exclusive biological processes is likely simplistic. Future studies combining ex vivo MRI and histological examinations of the postmortem human brain of patients with chronic pain would be helpful in testing this transcriptomic regional vulnerability model further.
Our study has some limitations worth noting. First, our casecontrol data sets are relatively small, which might have affected our statistical power to detect small differences (particularly in the context of stringent correction for the number of regions examined). Therefore, future studies with larger sample sizes would be important to assess the replicability of our findings. Second, in the absence of multimodal data, we calculated MS using only 5 parameters from the T1-weighted images of each participant. This might have reduced the precision in estimating MS. Future studies attempting to replicate our findings should ideally use multiple imaging modalities (ie, DTI). Third, although still a general limitation of the field and not of this specific work, the whole-brain gene expression data derive only from 6 postmortem adult brains (mean age 5 43 years) and include data in the right hemisphere from 2 donors, which led us to exclude MS changes in the right hemisphere for the transcriptomic association analyses. By using constitutive gene expression in a small cohort of 6 postmortem brains to infer associations with neuroimaging markers acquired in different cohorts, we are assuming that regional gene expression is a conserved canonical signature that generalizes well beyond the brain samples included in the AHBA. Although we focused our analyses on probes that were selected to maximize differential stability across donors, 6 postmortem brains are insufficient to make strong claims about the stability of gene expression across brains in humans. Third, we pooled data from 3 different cohorts of patients with different chronic pain syndromes that were collected using different protocols and setups. This aspect poses limitations for investigating condition-specific changes in MS, which are likely to exist and might be interesting to pursue. Moreover, even within the boundaries of a specific chronic pain syndrome, it is likely that different pathophysiological mechanisms are in play in different patients. 64,82,97 This within-group heterogeneity is an aspect we did not deal with in this study but that future studies should take into consideration. Fourth, the data sets have varied, limited clinical information available, making it difficult to assess the clinical significance of the MS phenotype. Moreover, patients in the different studies were assessed using different clinical tools and detailed characterizations of each cohort were not available (ie, whether different patients might have received different treatments before enrolment). This lack of detailed and comparable data raises the question of whether our patients' cohorts were matched for important clinical variables, such as pain duration or disability. Differences in these factors, if existent, might have contributed to accentuate differences in MS between different chronic pain conditions. Ideally, future studies should aim to recruit, assess, and test all patients under the same protocol to minimize methodological heterogeneity. Fifth, whether MS changes are permanent or might be reversible after treatment is a question that Figure 7. Gene set and cell-type enrichment analyses of the top genes associated with cross-condition changes in regional morphometric similarity during chronic pain. In panel A, we present the results of gene set enrichment analyses on the top genes positively (PLS11) and negatively (PLS12) associated with crosscondition changes in regional morphometric similarity during chronic pain, as implemented in the Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA) platform. For PLS11, we present the top 10 Gene Ontology (GO)-biological pathways terms for which we found significant enrichment, ranked by P-value after FDR correction. For PLS12, none of the GO terms survived correction, but we found significant enrichment for 4 terms from the Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways. Colour scale indicates -log(P FDR ). The full results of the FUMA analysis can be found in supplementary data (available at http://links.lww.com/PAIN/B509, http://links.lww.com/PAIN/B510, http://links.lww.com/PAIN/B511). In panel B, we present the results of a cell-type enrichment analysis, where we investigated whether PLS11 and PLS12 include overrepresentation of genes typically expressed in specific brain cell types. Enrichment was quantified using odds ratio (OR) and significance calculated with a Fisher exact test. Colour scale indicates OR. Higher ORs (in red) indicate higher enrichment in genes of a certain cell class. The asterisk denotes cell classes for which we found significant enrichment, after correcting for the number of cell types tested (*P FDR , 0.05). FDR, false discovery rate; OPCs, oligodendrocyte precursor cells.
June 2022 · Volume 163 · Number 6 www.painjournalonline.com e769 we did not examine but that should be investigated in future longitudinal studies, given that previous studies have reported structural changes in the brain of patients with chronic pain that reverted after treatment. 30,74,79 Finally, although our findings are suggestive of a potential contribution of neuroimmune responses and neural plasticity to changes in MS, our correlational approach does not allow us to infer causality. This could be investigated further in longitudinal studies examining whether pharmacologic modulation of either pathway might attenuate the changes in MS we report here.
In summary, our study describes a new pattern of cortical MS remodelling across 3 chronic pain syndromes and identifies factors related to the glial immune response and imbalances in neuronal plasticity as candidates for molecular and cellular mechanisms conferring vulnerability to divergent tails of these cortical changes. Altogether, our data indicate that cortical MS remodelling in chronic pain entails a shared component of disease mechanisms that goes beyond specific clinical syndromes boundaries and might involve disruption of multiple elements of the cellular architecture of the brain which is unlikely to be efficiently targeted by current one-size-fits-all treatments. Ultimately, these findings highlight that developing new effective therapeutic approaches to the brain pathology that accompanies chronic pain might require a multitarget approach modulating both glial function and neuronal plasticity.