How do morphological alterations caused by chronic pain distribute across the brain? A meta-analytic co-alteration study.

•In chronic pain, gray matter (GM) alterations are not distributed randomly across the brain.•The pattern of co-alterations resembles that of brain connectivity.•The alterations' distribution partly rely on the pathways of functional connectivity.•This method allows us to identify tendencies in the distribution of GM co-alteration related to chronic pain.

Introduction

Chronic pain is a severe and disabling condition (frequently associated with physical and psychological comorbidities) that negatively impacts on the quality of life (Dominick et al., 2012, Gatchel, 2004, van Hecke et al., 2013, Walker et al., 2014). The International Association for the Study of Pain (IASP) defines chronic pain as the “pain persisting over the healing phase of an injury” (Loeser and Treede, 2008). Typically, this condition is considered to be chronic when persisting or recurring repeatedly for > 3 to 6 months (Jacobsen and Mariano, 2001, Merskey and Bogduk, 1994).

Chronic pain has been found to be related considerably to functional and structural reorganization in the nervous system (Apkarian, 2011, Apkarian et al., 2013, Apkarian et al., 2011, Baliki et al., 2014, Baliki et al., 2011, Cauda et al., 2014a, Farmer et al., 2012, Hashmi et al., 2013, May, 2008); this may partly explain why patients continue to experience pain even after nociceptive inputs are no longer present.

A number of studies suggest that chronic pain is associated with regional gray matter (GM) alterations within the brain (Apkarian et al., 2004, Baliki et al., 2011, Geha et al., 2008, Maeda et al., 2013, Obermann et al., 2013, Riederer et al., 2012, Rodriguez-Raecke et al., 2009, Schmidt-Wilcke et al., 2005, Seminowicz et al., 2010, Seminowicz et al., 2011, Tu et al., 2010, Ung et al., 2014, Unrath et al., 2007, Wood et al., 2009, Yang et al., 2013), as well as with regional white matter (WM) abnormalities (Ceko et al., 2013, Chen et al., 2011, Ellingson et al., 2013, Farmer et al., 2015, Gerstner et al., 2011, Khan et al., 2014, Lieberman et al., 2014, Luchtmann et al., 2014, Mansour et al., 2013, Moayedi et al., 2012, Woodworth et al., 2015) and functional or structural connectivity changes within large-scale brain networks (Baliki and Apkarian, 2015, Baliki et al., 2011, Baliki et al., 2014, Cauda et al., 2014a). Regional GM alterations (i.e., decreased or increased morphometric values in GM volume, concentration, or density) affecting various brain regions have been repeatedly reported in patients with chronic pain and meta-analytic approaches have proved to be helpful in identifying reproducible trends among studies (Cauda et al., 2014a, Smallwood et al., 2013).

For instance, Smallwood et al. (2013) found that GM decreased morphometric values are not only present in regions traditionally considered to be involved in pain processing (Iannetti and Mouraux, 2010, Isnard et al., 2011, Mazzola et al., 2006, Mazzola et al., 2009, Melzack, 1999) – i.e., insula, anterior cingulate cortex (ACC), and thalamus – but also in areas that do not seem to be specifically involved in pain processing (Cauda et al., 2014b, Cauda et al., 2014c) – e.g., superior temporal gyrus (STG), inferior frontal gyrus (IFG), and superior frontal gyrus (SFG). Another recent meta-analytic study (Cauda et al., 2014a) supports these findings, as it identifies in patients with chronic pain the prefrontal cortex, the anterior insula, the cingulate cortex, the basal ganglia, the thalamus, the periaqueductal gray (PAG), the post and pre-central gyri and the inferior parietal lobule as common sites of GM alterations. Intriguingly, when considered under the perspective of functional brain networks, the spatial patterns of GM alterations appear to encompass a core set of networks (i.e., the salience, the attentional and the default mode networks), which is commonly targeted in all the chronic pain conditions (Cauda et al., 2014a).

The network analysis employed by the aforementioned studies is a statistical approach that has already been used successfully to study both the structural and the functional connectivity changes associated with various clinical conditions (Crossley et al., 2014, Deco and Kringelbach, 2014, Filippi et al., 2013, Fornito and Bullmore, 2015, Fornito et al., 2015, Stam, 2014).

It is well-recognized that chronic pain can lead to structural brain alterations, which are mainly associated with regional GM volumetric decrease. GM alteration distribution is rather difficult to investigate, mainly because there is a lack of methods specifically devoted to the computation as well as the identification of a GM co-alteration pattern, that is, a type of connectivity that, instead of looking for the transmission of information, looks for the distribution of GM alterations. Recently, however, a study by Cauda et al. (2015) has proposed a methodology capable of identifying the neural alterations' distribution across the brain. Authors have shown that, at least for some psychiatric disorders, it is possible to define a morphometric co-alteration network and infer a hierarchy of brain structures within the GM alterations' patterns.

The present study follows the same line of research and aims at investigating the co-alteration pattern across the brain of patients with chronic pain by applying the method already developed by Cauda et al. (2015). In particular, we used a meta-analytic approach and a network-based analysis in order to address the following issues: i) do neuronal alterations occur randomly across the brain in patients with chronic pain?; ii) if the answer to the previous question is negative, can we identify a co-alteration pattern in terms of co-altered brain areas (nodes) and distribution links between them (edges)?; iii) does this co-alteration pattern overlaps with the pattern of functional connectivity formed by the same nodes?; iv) on the basis of the topological properties of the GM co-alteration network and the comparison between the co-alteration and functional networks, may the nodes with a high degree value play a pivotal role in the neuronal alterations' distribution?

Materials and methods

Search and selection of studies

We accepted the definition of meta-analysis of the Cochrane Collaboration (Green et al., 2008) and conducted a systematic search on PubMed database in order to identify all the voxel-based morphometry (VBM) studies reporting regional GM decreased or increased values in patients suffering from chronic pain. We used a search strategy based on the following combination of terms: “VBM” OR “Voxel Based Morphometry” AND “chronic pain” OR each disorder indicated by the American Chronic Pain Association (ACPA) as part of the spectrum of chronic pain (http://www.theacpa.org/7/Conditions.aspx). Afterwards, we examined the reference lists of the identified papers, searching for studies not found in the online database with our search strategy (for more information about the search strategy, see the Supplementary Material, Fig. S1).

The retrieved papers were further analysed in order to ascertain that they met the following inclusion criteria: (1) to be an original study; (2) to report GM changes (increased or decreased morphometric values) in patients with chronic pain by means of a between-group comparison with healthy controls; (3) to report the location of GM changes in Talairach/Tournoux or in Montreal Neurological Institute (MNI) coordinate system; (4) to use a well-specified VBM analysis; (5) to perform a whole-brain analysis (i.e., field of view not confined to a restricted region of cortex). We also checked that the duration of pain in all the patient groups was > of 3 months. We included the results of studies reporting both modulated and unmodulated VBM data. Modulation is a step in the VBM processing algorithm that adjusts for volume changes induced by normalization (i.e., the normalized GM segments are multiplied with the Jacobian determinant from the deformation matrix). GM values for Jacobian modulated images indicate regional volume changes. Conversely, if the data are unmodulated, GM values indicate regional brain density/concentration changes. As a large majority of the studies included in our meta-analysis report modulated data, we will henceforth refer to the GM volumetric changes simply as “GM increase” or “GM decrease”. In order to ensure the quality of data selection we followed the “PRISMA Statement” international guidelines (Liberati et al., 2009, Moher et al., 2009). For more information about the selection procedure, see the Supplementary Material, Fig. S1 (PRISMA Flow Diagram).

We used the BrainMap database (Fox and Lancaster, 2002, Fox et al., 2005, Laird et al., 2005c) in order to extract the coordinates from all the original studies in standard Talairach space. The VBM studies that were not already included on the BrainMap database at the time of the search (22 papers) were added so that all the original studies used in our meta-analysis (51 papers) are now available on BrainMap (see Supplementary Material, Fig. S1) and the reported coordinates of GM alteration can be extracted in standard stereotaxic space. In order to convert the coordinates from MNI to Talairach space, BrainMap uses the Lancaster transform (Lancaster et al., 2007) that ensures the quality of coordinates transformation.

On the basis of the aforementioned inclusion criteria, we identified 55 eligible VBM studies (51 papers) published before 31 January 2014, with a total of 2493 subjects: 1197 patients with chronic pain (mean age ± SD 45.70 ± 1.63) and 1296 healthy controls (mean age ± SD 44.35 ± 1.66). GM decrease was reported in 47 studies (1036 patients, 316 foci) while GM increase was reported in 27 studies (584 patients, 142 foci). For a complete list of the papers, see the Supplementary Material, Table S1.

Among the 55 studies included in our meta-analysis, 16 studies (11 with voxel-level and 5 with cluster-level correction) reported results with an initial correction for multiple comparisons, whereas 39 studies reported an initial whole-brain uncorrected threshold (generally, p < 0.001). In this second case, 12 studies reported a cluster-level correction (RFT, SVC, Monte Carlo simulation, FWE with p < 0.05); 5 studies reported a cluster extent threshold derived from a priori information about the size of the small anatomical structures involved in pain processing; 9 studies reported a minimum cluster extent threshold; and 13 studies reported no cluster-level correction or cluster extent threshold. Overall, 28 studies reported results corrected for multiple comparison and 27 studies reported uncorrected results. More detailed information about the statistical analysis performed by the studies included in our meta-analysis may be found in the Supplementary Materials, Table S4.

We used tables to collect the data regarding patients' characteristics, MRI acquisition, VBM algorithm and the number of reported foci of GM alteration. We also checked the information regarding the clinical condition, comorbidity, and medication of the experimental group, as well as the overlap in terms of age and sex between the group of patients and the group of controls. For a complete list of demographic, clinical and technical details of the studies included in our meta-analysis see the Supplementary Material, Tables S2 and S3.

Anatomical likelihood estimation analysis

To derive the nodes for the subsequent analysis we first employed the anatomical likelihood estimation (ALE) analysis (Laird et al., 2009, Laird et al., 2005b, Turkeltaub et al., 2002) to estimate consistent GM alterations across studies. Two separate analyses were carried out in order to identify the brain regions presenting values of GM decrease and increase, respectively. Coordinates in the Talairach space were used as inputs to perform the ALE analysis.

ALE is a quantitative voxel-based meta-analytical method that allows to estimate consistent areas of activation or morphological alteration across studies (Laird et al., 2005b). Bearing in mind that there is variability between studies due both to small sample sizes (between-subject variance) and to different normalization strategies (between-template variance), we applied an algorithm capable of estimating the spatial uncertainty for each focus by considering the possible differences related to sample amplitude and normalization procedures. This method allowed us to calculate the above chance clustering between experiments (random effect analysis) rather than between foci (fixed effect analysis) (Eickhoff et al., 2012, Eickhoff et al., 2009).

Modeled anatomical effect maps were computed for all the studies included in the meta-analysis. Each focus was modeled as the center of a 3D Gaussian probability distribution. We subdivided the Talairach space in 2 mm3 volumes and applied the following probability function (product of the three one-dimensional probability densities):

where d is the Euclidean distance between the voxels and the focus taken into account and σ is the standard deviation of the one-dimensional distribution.

The standard deviation σ was obtained through the full width at half maximum (FWHM):

In meta-analytic studies, the FWHM is used to model the spatial uncertainty of the reported foci (coordinates) of activation or morphological alteration. We used the quantitative uncertainty model proposed by Eickhoff et al. (2009) to calculate the FWHM for the Gaussian probability distributions:

This model takes into account the random effect and provides a quantitative estimate of between-subjects and between-template variability.

By adding these Gaussian functions, we obtained a statistical map where the likelihood of morphometric alteration is estimated for each voxel as determined by the whole set of studies. This map was thresholded by using a permutation test (Laird et al., 2005a, Lancaster et al., 2007, Lancaster et al., 2000). We are aware of the recent indications about correction for multiple comparison in ALE calculations (Eickhoff et al., 2016, Eickhoff et al., 2017) but since at this step our analysis used the ALE only to derive the peaks with the 25% strongest brain alterations, to avoid an excessive false negatives rate, we employed a more liberal correction for multiple comparison: the false discovery rate (FDR) with a q < 0.05 and a minimum cluster size of k > 100 mm3.

The final output of the ALE analysis is a thresholded statistical map of morphological alteration where each voxel has an associated value comprised between 0 and 1 indicating the likelihood that the voxel presents anatomical alteration.

Co-alteration network analysis

Separate network analyses were conducted for both GM increase and GM decrease.

Nodes creation

A customized MATLAB® routine was developed to individuate the peak coordinates of GM alterations and to create spherical ROIs centered on these coordinates. The node creation is obtained from the ALE map using a peak detection algorithm that returns the set of local maxima. A local peak is supposed to be a voxel whose ALE value is larger than the its neighboring voxels. Subsequently we selected the voxels with a peak value higher than a given threshold, which took into account the 75 percentiles of the peak values distribution (i.e., only the 25% of peaks exhibiting the highest ALE values were selected). Then we calculated the distance between each peak and obtained a distance matrix for each peak, in which all peaks within a distance of 10 mm from other peaks were excluded to avoid overlaps between ROIs. Around each of these peaks we designed a spherical ROI (diameter of 10 mm), which was used for the subsequent analysis.

Two-mode matrix

The previously selected ROIs (or nodes) were used to build a two-mode matrix. In the two-mode matrix the rows are constituted by the experiments/studies and the columns by the nodes. Individual entries report whether or not a node is altered in a given study. In case the node is altered, the matrix value is set to 1, otherwise to 0.

We created a co-alteration matrix using the previously defined set of nodes with a radius of 10 mm as template. In a N × M matrix each row represents an experiment, while each column represents a node. As in Crossley et al. (2013), for each experiment we considered a node as altered if 20% or more of the modeled alteration map (MA), thresholded at 0.001, overlapped with the node defined in the template.

Co-alteration matrix

The co-alteration matrix was calculated from the two-mode matrix using the Jaccard Index (Jaccard, 1901) as a measure of co-alteration. For each pair of nodes, the Jaccard Index was calculated as the number of studies reporting alteration in both nodes (intersection: A ∩ B) divided by the number of studies reporting alteration in either one of them (union: A ∪ B):

We thus obtained a co-alteration matrix, which is a square matrix (the nodes in the rows are the same as the nodes in the columns) with individual entries (jab) reporting the corresponding Jaccard co-occurrence coefficients for each pair of nodes.

The Jaccard matrix was then probabilistically thresholded. The statistical correction for non-significant connections was based on the method described by Toro et al. (2008). The co-alteration between pairs of nodes is defined by their statistical dependence across studies. We examined whether the morphometric alterations located in two brain areas (nodes X and Y) are best modeled as independent or as dependent events. For our purpose, the null hypothesis H0 states that the two regions (nodes X and Y) are altered independently form each other, while the alternative hypothesis H1 states that there is a co-alteration relationship between them (i.e., nodes X and Y are altered dependently of each other).

The likelihood of the null hypothesis is defined as:

where B refers to a binomial distribution in which n is the number of contrast that activated the second region, m the number of contrast that activated the first region, N the total number of contrast, p = m/N and k the numbers of contrast that activated both regions.

The likelihood of the alternative hypothesis is defined as:

with

and

The likelihood-ratio test (λ = L(H) / L(H)) was used to evaluate the likelihood of the alternative hypothesis H1 with respect to the null hypothesis H0. The distribution was shaped by a function with one freedom degree. Only connections that resulted significant at p < 0.01 (corrected for FDR) were maintained.

Subsequently, the Jackknife analysis was conducted to evaluate the robustness of our findings. The Jackknife is a non-parametric method that allows to estimate the sample distribution of a statistic (Fan and Wang, 1996, Radua and Mataix-Cols, 2009, Radua et al., 2011, Shao and Tu, 1995, Wu, 1986). Given a sample dataset of N elements, the desired statistic is calculated by systematically leaving out an element from the sample dataset. The operation is repeated for N separate samples of size N–1, each sample representing the original sample dataset with an element not included.

The Jackknife analysis provided us the connections (co-occurrences) that remain statistically significant when studies are excluded one at a time, allowing us to rule out the possibility that some connections could be derived from a small subset of studies. The co-alteration matrix so obtained reports co-occurrences, indicating that there are brain areas tending to be altered together.

Filtering the co-alteration matrix with functional connectivity data

Since several studies (Raj et al., 2012, Ravits, 2014, Seeley et al., 2009, Zhou et al., 2012) suggested that neuronal alterations may propagate within the brain along the pathways of brain connectivity, to minimize false positives we further analysed the co-alteration matrix so as to isolate the pathways that are also expression of a verifiable functional connectivity profile. We adopted this filtering strategy because our method for calculating the co-alteration pattern is new and we preferred to be conservative in order to ensure that all the possibly identifiable “co-alteration distribution pathways” actually represent ascertainable functional connections, which can be confirmed by functional connectome data.

Thus, in order to filter the co-alteration matrix, we use a set of resting state data from 100 subjects. From the same set of nodes obtained from the co-alteration matrix we calculate a functional connectivity matrix using rs-fMRI data from 100 healthy subjects stored in the Human Connectome Project (HCP) database (2015 Q4, 900-subject release). We required and obtained the permission from the HCP to work with the unrestricted HCP data (released to Dr. Tommaso Costa on 28 August 2013). As from the use term of the HCP connectome data, we also obtained the authorization from the University of Turin ethical board (24 September 2017).

The data are minimally preprocessed and ICA-FIX denoised and the age range is 22–35. The matrix is obtained with the following steps: we create a mask of spherical ROIs of 10 mm of diameter from the set of nodes of the co-alteration matrix; the mask is used in a dual regression that estimates for each subject an individual map of the ROI mask. More specifically: we regressed the ROI mask on each subject dataset to obtain a set of time courses as well as on the same dataset to obtain a subject-specific set of spatial maps.

The final result is a set of 100 matrices, one for each subject, where each column corresponds to a time series associated to a specific ROI. From these matrices, we calculated the partial correlation (L-2 norm Ridge Regression with sigma = 0.01) from each node and mediated from all the partial correlation matrices to obtain a final group correlation matrix. The resulting group functional connectivity matrix was thresholded at p < 0.01(corrected for multiple correction) using a one sample permutation test (5000 permutations) using the FSL randomize program. The survived functional connections among the nodes were used to constrain further analysis: in fact the resulting connections between the nodes that failed to have functional features were eliminated from the unfiltered co-alteration matrix, so as to ascertain that all the remaining connections could effectively correspond to cerebral connectivity pathways liable to be involved in the distribution of neuronal alterations.

Correlation between the co-alteration matrix and the functional connectivity matrix

In order to find out how much of the pattern of GM alterations' distribution tends to overlap with the pattern of cerebral functional connectivity, we used the Mantel test to calculated the correlation (Spearman's rho – ρ) between the co-alteration matrix and the functional connectivity matrix. The Mantel test (Mantel, 1967) is generally used to calculate the correlation between two distance or similarity matrices and its significance is evaluated through permutation procedures (e.g., Monte Carlo test). Because the elements in a distance matrix are not independent of each other, we cannot simply calculate the correlation coefficient between the two sets of n(n–1)/2 distances and test its statistical significance. To overcome this problem, a permutation test must be used. In order to assess the significance of any apparent departure from a zero correlation, the rows and columns of one of the two matrices are randomly permuted 5000 times, and after each permutation the correlation is recalculated. The significance of the observed correlation is the proportion of the permutations that leads to a higher correlation coefficient. The null hypothesis (H0) states that there is no relation between the two matrices. If H0 is true, then permuting the rows and columns of the matrix should be tantamount to produce a larger or a smaller coefficient. We therefore used this test to calculate the correlation between the two matrices. P value was obtained with 5000 permutations.

The co-alteration network and its topological analysis

The resulting filtered co-alteration matrix was used to construct the co-alteration network and run further topological analyses. We used the graph-theoretical analysis to investigate the regional topological properties of the co-alteration network and look for the degree centrality as the main network measure. The degree of a node is defined as the number of edges (connections) directly linked to it (Fig. 1).

Fig. 1

Workflow pipeline.

Results

The distribution pattern of GM co-alterations

Very interestingly our results indicate that the morphological alterations related to chronic pain do not distribute randomly within the brain. The probability of a non-random pattern of alteration distribution was high for both GM decrease and GM increase.

The highest degree nodes of the co-alteration network were the anterior and posterior insulae, the anterior cingulate cortex (ACC) and posterior cingulate cortex (PCC), the secondary somatosensory cortex (S2), and the pars opercularis of the inferior frontal gyrus (BA 44). The insula and the cingulate cortex were identified by previous studies as high degree nodes in functional (rs-fMRI) and structural (DTI) normal brain networks (Crossley et al., 2014, Crossley et al., 2013). Our functional connectivity analysis provided further evidence that these nodes are key regions in both the normal functional connectivity networks and in the co-alteration network (see Fig. S2 of the Supplementary material).

ALE analysis

Significant clusters of GM alterations (to be used for the subsequent ROIs creation step) are located in several cortical and subcortical areas, including anterior and posterior cingulate cortex, anterior and posterior insula, anterior prefrontal cortex (BA 10), dorsolateral prefrontal cortex (DLPFC), thalamus, putamen, BA 44, S2, temporal cortex, premotor cortex (BA 6). The results of ALE analysis are illustrated in Fig. 2.

Fig. 2

GM anatomical likelihood estimation results.

MRI alterations (i.e., gray matter increase/decrease) identified meta-analytically in chronic pain patients. The illustration summarizes the results of the anatomical likelihood estimation (ALE) analysis of all the papers involved in this study. ALE maps were computed at a FDR corrected threshold of p < 0.05, with a minimum cluster size of k > 100 mm3. Colors from red to yellow show gray matter increases, colors from blue to green show gray matter decreases. Images are shown using the right-left radiologic convention and standard Talairach space. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Nodes of GM alterations

Using the peak detection algorithm previously described, we identified 44 nodes of GM decrease and 23 nodes of GM increase. The results are illustrated in Fig. 3 and in the Supplementary Material (Tables S6 and S7). The Supplementary tables S6 and S7 summarize the information regarding the nodes' coordinates in Talairach space and the corresponding anatomical labels.

Fig. 3

Spherical ROIs (nodes) representing regions of GM decrease (a) and GM increase (b).

Comparison between GM co-alteration and functional connectivity matrices

Our assumption that the distribution of the neuronal alterations within the brain rely on the pathways of functional connectivity has been verified: indeed the chronic pain-related morphological co-alterations distribute according to the functional connectivity patterns. In fact, we found strong (and very significant) correlations between the co-alteration matrix and the functional connectivity matrix (Fig. 4) for both GM increase (ρ = 0.69, p < 3.6557e-07) and GM decrease (ρ = 0.23, p < 2.9667e− 06), thus confirming that a great part of the GM co-alteration pattern is overlapped to the pattern of resting state brain functional connectivity.

Fig. 4

Functional connectivity and co-alteration networks for a) GM decrease and b) GM increase.

The co-alteration matrices for GM decrease and GM increase were calculated using the Jaccard index: 42 out of the initial 44 nodes (95.45%) resulted connected in the GM decrease co-alteration matrix (Table 1), while only 12 out of the initial 23 nodes (52.17%) resulted connected in the GM increase co-alteration matrix (Table 2).

Table 1

GM decrease (D) nodes that resulted connected within the “co-alteration matrix”.

Node ID	Talairach coordinates			Label	Degree
	x	y	z
D1	2.0	− 24.0	− 18.0	r Pons	2
D2	58.0	− 40.0	− 14.0	r BA 20	4
D3	− 58.0	− 28.0	− 8.0	l BA 21	4
D4	28.0	20.0	− 8.0	r BA 47	3
D5	20.0	8.0	− 6.0	r Put	3
D6	− 32.0	12.0	− 6.0	l ant INS	8
D8	− 36.0	34.0	− 4.0	l BA 47	1
D9	− 40.0	− 8.0	− 2.0	l ant INS	8
D10	− 22.0	54.0	− 2.0	l BA 10	2
D11	− 12.0	− 26.0	0.0	l Pulv	1
D12	− 22.0	8.0	0.0	l Put	3
D13	2.0	56.0	0.0	r BA 10 (m)	5
D14	16.0	58.0	0.0	r BA 10(m)	5
D15	56.0	− 10.0	4.0	r BA 22	3
D17	44.0	14.0	4.0	r ant INS	13
D18	− 30.0	48.0	4.0	l BA 10	2
D19	10.0	− 28.0	6.0	r Pulv	2
D20	− 10.0	− 18.0	6.0	l MD	2
D21	− 48.0	− 22.0	8.0	l BA 41	3
D22	40.0	− 2.0	8.0	r post INS	9
D23	− 10.0	− 28.0	10.0	l Pulv	1
D25	60.0	− 4.0	12.0	r S2(BA43)	10
D26	54.0	10.0	12.0	r BA 44	8
D27	− 10.0	− 20.0	16.0	l LDN	2
D28	− 6.0	48.0	16.0	l DLPFC (BA9)	7
D29	58.0	− 12.0	18.0	r S2 (BA43)	5
D30	6.0	58.0	20.0	r BA 10 (m)	3
D31	56.0	10.0	22.0	r BA 44	3
D32	− 4.0	36.0	22.0	l ACC(BA32)	7
D33	− 10.0	30.0	28.0	l ACC(BA32)	3
D34	2.0	30.0	30.0	r ACC (BA32)	11
D35	40.0	28.0	32.0	r DLPFC (BA9)	2
D36	− 6.0	38.0	34.0	l BA 6	5
D37	10.0	− 46.0	38.0	r PCC (BA31)	2
D38	− 10.0	46.0	40.0	l BA 8	5
D39	10.0	− 32.0	42.0	r PCC (BA31)	8
D40	10.0	− 12.0	52.0	r BA 6	6
D41	34.0	− 10.0	52.0	r BA 6	4
D42	16.0	− 40.0	60.0	r S1(BA3)	4
D43	− 18.0	0.0	62.0	l BA 6	6
D44	14.0	0.0	64.0	r BA 6	7

l: left; r: right; m: medial; BA: Brodmann Area; INS: insula; Put: Putamen; Pulv: pulvinar; MD: medial dorsal nucleus of the thalamus; S2: secondary somatosensory cortex; LDN: lateral dorsal nucleus of the thalamus; DLPFC: dorsolateral prefrontal cortex; ACC: anterior cingulate cortex; PCC: posterior cingulated cortex; S1: primary somatosensory cortex.

Table 2

GM increase (I) nodes that resulted connected within the “co-alteration matrix”.

Node ID	Talairach coordinates			Label	Degree
	x	y	z
I9	12.0	6.0	2.0	r lat GP	2
I10	8.0	− 26.0	4.0	r Pulv	2
I11	− 10.0	− 24.0	4.0	l Pulv	1
I12	42.0	− 74.0	6.0	r V +	2
I15	− 20.0	20.0	6.0	l Caud Body	2
I16	− 26.0	− 10.0	8.0	l Put	1
I17	28.0	− 14.0	12.0	r Put	1
I18	16.0	4.0	12.0	r Caud Body	1
I20	− 22.0	52.0	22.0	l BA 10	1
I21	− 46.0	4.0	32.0	l DLPFC	1
I22	− 38.0	− 30.0	52.0	l S1	2
I23	38.0	− 22.0	52.0	r M1	2

l: left; r: right; BA: Brodmann Area; lat: lateral; GP: globus pallidus; Pulv: pulvinar; V +: associative visual areas; Caud Body: body of the caudate nucleus; Put: putamen; DLPFC: dorsolateral prefrontal cortex; S1: primary somatosensory cortex; M1: primary motor cortex.

Characterization of the GM decrease co-alteration network

Fig. 6, Fig. 8 illustrate the co-alteration network overlaid on anatomical brain images. There is a network configuration with a core set of high degree areas and long-range intra- and inter-hemispheric connections. Our analysis revealed that within the co-alteration network there are nodes characterized by high topological values (i.e., areas having a degree that is much larger than the average) that may play a centrality position in the alterations' distribution. The nodes with the highest degree value were: the right anterior insula (13 connections), the right ACC (11 connections), the right S2 (10 connections), and the right posterior insula (9 connections). Also the right PCC, the right BA 44 and the left anterior insula resulted to have high degree values (8 connections).

Fig. 6

Gray matter co-alteration network. The nodes represent regions of gray matter decrease. The size and color of the nodes reflect their degree value. Small sizes correspond to low degrees.

Fig. 8

Gray matter co-alteration network for both GM decrease and GM increase. The size and color of the nodes reflect their degree values. Small sizes correspond to low degree values. r: right; l: left.

All the network's highly-connected nodes resulted to be connected to both ipsilateral and contralateral areas. For each of these nodes, the number of corresponding intra- e inter-hemispheric connections (i.e., supposed pathways for the alteration spread) are illustrated in Fig. 5; whereas the highest degree nodes with their nearest neighbors are shown in Fig. 10.

Fig. 5

Highest degree areas in the co-alteration network. For each highly connected node, the figure shows the total number of connections, the number of inter-hemispheric connections, and the number of intra-hemispheric connections.

Fig. 10

GM decrease. Graphical illustration of highly connected nodes with their intra- and inter-hemispheric connections and their first neighbors.

As can be seen in Fig. 10, the highest degree areas of the co-alteration network tend to be highly interconnected, forming a pattern that resembles the rich-club organization of the human connectome. Furthermore, these nodes clearly appear to be lateralized on the right side, located in the anterior and posterior insula, ACC, S2, BA 44, and PCC of the right hemisphere. Fig. 9 illustrates the hierarchical organization of the nodes according to their degree.

Fig. 9

Graphical illustration of nodes corresponding to GM decrease. Nodes with the same degree value are arranged at the same level. Highest degree nodes are at the top of the figure and lowest degree nodes are at the bottom. r: right; l: left.

Characterization of the GM increase co-alteration network

For GM increase, the network analysis revealed that there are only low degree nodes, being 2 the maximum number of connections per node. This pattern of organization suggests that GM increase may be a more local phenomenon (Fig. 7, Fig. 8). All the identified pathways, except one (rM1 – rV +), were inter-hemispheric, some of them connecting homologous brain sub-cortical structures (e.g., bilateral pulvinar, lateral globus pallidus, and putamen).

Fig. 7

Gray matter co-alteration network. The nodes represent regions of gray matter increase. The size and color of the nodes reflect their degree value. Small sizes correspond to low degrees.

Discussion

Our results indicate that: i) in chronic pain, GM alterations do not distribute randomly in the brain; ii) we can clearly identify a morphological co-alteration pattern in the brain of patients with chronic pain; iii) there is an important overlap between the co-alteration pattern and the pattern of brain functional connectivity; iv) within the co-alteration network it is possible to identify a core set of highly-connected nodes (i.e., areas exhibiting a degree that is much greater than the average), which are supposed to play a centrality position in the distribution of neuronal alterations.

The distribution of neuronal alterations in chronic pain

Our study provides evidence that in the chronic pain condition neuronal alterations distribute in a non-random fashion; on the contrary, they form patterns that strongly resemble those of brain connectivity. Although network-based pattern of disease propagation was previously shown for neurodegenerative brain disorders (Crossley et al., 2014, Raj et al., 2012, Seeley et al., 2009, Zhou et al., 2012), this is the first time that the same phenomenon has been suggested to be the case also for chronic pain. Intriguingly, we observed a substantial overlap of our co-alteration matrix with the pattern of cerebral functional connectivity. In fact, we found a significant correlation between the co-alteration matrix and the functional connectivity matrix for both GM increase (ρ = 0.69, p < 3.6557e− 07) and GM decrease (ρ = 0.23, p < 2.9667e− 06).

This result confirms that a great part of the GM co-alteration network is related to the pattern of cerebral connectivity and, consequently, that in chronic pain the distribution of neuronal alterations can rely on functional pathways, which is consistent with recent studies about the spread of neurodegenerative disorders across the brain (Crossley et al., 2014, Zhou et al., 2012). In fact, recent evidence concerning neurodegenerative disorders seems to support the hypothesis according to which the neuronal/synaptic toxicity spreading across the brain, as well as the associated activity-dependent dysregulation of local misfolded proteins and metabolic levels, strongly depend on structural, functional and metabolic cerebral connectivity patterns (Iturria-Medina and Evans, 2015). Similarly, the nodal stress mechanism may play an important role in the distribution of neuronal alterations in patients with chronic pain. Specifically, the higher degree (i.e., highly-connected) areas might be particularly exposed to high levels of stress, and consequently to structural alterations, due to their intense functional activity (Buckner et al., 2009, Crossley et al., 2014, Saxena and Caroni, 2011, Zhou et al., 2012).

In this study the nodes with the highest number of connections were located in higher-order associative areas, such as the insula, the cingulate cortex and the secondary somatosensory area. It has been proposed that, due to the convergence of heteromodal activity, brain hubs are the most active regions in “rest” (i.e., non-task) conditions and that their intense connectivity-dependent neuronal activity may significantly diffuse neurodegeneration across the brain (de Haan et al., 2012, Iturria-Medina and Evans, 2015). This could explain the high vulnerability of hub areas and their pivotal role in the neuronal alterations' distribution.

Interestingly, some high-degree nodes of the co-alteration network seem to be preferentially targeted by chronic pain. This finding has important implications for the persistence as well as development of the neuronal alterations' pattern across the brain. In fact, when a highly connected functional node is affected by the pathological condition, the effect on the whole network functioning will be far more serious than when a little connected functional node is affected. This is supposed to be so because the altered highly connected node is an important station for information exchange between a large number of functionally connected areas, so that its disruption has a great impact on the whole network (and specifically on the regions directly linked to the altered highly connected node) and, thereby, facilitates the distribution of morphological alterations.

The neurobiological basis of the GM morphological alterations detectable with structural MRI (i.e. VBM) in patients with chronic pain remains as yet unclear. These alterations may result from irreversible mechanisms (i.e. neuronal degeneration/apoptosis), or from fast-adjusting, reversible processes, such as dendrite spine and synapse turnover (Agostini et al., 2013, Apkarian et al., 2004, Trachtenberg et al., 2002).

It has been suggested that GM volumetric increases related to chronic pain may reflect some kind of supraspinal inflammatory mechanism (Schweinhardt et al., 2008, Pomares et al., 2017), or result from an increased activity of certain neuronal populations with consequent synaptic level changes, which resembles what occurs in learning-related neuroplasticity (Draganski et al., 2011, Pomares et al., 2017). On the other hand, regional GM volumetric decreases related to chronic pain are commonly interpreted in terms of atrophy/neurodegeneration, though this interpretation is not yet fully accepted. Recent evidence suggests that GM decreases can be partially restored with successful therapeutic interventions, thus leading to the alleviation of pain (Baliki and Apkarian, 2015, Ceko et al., 2015, Gwilym et al., 2010, Rodriguez-Raecke et al., 2009, Rodriguez-Raecke et al., 2013, Seminowicz et al., 2011). Erpelding et al. (2016) documented rapid treatment-induced GM changes and functional connectivity changes within brain areas involved in sensation, emotion, cognition and pain modulation in pediatric patients with complex regional pain syndrome. However, there are few chronic pain conditions that may be successfully relieved with effective interventions, as was demonstrated, for example, for hip osteoarthritis (Rodriguez-Raecke et al., 2013). In particular, pain relief after hip joint endoprosthetic surgery is accompanied by GM volume normalization, suggesting that GM changes in chronic pain are not indicative of irreversible damage or atrophy. Moreover, the study of Rodriguez-Raecke et al. (2013) provides evidence that nociceptive inputs and disease-related motor impairments may induce processing changes in brain regions, with structural consequences which seem to be reversible. So, the successful treatment of pain relief may depend on the local variations in synaptic density, while the persistence of chronic pain may be associated with atrophic processes (Baliki et al., 2011, Baliki and Apkarian, 2015).

However, an important aspect about the possibility of GM normalization after pain relief is related to the timing of measurement. Long lasting GM changes need more time to be restored and, therefore, early measurements may not reveal significant GM normalization (Rodriguez-Raecke et al., 2013). Furthermore, other factors (e.g. elevated levels of cortisol and cytokines related to stress, pharmacological therapy) may induce additional functional and structural brain changes that hamper and/or prevent the reversal of GM decrease (Erpelding et al., 2016). It seems that pain duration and maladaptive processes may underpin treatment resistance, influencing the possibility of reversal of pain-related structural alterations. These observations underlie the importance of early interventions to maximise the possibility of treating effectively chronic pain.

In general, the cellular/microscopic-level changes that could explain the MRI-observed alterations (i.e. regional GM increase/decrease) include: increase/decrease in neuronal or glial cell size, neuronal or glial cell genesis/degeneration/apoptosis, changes in spine size and density, dendritic atrophy, angiogenesis and endothelial cell proliferation, and changes in the blood flow or in the interstitial fluid (Agostini et al., 2013, Apkarian et al., 2004, Keifer Jr. et al., 2015, Schmidt-Wilcke et al., 2005). Reversal of GM alterations due to pharmacological and/or non-pharmacological treatments may depend on the restoration of one or more of these neurobiological mechanisms.

Finally, it is noteworthy to discuss a methodological caveat concerning the relationship between our co-alteration network analysis and anatomical covariance (Evans, 2013, Mechelli et al., 2005). Anatomical covariations can be considered as “the covariance of morphological metrics derived from morphological MRI” (Evans, 2013). So it seems that the morphological co-alterations examined in this study may be conceived of as a type of anatomical covariance. However, methodologically speaking, the two approaches are utterly different, as anatomical covariance is always derived from single-subject data, while our meta-analytic method works on raw data obtained from a statistical comparison between pathological and healthy subjects. This is the reason why we prefer to avoid the expression “anatomical covariance” and, instead, call our approach co-alteration network analysis.

The GM decrease co-alteration network

Initially, most of the identified GM decrease co-alteration network nodes (42/44; 95.45%) were found to be interconnected in the GM co-alteration matrix, thus indicating that they tend to be altered together as elements of a common network structure. In order to ascertain how many nodes could be parts of functional networks as well, we decided to filter this matrix with functional connectivity data. The co-alteration network (composed of 41 nodes) exhibits long-range inter- and intra-hemispheric connections, as well as a number of highly interconnected brain areas. Within this network we were able to identify a core set of highly connected nodes, of which the most important are the right anterior insula, ACC, S2, posterior insula, PCC, BA 44, and left anterior insula. These areas form a core alteration pattern and, notably, have been equally identified as key hubs of the healthy structural and functional connectome (Crossley et al., 2014).

It is worth noting that in the co-alteration network the highly connected nodes appear to be significantly lateralized, with most of them located in the right hemisphere (Fig. 8, Fig. 9). The important role of the right hemisphere in pain processing has been highlighted by previous studies (Jensen et al., 2016, Ostrowsky et al., 2002, Symonds et al., 2006) and our findings provide further evidence to support this view. Another interesting finding is that the highest connected nodes of the co-alteration network closely resemble the rich-club organization of the human connectome (van den Heuvel and Sporns, 2011).

Recent studies that have investigated the pattern of structural brain alterations in neurodegenerative diseases by confronting the spatial distribution of the GM atrophy with the normal structural (DTI) or functional (rs-fMRI) connectivity profiles (Crossley et al., 2014, Zhou et al., 2012). These investigations have found that the regions showing higher structural/functional connectivity profile in the healthy brain networks as well as the regions showing shorter functional pathways to the network epicenters are the more vulnerable to pathological alterations. In agreement with these studies, our results have important clinical implications, as they suggest the possibility that the analysis of healthy structural and/or functional intrinsic connectivity patterns may help to track the distribution of pathological alterations within the brain and predict the regional vulnerability to disease. Since we found that in chronic pain the highest degree areas of the co-alteration network correspond to the highest degree areas of the normal functional connectivity network, it is reasonable to hypothesize that the nodal stress mechanism might play a central role in the development of GM neuronal alterations.

The GM increase co-alteration network

As to the GM increase data, only 12 out of the initial 23 nodes (52.17%) resulted connected in the GM co-alteration matrix, thus suggesting that GM increase may be a more localized and limited occurrence in patients with chronic pain. GM increase may also be a more variable phenomenon and, thereby, less easily detectable in VBM studies. The identified GM increase co-alteration network was characterized by a pattern with few and sparsely connected nodes. Notably, no node with more than two connections was identified.

The insular cortex

In our analysis the right anterior insula emerged as the region with the highest degree value. More specifically, the right anterior insula was found to be connected inter-hemispherically with the contralateral anterior insula, ACC and BA 6, and intra-hemispherically with ipsilateral ACC, PCC, S2, BA 44, putamen, BA 47, BA 44 and DLPFC. Interestingly, several regions connected to the anterior insula are also highly-connected nodes of the co-alteration network. These results are consistent with other studies that found the insula to elicit painful sensations when stimulated, as well as to be involved in pain processing, both in its anterior and posterior parts, which also exhibited a right-side lateralization (Jensen et al., 2016, Mazzola et al., 2009, Ostrowsky et al., 2002, Segerdahl et al., 2015).

Our connectivity model of neural co-alteration is also supported by animal studies. In fact, the antero-inferior division of the insula has been found to be strongly connected with the rostral anterior cingulate cortex (Vogt and Pandya, 1987, Vogt et al., 1987, Vogt, 1993, Vogt et al., 2004), and tract tracing studies in primates have shown that the insula is connected with the primary (S1) and secondary somatosensory cortex (S2), parietal operculum, prefrontal and motor cortex, supplementary motor area, orbitofrontal cortex, prefrontal and orbitofrontal areas, frontal operculum, anterior cingulate cortex, superior temporal gyrus, temporal pole, primary auditory and auditory association cortices, visual association cortex, olfactory bulb, amygdaloid body, hippocampus, and entorhinal cortex (Flynn et al., 1999, Mesulam and Mufson, 1982, Mufson and Mesulam, 1982). Most of these connections have been also confirmed by human studies of structural (DTI) and functional (rs-fMRI) connectivity (Cauda et al., 2011, Ghaziri et al., 2017).

Overall, our findings suggest that the insula is a preferentially targeted network node in various chronic pain conditions and appears to play an important role in the development of GM alterations across the brain.

The cingulate cortex, S2, and BA 44

The ACC node was located in the dorsal ACC (BA 32). This region is thought to be involved in a wide range of cognitive and emotional functions and has significant connections with the insula (Menon, 2015). As we have seen, both these areas are key regions of the salience network (Beissner et al., 2013, Menon, 2015, Menon and Uddin, 2010, Seeley et al., 2007), which supports the quick identification of external and/or internal stimuli that are relevant for the organism's survival (Menon, 2015). Chronic pain may disrupt the synchronization of this survival mechanism, affecting a common set of cerebral regions involved in pain processing, self-monitoring and salience detection (Cauda et al., 2010, Hemington et al., 2016, Kucyi and Davis, 2015, Legrain et al., 2011, Mouraux et al., 2011, Torta and Cauda, 2011).

S2 is involved in a series of important functions, including sensorimotor integration, attention, learning, memory, as well as the elaboration of nociceptive and pleasant inputs that are salient for further higher-order cognitive processing (Chen et al., 2008, Ferretti et al., 2003, Hamalainen et al., 2002, Ploner et al., 1999, Schnitzler and Ploner, 2000, Timmermann et al., 2001). S2 has significant connections with S1 and is supposed to play a relevant role in encoding the intensity of painful sensations (Lockwood et al., 2013).

PCC is densely connected with a high metabolic activity profile; its cognitive functions are principally devoted to internally directed cognition (Buckner et al., 2008, Hagmann et al., 2008, Leech and Sharp, 2014, Raichle et al., 2001). It is a key region of the default mode network, which is thought to be involved in retrieving autobiographical memories, planning the future and free wandering, as well as in integrating the sense of self-location and body ownership (Guterstam et al., 2015). Importantly, PCC seems to be involved in fine-tuning the meta-stability of intrinsic connectivity networks, thus allowing to control the focus of attention across time by regulating the variability of neuronal activation within the networks (Leech and Sharp, 2014).

The involvement of the right BA 44 in pain processing is less clear. In fact, its functions have been mainly associated with the elaboration of prosodic information (Hesling et al., 2005a, Hesling et al., 2005b, Wildgruber et al., 2005), generation of melodic phrases (Brown et al., 2006), and motor response inhibition (Bernal and Altman, 2009, Rubia et al., 2003). However, BA 44 seems to be also involved in action monitoring (Rubia et al., 2003) and lesions in this area have been associated with anosognosia for hemiplegia (Berti et al., 2005). Moreover, it has been suggested that the right frontal operculum (right IFG) may play a role in body ownership, in the resolution of conflicting signals between internal and external representations of body-related events (Tsakiris et al., 2007), and in the awareness of limb functioning (Kortte et al., 2015).

Relationship between the co-alteration matrix and diffusion matrix

At first glance the patterns of structural co-alterations and of alterations' spread seem to be unrelated. The former is typically derived from structural data and generally associated with a static event; the latter is instead studied by means of longitudinal investigations and frequently thought of as the result of causal events having a specific development in time. However, we propose to show that a deep mathematical relationship lies at the root of both patterns. In fact, the dynamics of how neuronal alterations spread across the brain can be described in mathematical terms with a Laplacian matrix, which, in turn, can be obtained from co-alteration data of meta-analytical origin.

If we consider two groups of neurons structurally connected (brain areas) – A1 that is altered and A2 that is not altered – then the alteration factor, spreading from A1 to A2, is the product of the concentration of the alteration factor x1 and the strength of the connectivity c12 between A1 and A2. As a result, at a certain time the concentration of the alterations in A2 will grow by a factor of βc12(x2 − x1)δt, in which β is the diffusion constant defining the pace of the alterations' spread. This system can be mathematically explained as follows:

Starting from this model Abdelnour et al. (2014) have showed the possibility to build a network of brain nodes by using the following equation:

in which L is the Laplacian matrix, expressed as:where D is the degree matrix, a diagonal matrix reporting the number of edges linked to every node, and A is the adjacency matrix, a square matrix reporting whether or not couples of nodes are adjacent or linked.

The Laplacian matrix can be proved to be tantamount to the heat equation, which, in turn, is a case of the diffusion equation generalized to complex networks (Kondor and Lafferty, 2002). This equation has the following solution:

The formula shows the development of a diffusion process, which begins from an initial stage x0.

Notably, to resolve the diffusion equation we need the Laplacian matrix, which can be derived from the co-alteration matrix of meta-analytical data. In fact, the co-alteration matrix describes the relations among the different altered nodes of a complex network (Crossley et al., 2013). This matrix can be constructed in several ways; to estimate the relationship between the altered nodes we used the Jaccard Index and to identify the significant relations between them we applied the same statistical method proposed in Toro et al. (2008). As a consequence, we can construct a square and symmetric matrix, which is an adjacency matrix A that reports in its rows and columns information about the co-occurrences of the altered network nodes. Thus, from this connection matrix we can obtain the degree matrix D as well as the Laplacian matrix L and in principle describe the diffusion of alterations within the complex network.

We can see now the mathematical relationship between the pattern of anatomical co-alterations and the pattern of alterations' spread in that the diffusion matrix of the alterations can be constructed from co-alteration data. In other words, the Laplacian matrix in Eq. (1) allows to describe the development of the alterations' spread. However, to do so we need to know the initial stage, that is, the starting configuration after which alterations begin to propagate (start condition x0). After knowing this onset we could describe how the alterations propagate by affecting sequentially the nodes. From data of meta-analytical source, however, it is not possible to recognize the point(s) of onset of the alteration process (i.e., the starting nodes). This is why with our methodology we cannot describe the exact progression of the pathological spread; in the article we therefore discuss the alterations' diffusion using the term “distribution” rather than the terms “spread” or “propagation”. Despite this limitation, we can maintain that, given the mathematical considerations stated above, alterations can develop from one node to another, thus allowing us the interpretation that the co-alteration pattern obtained from examining our meta-analytical data has a network-like architecture.

Limitations

Our findings indicate that the co-alteration network of chronic pain follows closely the brain patterns of functional connectivity. This intriguing result, however, could also be a limitation, as our co-alteration network has been compared with a model of healthy brain connectivity only based on one MRI modality (i.e., resting state functional connectivity). In fact, although scientific data suggest that the baseline functional pathways reflect to some extent the patterns of structural connections (Behrens and Sporns, 2012, Deco and Corbetta, 2011, Wang et al., 2013) the exact relationship between structural and functional connectivity is as yet unclear. Therefore further analyses on how structural and functional connectivity relates to each other are needed in order to understand whether or not our model is able to explain thoroughly the distribution of neuronal alterations in patients with chronic pain.

Another issue regards the high heterogeneity of chronic pain conditions, whose integration may be particularly challenging. Leaving aside the well-known aetiology-related issue, there are also differences in the methodological approaches used to study brain changes. Moreover, pain processing relies on distributed networks integrating various sensory, emotional and cognitive aspects, which, according to their different involvement, may lead to discrepancies between studies.

Therefore, differences in the brain location of GM changes between studies can be associated with a variety of factors, related to pain aetiology, somatotopic localization of pain, pain duration, pharmacological therapy, and co-morbidities (e.g. depression). For example, painful conditions with alteration of the nervous system (i.e. neuropathic pain) may differ from conditions in which the alteration of the nervous system does not occur. Differences can also regard the contribution of inflammatory, autoimmune or vascular components. However, with the help of a meta-analytic approach, common trends between various conditions can be identified. As previous meta-analyses suggest, chronic pain could lead to common structural and functional brain changes, independently of the possible condition-specific differences as well as of the aetiology or other condition-related factors (Cauda et al., 2014a, Smallwood et al., 2013).

In our case it was not possible to run separate analyses for different clinical conditions, as studies for each condition were not enough to allow appropriate statistical analysis. This is why we decided to include in our meta-analysis all the studies reporting GM alterations in chronic pain; in fact the inclusion of a substantial number of studies was necessary to achieve a significant statistical power.

Unfortunately, in the literature most of the studies about chronic pain (here included) do not have correction for multiple comparisons; instead they often use a more liberal thresholding or correction method like that applied by the study of Mordasini et al. (2012) relying on the study of Rüsch et al. (2003). This study suggests a more liberal thresholding in case of a clearly defined a priori hypothesis, based on previous knowledge regarding brain structures that are involved pain processing and modulation (e.g. midbrain, thalamus, putamen, amygdala, somatosensory cortex, insular cortex, cingulate cortex, prefrontal cortex, caudate nucleus, nucleus accumbens, amygdala). This thresholding strategy was mentioned in 16 studies, of which 13 used an initial uncorrected threshold for the whole brain analysis. Such a priori knowledge about the aforementioned brain structures was used to define an appropriate cluster extent threshold (5 studies) or to define a more liberal thresholding or correction method (11 studies), thus allowing to look for significant changes even in small structures that are known to be involved in pain processing an modulation (Table S4). Previous studies have shown that uncorrected levels larger than p < 0.001 can be included into the statistical parametric mapping analysis of gray matter volume differences if they match an priori hypothesis regarding the anatomical location of the findings (in our case, anatomical structures within the pain processing and modulation system) without losing protection against false-positive results (Mordasini et al., 2012, Rüsch et al., 2003).

To sum up, among the 55 studies included in our meta-analysis, 16 studies (11 with voxel-level and 5 with cluster-level correction) reported results with an initial correction for multiple comparisons, whereas 39 studies reported an initial whole-brain uncorrected threshold (generally, p < 0.001). In this second case, 12 studies reported a cluster-level correction (RFT, SVC, Monte Carlo simulation, FWE with p < 0.05); 5 studies reported a cluster extent threshold derived from a priori information about the size of the small anatomical structures involved in pain processing; 9 studies reported a minimum cluster extent threshold; and 13 studies reported no cluster-level correction or cluster extent threshold. Overall, 28 studies reported results corrected for multiple comparison and 27 studies reported uncorrected results. More detailed information about the statistical analysis performed by the studies included in our meta-analysis may be found in the Supplementary Materials, Table S4.

In the light of the previous considerations, we could not avoid the inclusion of studies reporting more liberal thresholded results, as the number of studies with fully corrected results was too limited for a significant statistical co-alteration analysis. So, although it is common practice in this field of research the use of uncorrected data for multiple comparisons (Dai et al., 2015, Pan et al., 2015, Smallwood et al., 2013, Yuan et al., 2017), we cannot rule out that this fact may have introduced a bias or inflated the statistics. To thoroughly address this issue, we hope that in the future, when a sufficient number of studies with corrected data will be present in the literature, it will be possible to provide a further meta-analysis that will be able to confirm our conclusions.

Finally, we cannot completely exclude that possible partial volume effects may have contaminated the results of the functional connectivity analysis between the nodes. In fact, the positioning of the nodes was driven by the anatomical data instead of the functional segmentation. Functional segmentation would have positioned the nodes in the areas of maximum functional homogeneity, thus minimising possible partial volume effects. In contrast, the anatomical positioning cannot minimize these effects from the beginning. However, the choice of the anatomical positioning was based on constraints related to the co-alteration analysis, specifically we wanted i) to position the nodes in the areas of maximum anatomical homogeneity and ii) to maximise the number of foci in each node.

Moreover, artefacts associated with partial volume effects are not likely to cause false positives in the correlation analysis between co-alteration and functional connectivity data. More probably, these artefacts could reduce the correlation value between the two types of connectivity. So, this consideration leads us to think that such artefacts, if present, did not alter significantly the similarity result between anatomical and functional matrices.

Conclusions

Our analysis provides further support to the hypothesis that neuronal alterations distribute according to the logic of a network-like propagation, not only in neurodegenerative disease but also in other pathological conditions such as chronic pain. We were able to answer all the four questions put forward in the introduction. In fact, we found out that i) in patients with chronic pain GM alterations do not distribute randomly, but ii) they rather form a symptom-related pattern of structural co-alterations, which iii) strongly mirrors the pattern of brain functional connectivity. Finally, iv) within this pattern or co-alteration network a set of highly connected nodes can be clearly identified. In particular, while areas showing GM increase seem to be a more local phenomenon, areas showing GM decrease tend to form a more densely connected network, with long-range inter- and intra-hemispheric connections. We can therefore hypothesize that these highly connected network nodes might play a pivotal role in the development and distribution of neuronal alterations in patients suffering from chronic pain; this suggests the nodal stress mechanism as the principle factor involved in this pathological process.

Another interesting finding is that the areas of the co-alteration network with the highest degree of connections are also parts of the salience network and appear to be more localized on the right-side, thus supporting the right lateralization in the processing of painful stimuli. This, in turn, suggests greater susceptibility of right-side highly connected nodes to alterations caused by chronic pain. These findings might have important implications for the emerging field of pathoconnectomics, paving the way for better strategies to track and predict symptom-related patterns of pathological alterations within the brain (Deco and Kringelbach, 2014, Filippi et al., 2013, Fornito and Bullmore, 2015, Fornito et al., 2015). What is more, they provide an additional argument in recommending the treatment of chronic pain as early as possible in order to prevent the development of GM pathological changes.