Structural Network Analysis Using Diffusion MRI Tractography in Parkinson's Disease and Correlations With Motor Impairment

Functional impairment of spatially distributed brain regions in Parkinson's disease (PD) suggests changes in integrative and segregative network characteristics, for which novel analysis methods are available. To assess underlying structural network differences between PD patients and controls, we employed MRI T1 gray matter segmentation and diffusion MRI tractography to construct connectivity matrices to compare patients and controls with data originating from two different centers. In the Dutch dataset (Data-NL), 14 PD patients, and 15 healthy controls were analyzed, while 19 patients and 18 controls were included in the Canadian dataset (Data-CA). All subjects underwent T1 and diffusion-weighted MRI. Patients were assessed with Part 3 of the Unified Parkinson's Disease Rating Scale (UPDRS). T1 images were segmented using FreeSurfer, while tractography was performed using ExploreDTI. The regions of interest from the FreeSurfer segmentation were combined with the white matter streamline sets resulting from the tractography, to construct connectivity matrices. From these matrices, both global and local efficiencies were calculated, which were compared between the PD and control groups and related to the UPDRS motor scores. The connectivity matrices showed consistent patterns among the four groups, without significant differences between PD patients and control subjects, either in Data-NL or in Data-CA. In Data-NL, however, global and local efficiencies correlated negatively with UPDRS scores at both the whole-brain and the nodal levels [false discovery rate (FDR) 0.05]. At the nodal level, particularly, the posterior parietal cortex showed a negative correlation between UPDRS and local efficiency, while global efficiency correlated negatively with the UPDRS in the sensorimotor cortex. The spatial patterns of negative correlations between UPDRS and parameters for network efficiency seen in Data-NL suggest subtle structural differences in PD that were below sensitivity thresholds in Data-CA. These correlations are in line with previously described functional differences. The methodological approaches to detect such differences are discussed.


INTRODUCTION
Neuronal degeneration in the substantia nigra resulting in dopamine deficiency in the basal ganglia is a major characteristic pathophysiological change in Parkinson's disease (PD). As the basal ganglia are involved in multiple corticosubcortical networks, PD can be viewed as an extended network disease of the brain (1). In the present study, we aimed to explore structural network differences in PD by combining gray matter (GM) segmentation, diffusion MRI tractography, and complex network analysis in two independent datasets.
The concept that PD symptoms and signs arise from functional impairment in coherent basal ganglia-cortical networks became generally acknowledged after acceptance of the model of segregated circuits (e.g., motor, oculomotor, and limbic), originating from the cortex via the basal ganglia and the thalamus back to particularly frontal cortical regions (2). The prominence of dopamine depletion in the posterior putamen (3) (which is part of the motor circuit in the model) is consistent with the prominent motor features in early PD stages. Further elaborations of the basic model included altered direct and indirect pathways within the basal ganglia (4), leading to an increased basal ganglia output to the thalamus and a subsequent decreased excitatory output back to the cortex. A physiological characteristic of PD-related changes in basal ganglia circuits is a more synchronous firing pattern (4), pointing at reduced segregation of neuronal activities of the basal ganglia loops in PD.
In vivo imaging of the human brain has complemented and extended insights gained from animal research and human histological examination of these cerebral pathways, both in health and in disease. Functional imaging with positron emission tomography (PET) using L-[ 18 F]fluorodopa has enabled the spatial identification and quantitative assessment of striatal dopamine deficiency, a key feature of PD (5), which is particularly pronounced in posterior parts of the putamen (6,7). Application of the PET tracer [ 18 F]fluorodeoxyglucose (FDG), a marker of regional cerebral metabolism, allowed for identification of a characteristic PD-related pattern of relative decreased metabolic activity in parietal, visual, and lateral premotor and prefrontal cortices and relative increases in the pons, thalamus, pallidum, dorsal putamen, primary motor cortex, and supplementary motor areas (8,9). Functional interactions between such regions can be identified with resting state fMRI and targeted activation paradigms. For example, resting state fMRI enabled the demonstration of reduced functional connectivity between the posterior putamen and inferior parietal cortex together with increased anterior putamen-inferior parietal coupling, a remapping considered to reflect reduced spatial segregation between different corticostriatal loops in PD (10). By applying an fMRI visual optic-flow paradigm, mimicking the perception of forward locomotion, we previously found that interruption of such a gait-supporting stimulus failed to activate the (pre-)supplementary motor area in PD, while functional connectivity between this region and the visual motion area V5 was enhanced in patients, a result which is consistent with the increased interference of perceptual stimuli with motor intentions in PD (11).
The functional connectivity and inferred interactions between spatially distributed brain regions in fMRI data are based on the temporal correlations of signal changes in such regions, either evoked in hypothesis-driven experiments, or spontaneously occurring (12,13). A basic network model of small-world architecture, describing the dynamic consequences of local and remote interconnectivity features (14) has strongly encouraged the development of whole-brain "complex network" analysis (15,16). Applications of the latter in neuroimaging data made it possible to reach higher levels of exploring the brain's neuronal organization that underlies the balance between functional segregation and integration. This coincides with regional modular integration in specialized brain regions as well as global integration based on remote whole-brain interconnections of regional modules. While "connectivity" in the obtained datasets may concern networks of either actual white matter connections or functional associations (13,15), these networks can generally be mathematically defined as a collection of nodes and edges (links) between pairs of nodes (15). This approach enables the construction of connectivity matrices and the subsequent calculation of various measures, such as global efficiency (E glob ) and local efficiency (E loc ) that represent network properties concerning integration and segregation, respectively (16,17).
In PD, the wide spectrum of specific motor (18) and nonmotor (19) symptoms provides a challenge to capture differences in the integrative and segregative properties of the PD-affected cerebral networks. Complex network approaches using resting state fMRI have shown decreased whole brain E glob in PD (20,21) as well as decreased motor network E glob in preselected cortical and basal ganglia regions (22). E loc differences in distinct nodes have also been reported without differences in average whole-brain E glob and E loc (23). In mildly cognitively impaired PD patients, Pereira et al. demonstrated reduced E glob compared to controls, deriving network parameters from correlation matrices based on (structural) cortical thickness measures (24). While correlation matrices on fMRI data and regional cortical thickness measurements use edges between nodes without taking into account the existence of actual white matter connections, complex network analysis on white matter tracts that interconnect distinct GM regions requires the reliable identification of such structural edges, e.g., based on diffusion MRI tractography. With diffusion MRI, E glob reduction has been shown in PD (25)(26)(27). The number of complex network studies that have been published on PDassociated differences in structural cerebral connectivity however remains limited, in such a way that it remains difficult to draw unequivocal conclusions concerning changes in specific networks (28).
In the present PD imaging study, we performed complex network analysis on structural white matter connectivity. This required the initial steps of (i) GM segmentation to provide regions of interest (ROIs) and (ii) diffusion MR tractography rendering streamlines. From these data, connectivity matrices were constructed serving the model of nodes and edges, to be used in subsequent calculations. This was performed in two independent datasets, each consisting of PD patients and control subjects. We hypothesized finding reductions of both E glob and E loc in PD patients and negative correlations of these measures with clinical motor scores in interconnected brain regions previously implicated in functionally impaired networks in PD.

MATERIALS AND METHODS
Two MRI datasets of PD patients and controls were acquired. One dataset [Dutch dataset ("Data-NL")] was acquired at the University Medical Center Groningen, the Netherlands. The other [Canadian dataset ("Data-CA")] was acquired at the Pacific Parkinson's Research Center of the University of British Columbia, Vancouver, Canada. Unless specified otherwise, the same procedures were followed with both datasets.

Data-NL
For the first dataset, 15 PD patients and 16 healthy controls (HCs) were included. During processing of the imaging data, two subjects had to be excluded from further analysis due to suboptimal quality of the data, leaving 14 patients and 15 agematched controls for the final analyses. PD patients fulfilled the UK Brain Bank criteria for PD (29,30). HCs were required not to have first-degree family members with parkinsonism or dementia. All subjects underwent MRI scanning. Before MRI, antiparkinson medication and benzodiazepines were withheld for at least 12 and 24 h, respectively. Part 3 of the Unified Parkinson's Disease Rating Scale (UPDRS) was applied in the PD group. Table 1 shows some demographics of both datasets.

Data-CA
In Data-CA, 19 subjects with PD and 18 age-matched HCs were enrolled. All PD patients were diagnosed and under the care of a tertiary-care movement disorders specialist. Exclusion criteria included atypical parkinsonism, other neurological or psychiatric conditions, and use of antidepressants, hypnotics, or dopamine blocking agents. In this dataset, part 3 of the UPDRS was applied in the PD group as well. Both motor assessment and image acquisition were done after withdrawal of L-DOPA for 12 h and dopamine agonists for 18 h.

Image Processing
First, both T1-weighted and diffusion-weighted images were processed to render the ROIs and the streamlines, which were subsequently combined to form the connectivity matrices. From these matrices, E glob and E loc were calculated. See Figure 1 for a graphical overview of the image processing pipeline.

GM Segmentation
For GM segmentation (providing nodes for the complex network analysis), the T1 images were processed using the FreeSurfer pipeline (31) (RRID:SCR_001847), which included defining ROIs according to the Desikan-Killiany atlas (32). FreeSurfer version 5.3 was applied on a cluster of processors running Linux (Ubuntu 12.04.5 LTS, CPU model AMD Opteron Processor 6272). Default parameters were used. The output of every subject was checked visually by viewing the subcortical segmentation and the white and pial surfaces overlaid on coronal, sagittal, and axial T1 slices. If necessary, adjustments were made to the output followed by rerunning parts of the FreeSurfer pipeline. As a result, data from all subjects were suitable for further analysis. The resulting files containing the ROIs (the aparc+aseg.mgz files) and the T1 volumes (T1.mgz files) were converted to the nifti file format and stored for this purpose. It should be noted that using these methods, the ROIs were in the individual subject space, reducing the need for non-linear spatial transformations of the diffusionweighted images.
In the ROI file, ROIs that were of no interest, for example, the ventricles and the cerebellum (which was not fully covered in all acquisitions), were excluded from the analysis, resulting in 85 ROIs for each subject. These ROIs comprised all cortical ROIs from the Desikan-Killiany atlas (33) as well as the thalamus proper, caudate, putamen, pallidum, hippocampus, amygdala, accumbens area and ventral DC (all bilateral), and brainstem.
The putamina were separated by the vertical plane traversing the anterior commissure in an anterior part and a posterior part, taking into account that the posterior part is more severely affected in PD (3,6,7). This was achieved by (a) applying the talairach.xfm transform as saved by the FreeSurfer output to the T1.mgz file (i.e., resulting in the anterior-posterior direction of the voxels in this new file being parallel to the anteriorposterior direction of the brain), (b) using this file to construct a plane through the anterior commissure and orthogonal to the anterior commissure-posterior commissure line with a thickness of 2 mm, (c) applying the inverse of the talairach.xfm to this plane, and (d) separating the putamina in native space according to this plane.

White Matter Processing
The white matter analyses were performed using software running under MATLAB, MathWorks, Inc., Natick, Massachusetts, USA. In order to define the streamlines (providing edges for the complex network analysis) for each subject, the diffusion-weighted images were corrected for motion, eddy currents, and echo planar imaging (EPI) distortions, followed by tractography (see pipeline below). The DWI data were processed in ExploreDTI (34) (RRID:SCR_001643) running under Linux (Ubuntu 10.04 LTS, CPU model Intel Core i7-2600). The pipeline consisted of the following steps: 1. To aid in later registration, the FreeSurfer T1 nifti files were masked using ExploreDTI (and checked visually afterwards), applying a kernel size of morphological operators of 5 and a threshold of 0.05. 2. For each subject in Data-CA, the three DWI sessions were concatenated resulting in one file with three b = 0 s/mm 2 volumes and 96 b = 700 s/mm 2 volumes. 3. Diffusion gradient directions were corrected for the image angulation as described in the Philips PAR header. 4. The DWI volumes of all subjects were corrected for motion, eddy currents, and EPI distortions using default parameters (33,(35)(36)(37). As a part of this, the (originally distorted) DWI volumes were aligned with the masked T1 file. 5. Tensor estimation was performed using the REKINDLE method (38). 6. The processed diffusion data were checked visually for each subject by (a) viewing all three planes of the corrected diffusion-weighted images in a movie loop, (b) viewing the axial slices of the color-coded fractional anisotropy (FA) map, and (c) viewing all three planes of the color-coded FA map overlaid on the T1 volume. Data-CA did not raise concerns in this regard, but Data-NL showed many erroneous motion correction results, due to the high b-value employed (4,000) and the accompanying low SNR. Therefore, the correction process for this dataset was rerun using adjusted motion correction settings. The number of histogram bins was set to 16 and the number of data samples to 8,000. Also, the "Scales" parameters for the Elastix program were adjusted so that rotation, scaling, and shearing were constrained (set to 10 6 , whereas the Scales parameters for translation were kept at 1). After rerunning all subjects' scans with these settings and using the quality assurance program as described above, all datasets were suitable for further analysis, except for two datasets which were excluded as the correction results remained poor. 7. Constrained spherical deconvolution (CSD) whole-brain tractography (39)(40)(41) was performed using the following parameters: step size 1 mm, fiber orientation distribution (FOD) threshold 0.1, angle threshold 30 • , length range 50-500 mm, and seed point resolution 2 × 2 × 2 mm. The maximum harmonic order was set to 8 for Data-NL and to 4 for Data-CA. This resulted in one set of streamlines for each subject.

Complex Network Analysis
Connectivity matrices were created for each subject by combining the ROIs (nodes) and the streamlines (edges). The matrices were constructed with 87 rows and 87 columns (for 87 ROIs) in which the intersection of two ROIs was labeled by the mean FA if at least one streamline was found between these ROIs. Zeros were placed both at intersections without a streamline and on the main diagonal. Similarly, connectivity matrices weighted by the number of streamlines (NOS), and binary connectivity matrices (i.e., the distinction between connections present yes or no) were created and used in subsequent analyses. For the network analyses, an initial threshold was applied by discarding all connections present in <50% of subjects (42). That is, for every possible connection out of a total of 87 × 86 = 7,482 connections, only those connections that were found using tractography in at least 50% of all subjects would be used in subsequent analyses. In this way, we aimed to reach a balance between the number of false-positive streamlines (also called spurious streamlines, i.e., following trajectories without the presence of actual underlying fibers) and false-negative streamlines (i.e., underlying trajectories not being found by the streamline algorithm). This was done independently for each dataset. Next, the FA weighted connectivity matrices were normalized per dataset; i.e., all FA values in these matrices were divided by the maximum value of all FA matrices within the corresponding dataset. The same was done with the NOS matrices. This resulted in values between 0 and 1, which is a prerequisite for the subsequent complex network analysis, while at the same time, relative differences between subjects remained. The E glob and E loc were calculated using the Brain Connectivity Toolbox that accompanied a 2010 paper by Rubinov and Sporns (16). E loc was calculated using the efficiency_bin (for binary matrices) and efficiency_wei (for weighted matrices) functions from the toolbox. E glob was calculated using the formula presented in Appendix A of the paper; custom code and the distance_wei function from the toolbox were used. For each connectivity matrix (of each subject), E glob and E loc were calculated across the entire brain, providing whole-brain values. Whole-brain E glob and E loc were constructed by taking the mean of the nodal values of E glob and E loc , respectively. Nodal values for E glob and E loc were also stored separately.

Statistical Analysis
Within each dataset and for the different weightings separately, whole-brain E glob and whole-brain E loc values were compared between PD patients and HCs using a Student t-test (alpha 0.05) and related to UPDRS motor scores within the PD group using Pearson's correlation. The same analyses were performed for nodal E glob and nodal E loc and then for each node separately. Here, a false discovery rate (FDR) of 0.05 was applied to correct for multiple comparisons, for each measure and weighting separately.

RESULTS
A qualitative assessment of the FA weighted connectivity matrices revealed a clear consistency between the patterns of the four groups (Figure 2). These matrices were constructed after thresholding and normalization (during which the Data-NL FA values were divided by 0.666 and the Data-CA values by 0.661, which were the respective maxima of the corresponding datasets) and averaging per group. The similarity is illustrated by the robust cross-callosal connections between left and right occipital regions (via the splenium) and between left and right (pre)motor areas in the two datasets. Furthermore, the mean FA per connection was distributed similarly in both datasets (higher pixel values show at similar locations). Also see Supplementary Figure 1 which can be enlarged, in such a way that the labels of all separate ROIs can be read, allowing for a more detailed assessment of the connectivity matrices.

Data-NL
Comparing FA weighted network parameters between the PD and HC groups in Data-NL [mean (std)], whole-brain E glob [PD 0.29 (0.02); HC 0.29 (0.02)], and whole-brain E loc [PD 0.40 (0.03); HC 0.40 (0.03)] did not yield significant differences. Neither were such differences found at regional levels when comparing the nodal E glob and nodal E loc values. However, we did see significant correlations between network measures and clinical PD parameters. Whole-brain E glob and E loc negatively correlated with the UPDRS motor score (E glob p = 0.02 and E loc p = 0.03), meaning that lower efficiency scores accompanied worse motor symptoms (see Figure 3). At the nodal level, significant negative correlations (using an FDR of 0.05) between E glob and E loc with the UPDRS motor score resulted in two distinct patterns of regional cortical involvement: particularly, the right superior parietal cortex and right posterior cingulate cortex showed strong negative correlations between the UPDRS motor score and E loc , while for E glob , the strongest negative correlation with the UPDRS motor score concerned the sensorimotor cortex (Figure 4). At more lenient FDR values (e.g., up to 0.2), the number of regions in these patterns considerably increased, with an overlap between the two patterns. Analysis of the NOS and binary matrices did not result in any significant differences or correlations.

DISCUSSION
The complex network analysis on connectivity matrices constructed using GM segmentation and DWI tractography did not result in significant differences when comparing the PD patients and HC, either in Data-NL or in Data-CA. This suggests that PD did not lead to obvious changes in anatomical interconnection. However, while the connectivity matrices of the four groups indeed showed similar patterns, the negative correlations between both E glob and E loc and the UPDRS motor score in PD of Data-NL (at whole-brain and nodal levels) point to subtle disease-related reductions of integrative and segregative properties of the brain's structural network in PD. To possibly explain why these correlations were identified only in Data-NL and not in Data-CA, clinical explanations and methodological issues will be addressed.

Clinical Neuroscience Considerations
No differences were found between patients and HCs, whereas decreases of efficiency values have previously been shown with fMRI (20)(21)(22) and structural MRI (24)(25)(26). At the whole-brain level, differences in E glob and E loc provide a general indication of disease-associated differences in the integration and segregation capacities of the entire cerebral network. At the nodal level, E glob and E loc differences may provide insight into the changed contributions of specific brain regions to functional integration and regional specialization. It should be considered that, while structural network differences may be subtle in PD, functional network differences identified with fMRI may be more robust, preceding structural changes. We did, however, find that E glob /E loc of structural interconnections correlated negatively with UPDRS scores in PD, which is generally consistent with a decrease of efficiency values in previous PD studies. The observed PD-related effect on posterior parietal nodal level E loc suggests impaired efficiency of posterior parietal processing in PD, which is consistent with the enhanced impact of visual cues on motor actions (43,44), including gait (45). It has been argued that reduced dorsal visual pathway function due to less elaborate early-stage visuomotor processing results in functional shortcuts to remote mediofrontal motor regions (11). The impact of PD pathology on posterior cortical regions can also be inferred from parietal atrophy (46) and the PD-related profile of cerebral metabolism which is characterized by reduced parietal metabolism (8,9). With structural connectome-based analysis, reduced E loc in particularly parietal regions were very recently identified in a large PD patient group (n = 65) compared to 65 matched HCs (47), thus providing support for the observed parietal correlation in our Data-NL patient group. Complementary to the PD-associated reduction of parietal E loc , one might speculate that the reduced E glob (at the nodal level) of the motor cortex is based on an increase of non-specific network input from remote regions, in addition to impairment of the specific input from basal ganglia/thalamus loops and premotor regions. Such a model would fit the previously described decrease of network efficiency centered onto the motor cortex (20,22,48) together with relatively increased metabolism at this motor target site (8,9).

Methodological Considerations
In Data-NL, we identified negative correlations between the above-described parameters for network efficiency and UPDRS motor scores, which was not found in Data-CA. An important difference in data acquisition between the two groups concerned the employed b-value, which was higher in Data-NL (4,000), leading to a lower SNR than in Data-CA (b-value 700). In Data-NL, this resulted in suboptimal tensor estimations and accompanying FA values which tended to be lower (illustrated in Figure 2). Still, CSD tractography benefits from this b-value together with 60 diffusion gradient directions, compared to the more usual b-value of around 1,000 and fewer diffusion gradient directions, especially when it comes to resolving fiber crossings (49). The b-value of 700 in Data-CA thus assured FA estimations to be more reliable; however, together with the limited number (32) of diffusion directions (and the accompanying maximum harmonic order of 4), it does not serve the best interests of CSD. For example, crossing fibers are better separable at a maximum harmonic order of 8. As a result, more connections were found in Data-NL reflected by more dense connectivity matrices (see Figure 2). Ideally, probably an intermediate bvalue (or multiple b-values) and at least 60 diffusion gradients should be employed. Nonetheless, the patterns in the connectivity matrices were very similar between datasets, suggesting that these matrices may be constructed robustly, despite the applied (very different) acquisition parameters, and the adjustments in the motion correction tuning parameters (i.e., restricting the rotation, sharing, and scaling) that had to be applied in Data-NL.
Assessing tracts between cerebral regions implies demarcation of such regions. In the present study, the latter was performed by defining ROIs according to the Desikan-Killiany atlas (32). It may be evident that improvement in defining nodes and edges is a topic of further development, although to date, no simple solutions exist. Whereas ideally nodes and edges would form a complete (i.e., connections at all locations of underlying biological connections and nowhere else), weighted (i.e., connectivity measures correlating with "connection strength"), and signed (i.e., reflecting the direction of information processing) connectivity matrix, we are far from quantifying the brain's network in this way when constructing diffusion MRI tractography-based connectivity matrices (50)(51)(52). Several ways of defining the nodes are available, for example, based on atlases, randomness, functional information, and voxel (50), each with its own strengths and weaknesses. Still, the presently applied Desikan-Killiany atlas has proven to be successful in the past, e.g., (53). Furthermore, we chose to use the FA and NOS weighted and binary (i.e., unweighted) connectivity matrices. Next to the problem of false-negative and false-positive streamlines, there are a number of other methodological considerations on why these may not be suitable for quantifying "connectivity health" (51). In general, we are looking at streamlines drawn in the diffusion field rather than tracking the actual fibers (albeit with the assumption that they somehow correlate) and at best infer indirect measures of the integrity of these underlying fibers. The impossibility of showing differences using these techniques does not mean that integration and segregation functions of the brain do not differ between groups. However, to date, the weightings we employed are the most widely used way of weighting the connectivity matrices. Also, in comparable analyses, intrasubject variability has been shown to be lower than intersubject variability, suggesting that biological processes are being measured (54).
A general issue to consider with regard to (i) the absence of differences when comparing PD and HC and (ii) the observed anticorrelation between network efficiency parameters in only Data-NL and not in Data-CA is the fact that PD is a heterogeneous disease, not only with regard to the spectrum of symptoms but also regarding disease progression (18,19,55). Together with the limited number of subjects, disease heterogeneity may have compromised sensitivity to find differences within and between the datasets. On the other hand, one might argue that PD patients and controls did indeed not differ with respect to the integrative and segregative properties of their brains' structural networks. In this respect, it should be kept in mind that functional changes may likely develop before the onset of structural network changes, which might contribute to the discrepancy between previous functional studies and our structural graph analyses. Including larger patient numbers and more severely affected patients has recently shown to yield more robust results (47). The particular strengths of the present study are (i) the inclusion of two independent datasets, each obtained with a different acquisition scheme, and (ii) the use of CSD tractography, which better enables resolving crossing fibers than DTI tractography and thus results in less falsenegative streamlines.

CONCLUSION
From a qualitative point of view, the connection matrices among the groups in two data sets showed a general similarity, which reinforces the confidence in the results of subsequent quantitative analyses. The negative correlation between E glob /E loc and the UPDRS scores in Data-NL suggested specific PD-related network differences for posterior parietal and sensorimotor cortical regions. Particularly, the observed parietal effect reinforces current reports on the involvement of this posterior brain region in PD. In this way, the present study provides leads to extend the search for biomarkers in PD using diffusion MRI complex network analysis. This includes the need to study longitudinal data, increase the sample size, and optimize acquisition parameters. Moreover, analyses need further improvements and validation before certainty can be reached about reduced integrative and segregative capacities of the brain affected by PD.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the University Medical Center Groningen and the University of British Columbia Ethics Board. The patients/participants provided their written informed consent to participate in this study.