Brain white matter pathways of resilience to chronic back pain: a multisite validation

Chronic back pain (CBP) is a global health concern with significant societal and economic burden. While various predictors of back pain chronicity have been proposed, including demographic and psychosocial factors, neuroimaging studies have pointed to brain characteristics as predictors of CBP. However, large-scale, multisite validation of these predictors is currently lacking. In two independent longitudinal studies, we examined white matter diffusion imaging data and pain characteristics in patients with subacute back pain (SBP) over six- and 12-month periods. Diffusion data from individuals with CBP and healthy controls (HC) were analyzed for comparison. Whole-brain tract-based spatial statistics analyses revealed that a cluster in the right superior longitudinal fasciculus (SLF) tract had larger fractional anisotropy (FA) values in patients who recovered (SBPr) compared to those with persistent pain (SBPp), and predicted changes in pain severity. The SLF FA values accurately classified patients at baseline and follow-up in a third publicly available dataset (Area under the Receiver Operating Curve ~ 0.70). Notably, patients who recovered had FA values larger than those of HC suggesting a potential role of SLF integrity in resilience to CBP. Structural connectivity-based models also classified SBPp and SBPr patients from the three data sets (validation accuracy 67%). Our results validate the right SLF as a robust predictor of CBP development, with potential for clinical translation. Cognitive and behavioral processes dependent on the right SLF, such as proprioception and visuospatial attention, should be analyzed in subacute stages as they could prove important for back pain chronicity.


New Haven data set
Magnetic resonance imaging.Participants underwent an anatomical T1-weighted scan and two consecutive 2.5-minute-long diffusion tensor imaging (DTI) scans in the same session.A Siemens 3 T Trio B magnet (Erlangen, Germany) equipped with a 32-channel head coil was used to acquire the images.The 3D magnetization prepared rapid gradient echo (MPRAGE) T1-weighted acquisition sequence was as follows: repetition time/echo time (TR/TE) = 1900/2.52ms, flip angle = 9°, matrix size = 256 x 256, number of slices = 176, image resolution = = 1  1  1 mm 3. The diffusion-weighted images were acquired using spin-echo echo planar imaging (SE-EPI) sequence using the following scan parameters: TR/TE = 2200/84.0ms, flip angle = 90°, matrix size = 110  110  64, multiband acceleration factor = 4, image resolution = 2  2  2 mm 3 .Diffusion gradients were applied along 64 directions with a b-value of 1000 s/mm 2 .For each set of DTI data, one volume with no diffusion weighing (i.e., b=0 s/mm 2 ) was acquired at the beginning of the scan.

Mannheim data set
Participants.A trained psychologist interviewed all participants to assess comorbid mental disorders using the German version of the Structured Clinical Interviews (SCID I) for the Diagnostic and Statistical Manual of Mental Disorders (DSM IV)(1) Supplementary Table 2).In addition, all subjects completed the German versions of the Chronic Pain Grade (CPG), (2) the Örebro Musculoskeletal Pain Questionnaire (OMPQ yellow flag), (3) the Hospital Anxiety and Depression Scale (HADS),(4) and the Perceived Stress Scale (PSS).(5) Magnetic resonance imaging.A 3 Tesla Tim TRIO whole body scanner (SIEMENS Healthineers, Erlangen, Germany), equipped with a 12-channel head coil was used to acquire the images.Shimming of the scanner was done to account for maximum magnetic field homogeneity.Participants underwent an anatomical T1-weighted MPRAGE imaging scan with the following sequence: TR/TE = 2300/ 2.98 ms, flip angle = 9°, matrix size = 240  256, number of slices = 192, image resolution = 1  1  1 mm3.The diffusion weighted images were acquired using SE-EPI sequence with the scan parameters as follows: TR/TE = 7400/85 ms, matrix size = 220  256  30, GRAPPA = 2, image resolution = 2  2  2 mm 3 .Diffusion gradients were applied along 30 directions using a b-value of 1000 s/mm2.One volume with no diffusion weighing was acquired at the beginning of the scan.

Open Pain (Chicago) data set
Magnetic resonance imaging.Part of this data has been already published by Mansour et al.( 6) Therefore, we will briefly describe acquisition parameters.MPRAGE type T1anatomical brain images were acquired with a 3T Siemens Trio whole-body scanner with echo-planar imaging capability using the standard radio-frequency head coil with the following parameters: TR/TE = 2,500/3.36ms,flip angle = 9°, matrix size = 256 × 256; number of slices = 160, image resolution = 1 × 1 ×1 mm 3 .DTI images were acquired on the same day using SE-EPI with TR/TE =9000/83, flip angle = 90 degrees, in-plane matrix resolution = 112 × 130; 73 slices; image resolution = 2 × 2 × 2 mm 3 .Images had an isotropic distribution along 60 directions using a b value of 1000 s/mm 2 .For each set of diffusiontensor data, 8 volumes with no diffusion weighting were acquired at equidistant points throughout the acquisition.
Preprocessing of DTI Data.Preprocessing of all data sets was performed employing the same procedures and the FMRIB diffusion toolbox (FDT) running on FSL version 6.0.( 7) First, diffusion-weighted data were visually inspected for obvious artifacts or missing parts of the brain.Next, the data were corrected for eddy currents and head motion by employing affine registration to the no diffusion volume using eddy_openmp from the FSL toolbox.Eddy current corrects for image distortions due to susceptibility-induced distortions and eddy currents in the gradient coils.(8) We do note, however, that as we did not acquire data in the phase-opposite direction, the susceptibility-induced distortions may not be fully corrected.Brain images were then skull stripped and a diffusion tensor model was fit, using FMRIB Diffusion Toolbox (FDT) part of FSL (9), at each voxel to calculate the fractional anisotropy.(10,11) FA reflects the degree of water diffusion within a voxel with values ranging between 0 and 1 where large values indicate directional dependence of Brownian motion due to white matter tracts and small values indicate more isotropic diffusion and less directionality.(12) Harmonization of DTI data using neuroCombat.Because the 3 data sets originated from different sites using different MR data acquisition parameters and slightly different recruitment criteria, we applied neuroCombat (13) to correct for site effects and then repeated the TBSS analysis shown in Figure 1 and the validation analyses shown in Figures 5 and 6.First, the FA maps derived using FDT toolbox were pooled into one TBSS analysis where registration to a standard template FA template (FMRIB58_FA_1mm.nii.gzpart of FSL) was performed.Next, neuroCombat was applied to the FA maps as implemented in Python with batch (i.e., site) effect modeled with a vector containing 1 for New Haven, 2 for Chicago, and 3 for Mannheim originating maps respectively.The harmonized maps were then skeletonized to allow for TBSS.Questionnaire; PRSS, pain related self-statements; HADS, Hospital Anxiety and Depression Scale; HC, healthy control; CBP, chronic back pain; SBP, patients with subacute back pain; SBPp/r, patients with SBP with persistent pain or recovered pain after 6 months; FU, follow-up after 6 months; BL, baseline assessment; SD, standard deviation.Values show the group mean and standard deviation in parenthesis.a For the Chronic Pain Grade the median and interquartile range are shown.

Fig. S1 .
Fig. S1.Schematic representation of the sequential steps from left to right performed in structural connectivity data preparation, correction for confounders, and machine learning model building and testing.The dashed rectangle indicates that the combination of the data during model training was bootstrapped 50 times and validation and testing were repeated accordingly.

Fig. S2 .
Fig. S2.Head displacement (translation (A) and rotation (B) parameters) estimated by eddy current correction per each group in the Discovery data set.

Fig. S3 .
Fig. S3.Distribution plots of mean diffusivity (MD) extracted from the RSLF mask shown in Fig 1.The right SLF MD is significantly increased (p < 0.05) in the SBPr compared to SBPp patients in the New Haven data (A).The RSLF MD is not significantly different for the Mannheim data (p = 0.12) (B) nor for the Chicago data neither on visit 1 (p = 0.44) (C) nor on visit 4 (p = 0.38) (D).The MD values were compared between groups after accounting for age, sex, and motion parameters.The values on the plot are scaled up by 10 -4 for graphical depiction.

Fig. S4 .
Fig. S4.Distribution plots of radial diffusivity (RD) extracted from the RSLF mask shown in Fig 1.The right SLF RD is significantly decreased (p < 0.05) in the SBPr compared to SBPp patients in the New Haven data (A).The RSLF RD is not significantly different for the Mannheim data (p = 0.24) (B).The RSLF RD is not significantly different in the Chicago data neither on visit 1 (p = 0.78) (C) nor on visit 4 (p = 0.70) (D).The RD values were compared between groups after accounting for age, sex, and motion parameters.The values on the plot are scaled up by 10 -4 for graphical depiction.

Fig. S5 .
Fig. S5.Distribution plots of axial diffusivity (AD) extracted from the RSLF mask shown in Fig 1.The RSLF AD is not significantly different (p = 0.28) in the SBPr compared to SBPp patients in the New Haven data (A).The RSLF AD (p = 0.62) is not significantly different for the Mannheim data (B) nor for the Chicago data neither on visit 1 (p = 0.15) (C) nor on visit 4 (p = 0.36) (D).The AD values were compared between groups after accounting for age, sex, and motion parameters.The values on the plot are scaled up by 10 -3 for graphical depiction.

Fig. S6 .
Fig. S6.The receiver operating characteristic curve shows the area under the curve (AUC=0.53)when using the mask from the New Haven data to classify the SBPr and SBPp patients from Mannheim using the 30% recovery criterion.

Fig. S7 .
Fig. S7.Distribution of the values of the area under the curve (AUC) when shuffling the label of the groups (recovered, persistent) using the FA values from the right SLF before (A-C) and after (D-F) neuroCombat harmonization.The distributions for the Mannheim data are shown in A and D; the distributions of the Chicago data sets are shown in B and E and C and F respectively.The labels were shuffled 10,000 times; the p-values represent the probability of obtaining an AUC larger than the one obtained from the non-shuffled labels.

Fig. S8 .
Fig. S8.Correlation between FA values and pain severity with higher FA values (greater structural integrity) associated with the greater reduction in pain (percentage change) at baseline (A) and follow-up (B) in the Chicago data.

Fig. S9 .
Fig. S9.TBSS analysis and RSLF FA validation after applying neuroCombat after pooling subjects from the 3 sites into one analysis.(A) A whole brain comparison over the white matter skeleton between SBPr and SBPp patients at baseline and distribution of FA values for each group (New Haven data set).Results of unpaired t-test (p < 0.05, 10,000 permutations) showing significantly increased FA (fractional anisotropy) in SBPr (recovered) compared to SBPp (persistent) patients within the right superior longitudinal fasciculus (SLF) in the New Haven data set.(B-D) ROC shows the AUC when using the mask in A to classify the SBPr and SBPp patients at each of the respective sites.*, p < 0.05; **, p < 0.01.

Fig. S10 .
Fig. S10.Plot of FA values before (black) and after (red) applying harmonization using neuroCombat calculated as mean skeletal FA (A) and right SLF FA (B).Each subject's (x-axis) FA value (y-axis) is plotted for the 3 data sets Chicago, New Haven, and Manheim from left to right.FA values significantly changed when we compared the mean skeletal FA or right SLF FA before to after neuroCombat (paired t-test, p < 10 - 7 ).

Fig
Fig. S11.(A) Three-dimensional illustration of the left (L) SLF tracts connecting the inferior parietal lobe with prefrontal operculum (left), inferior parietal lobe with middle frontal gyrus (middle), and superior parietal lobe and primary motor cortex (right).Upper two rows show SLF tracts connecting the two regions in 3-D with FA values along those tracts shown in blue-to-red.The first row shows two representative SBPr patients and the second row shows 2 representative SBPp patients.(B) Plot of average FA ± SEM along the tracts depicted in A for 15 SBPr and 12 SBPp patients.The histogram plot shows the average number of tracts for these same patient groups.Abbreviations: IPL, inferior parietal lobe; PFO, prefrontal operculum; MFG, middle frontal gyrus; SPL, superior parietal lobe; M1, primary motor cortex.

Fig. S12 .
Fig. S12.Scatterplot of the testing AUC vs. the model validation AUC for all the machine learning models applied to the structural connectivity data.A combination of 40% of the New Haven data and 100% of the Open Pain data was used for model training and validation Testing was performed on 60% of the New Haven data (A), and the Mannheim data (B) Each point is the result of averaging the performance of the models (i.e., AUC) across 50 iterations as the combination of the data in model training stage was bootstrapped 50 times.The green circles were the final models considered in the result sections with an initial validation AUC ≥ 0.75.

Table S1 . New Haven sample characteristics
McGillPain Questionnaire; PCS, Pain Catastrophizing Scale.Values show the group mean and standard deviation in parenthesis.

Table S7 . Summary of the type of data combinations expressed in % of subjects from each site used to build and test the brain connectivity-based machine learning models to classify recovered and persistent SBP patients.
The combination of data sets was bootstrapped 50 times and the training and testing was repeated accordingly.