Partial Similarity Reveals Dynamics in Brainstem-Midbrain Networks during Trigeminal Nociception

Imaging studies help us understand the important role of brainstem and midbrain regions in human trigeminal pain processing without solving the question of how these regions actually interact. In the current study, we describe this connectivity and its dynamics during nociception with a novel analytical approach called Partial Similarity (PS). We developed PS specifically to estimate the communication between individual hubs of the network in contrast to the overall communication within that network. Partial Similarity works on trial-to-trial variance of neuronal activity acquired with functional magnetic resonance imaging. It discovers direct communication between two hubs considering the remainder of the network as confounds. A similar method to PS is Representational Similarity, which works with ordinary correlations and does not consider any external influence on the communication between two hubs. Particularly the combination of Representational Similarity and Partial Similarity analysis unravels brainstem dynamics involved in trigeminal pain using the spinal trigeminal nucleus (STN)—the first relay station of peripheral trigeminal input—as a seed region. The combination of both methods can be valuable tools in discovering the network dynamics in fMRI and an important instrument for future insight into the nature of various neurological diseases like primary headaches.


Introduction
Functional magnetic resonance imaging (fMRI) can identify hubs of neuronal networks and as such simultaneously activate areas specific to an investigated task, but allows little to no access to the connectivity within these networks and particularly no insights into their dynamic during tasks. Nevertheless, exactly these network dynamics are probably an important key to understanding the pathophysiological mechanisms of certain diseases. A good example is cycling primary headaches like migraine and cluster headaches where the network dynamics between the hypothalamus and certain brainstem nuclei change during the cycling phases of the disease, i.e., from interictal to ictal phases, leading to altered nociceptive processing [1][2][3]. However, further analysis of such dynamics is limited since causal insights into connectivity like Granger causality are not reliable for most fMRI experiments due to the poor temporal resolution [4][5][6]. The most promising alternative, dynamical causal modelling (DCM) [7], is a model with numerous limitations and constrains. The repetition time needs to be faster than 2 s and a strong a priori knowledge of all the regions of interest, their possible We tested our approach using simulations presented in the Supplementary Document 2 and in brainstem-optimized fMRI data where 29 subjects underwent 15 trials of strong nociceptive trigeminal input into the left nostril as a painful condition and air puffs in the control condition. This leads to well-known activation of the trigeminal system [3,26,27], with the spinal trigeminal nucleus being the first relay station of the central nervous system for incoming trigeminal input [28]. Combining Representational Similarity with Partial Similarity, we were able to more reliably describe the trigeminal nociceptive network and additionally unravel direct connectivity between individual hubs of this network which may otherwise be hidden.

Simulations
To prove the translation of using partial correlations for beta-time series analysis, we ran simulations to evaluate the complementary information revealed by Ordinary and Partial Similarity. In summary, we simulated several amounts of global network signals mixed into 100,000 voxel with 435 time points, i.e., betas. The simulations include: (i) no direct connectivity between a seed and a target voxel, (ii) several levels of direct connectivity between a seed and a target voxel, and (iii) a situation where the seed and the target voxel are influenced by one or more regions of a global network signal which has excitatory influence on the seed but inhibitory influence on the target voxel. We present these simulations in the Supplementary Document 2 entitled "Simulations on the effect of a global network signal on Ordinary and Partial Similarity".

Subjects and Experimental Design
Twenty-nine healthy volunteers participated in our experiment on trigeminal nociception which has been repeatedly shown to prompt robust results (see Schulte et al. [29][30][31]). The current data was partly already published [29][30][31]. In short, participants received 15 trials of gaseous ammonia as nociceptive trigeminal input and 15 air puffs as a control condition. Each trial was rated for intensity and unpleasantness. The standardized experiment also included 15 trials of rose odor and 15 trials of visual stimulation not analyzed for this work. Written informed consent was obtained from all participants and the study was conducted according to the Declaration of Helsinki and approved by the Ethics Committee in Hamburg, Germany (PV 4522). The datasets generated and analyzed during the current study are not publicly available due to national data protection acts. Data are available from the corresponding author upon reasonable request.

Image Acquisition
MR-Images were collected on a 3T scanner (TRIO, Siemens, Munich, Germany) using a 32-channel head coil. Functional images were acquired with high-resolution EPI optimized for the human brainstem (38 axial slices, 1.25 × 1.25 × 2.5 mm 3 , TR 2.61 s, TE 27 ms, FOV 216 mm 2 , GRAPPA accelerated, 2 saturation pulses) [29]. The measured volume was restricted to the brainstem including the Thalamus and parts of the Cingulate Gyrus as the top border and the foramen magnum as the bottom border. The cerebellum was only partly within the measured volume. High-resolution structural images were obtained using an MPRAGE sequence with 1 mm 3 isotropic resolution [29].

Preprocessing
Functional and structural images were denoised using a spatially adaptive non-local means filter [32] as implemented in the CAT12 toolbox (http://www.neuro.uni-jena.de/cat/). Functional images were then realigned and slice time corrected using SPM12 (Wellcome Trust Center for Neuroimaging, London, UK). The average EPI was further co-registered to the participants' structural image. To realize our analytical approach, we calculated one General Linear Model (GLM) where each painful and each control trial was individually modeled with an HRF and included them as regressors of interest in a trial-by-trial GLM [13,15], while visual and odor condition were included as single, condition-wise regressors. Furthermore, we included models of button presses, the 6 movement regressors obtained in the realignment step within the preprocessing, and the discrete cosines (DCT), used in the standard SPM analyses as a high pass filter (128 Hz), in the GLM analysis. For physiological noise correction we additionally included 18 to 20 regressors extracted from the subjects' breathing and pulse signals with the approach described by Deckers and colleagues [33]. Resulting trial-specific beta-images were then normalized to MNI space with an isotropic voxel size of 2 mm 3 using the segmentation of the participants' structural image as implemented in SPM12 and smoothed by a 2 mm 3 isotropic Gaussian kernel and are the basis of all further analyses.
We then compared the trial-by-trial GLM results [13] to previously published results on condition-wise GLMs [2,29,34], calculating a contrast image with all trials on the subject level which entered a group level one-tailed t-test with an uncorrected threshold of p < 0.0001.
Trial-by-trial beta-images were further subject-wise subtracted by their mean and concatenated across subjects (Figure 1). We used a mask for gray and white matter gained from segmenting the average of the spatially-adaptive non-local mean corrected and normalized structural images. Furthermore, all voxels with zero variance were rejected, leaving about 100,000 (precisely 96,973) voxels for the proposed analysis.
Finally, the weight w used in Formula 2 can then be estimated with simple linear regression for x and y, respectively, using = min ( − 〈 , 〉 ) .

Application of Similarity and Partial Similarity
In the presented analyses we chose the left spinal trigeminal nucleus (STN) as a seed region as it is the first relay station in the central nervous system for trigeminal nociceptive input [28]. As the center of the seed, we used MNI coordinates (−4, −46, −53) from a published independent sample [34] and constructed a sphere around it with a radius of 8 mm. We further excluded the cerebellum from this sphere by masking it with the template provided by the SUIT toolbox for SPM 12 [38] normalized to MNI space.
The target area is constructed as a sphere with a radius of 8 mm (like the seed) and wanders around the entire volume using a searchlight approach [39] such that each voxel is the center of the target in exactly one iteration. The spatial average is used for calculating Similarity and Partial Similarity.
To perform Partial Similarity, the activity of the volume of no interest is first reduced along the variance into its first 15 components explaining a minimum of 30% of the variance (100 randomly chosen voxel were used for this estimation) by PCA as described above. Due to the high homogeneity between nearby areas, we defined an additional space with a radius of 15 mm around seed and target regions which were excluded from the analyses to prevent partial correlation coefficient of only −1 and 1. The proposed analyses are sketched in Figure 1. The searchlight approach used for the analysis involves the STN as seed (S, blue) whereas the target (T, red) represents a sphere wandering (represented by the black arrow) through the volume in a searchlight approach. Any signal that is detected in all volumes (such as movement artefacts) is used for correction in the Partial Similarity analysis. The plot bottom right shows three possible time courses for seed, target and general signal shared by the whole volume (V).

Statistics
The p-values stemming from the calculated Spearman correlations and partial correlations are corrected for multiple comparison using the Benjamini-Hochberg [40] and the Benjamini-Yekutieli [41,42] procedure for controlling the false discovery rate (FDR) of a family of hypothesis tests [43]. All reported statistical tests are one-tailed. The searchlight approach used for the analysis involves the STN as seed (S, blue) whereas the target (T, red) represents a sphere wandering (represented by the black arrow) through the volume in a searchlight approach. Any signal that is detected in all volumes (such as movement artefacts) is used for correction in the Partial Similarity analysis. The plot bottom right shows three possible time courses for seed, target and general signal shared by the whole volume (V).

Plausibility Checks
To gain first insights into the meaning of the variance of trial-by-trial activity we ran a Principal Component Analysis (PCA) on the preprocessed, concatenated and standardized data. A PCA finds uncorrelated components of the data, which are sorted by their explained variance, and is used to reduce high dimensional data. The first component is therefore the most meaningful compartment of the data. Its weights explain which voxels contribute to which extent to the component and give important insights into the data quality and hint its meaning [35][36][37].

Formulation of Similarity and Partial Similarity
As a measure of connectivity between a seed area x and a target area y, we calculated Spearman correlation coefficients between their spatial means referred to as (representational) Similarity [21]. Spearman correlation can be understood as a special case of classical Pearson's correlation using the rank rg of the sorted entries of x and y and can be formulated as where COV is the covariance matrix and s is the standard deviation. We further introduce Partial Similarity, which is also based on the Spearman correlation but corrects for matrix C by subtracting a weighted, linear product of C before calculating the correlation coefficient, i.e., performs a partial correlation. This can be formulated by where x and y are again the spatially averaged beta-estimates of seed and target and w weights the vectors in C. In C we find the first components from the transform along the variance of the standardized beta-images of V (Volume of No Interest) by a PCA. This transformation is used for data reduction which becomes necessary because the number of concatenated beta-images limits the number of components, which can be used as controlling variables within the partial correlation, i.e., more controlling variables than the number of elements in x (and thereby y) are not possible. C is formulated as where V are the standardized entries of all voxels with the Volume of No Interest (sketched in Figure 1).
The first component C 1 of V is constructed by finding its weight b 1 according to and the further components C i iteratively by finding their weights b i using an altered V in Formula 4 with the already calculated components subtracted from V. Finally, the weight w used in Formula 2 can then be estimated with simple linear regression for x and y, respectively, using

Application of Similarity and Partial Similarity
In the presented analyses we chose the left spinal trigeminal nucleus (STN) as a seed region as it is the first relay station in the central nervous system for trigeminal nociceptive input [28]. As the center of the seed, we used MNI coordinates (−4, −46, −53) from a published independent sample [34] and constructed a sphere around it with a radius of 8 mm. We further excluded the cerebellum from this sphere by masking it with the template provided by the SUIT toolbox for SPM 12 [38] normalized to MNI space.
The target area is constructed as a sphere with a radius of 8 mm (like the seed) and wanders around the entire volume using a searchlight approach [39] such that each voxel is the center of the target in exactly one iteration. The spatial average is used for calculating Similarity and Partial Similarity.
To perform Partial Similarity, the activity of the volume of no interest is first reduced along the variance into its first 15 components explaining a minimum of 30% of the variance (100 randomly chosen voxel were used for this estimation) by PCA as described above. Due to the high homogeneity between nearby areas, we defined an additional space with a radius of 15 mm around seed and target regions which were excluded from the analyses to prevent partial correlation coefficient of only −1 and 1. The proposed analyses are sketched in Figure 1.

Statistics
The p-values stemming from the calculated Spearman correlations and partial correlations are corrected for multiple comparison using the Benjamini-Hochberg [40] and the Benjamini-Yekutieli [41,42] procedure for controlling the false discovery rate (FDR) of a family of hypothesis tests [43]. All reported statistical tests are one-tailed.

Simulations
The first simulation shows that a global network signal induces a strong connectivity between seed and target as measured by Ordinary Similarity analysis which can reach correlation values close to 1 when the influence of the global network signal is extreme ( Figure S2). The Partial Similarity analysis shows, as expected, that there is no direct connectivity between seed and target ( Figure S2). As the maximally measured correlation coefficients in the actual data reaches levels of around r = 0.2 and the explained variance of the controlling variables is around 30%, the global network signal might have an influence of up to 60% in our real fMRI beta time series data.
The second Simulation shows the complementary information of Ordinary and Partial Similarity: While the Ordinary Similarity reveals connectivity induced either by the global network signal or direct connectivity or a mixture of both, Partial Similarity only reveals the direct connectivity component.
Simulation 3 tests to mislead Ordinary Similarity by adding a global network signal to the seed but subtracting it from the target. Here, Ordinary Similarity can reveal even negative correlations or near-zero correlations when the magnitude of the global signal and the direct connectivity is similar. Details of the simulations can be found in the Supplementary Document 2.

Plausibility
Trial-by-trial GLM analysis confirmed previously published results [2,26,27,29,34,[44][45][46], i.e., significant activation in left STN, rostral Pons, Cerebellum, periaqueductal gray (PAG), Thalamus, and bilateral Insula (Figure 2A). Due to the restricted measured volume, no activation could be shown in other areas. The first component of the PCA on the concatenated trial-by-trial beta-images from the GLM explains 9.85% of the overall variance (voxel × trials) and reveals a plausible network within the brainstem during trigeminal nociception. The network includes previous reported hubs of trigeminal nociception in the brainstem and midbrain such as the STN, rostral Pons, Cerebellum, PAG, bilateral Thalamus and bilateral Insula ( Figure 2B). This unsupervised analysis underlines the potential of analyzing trial-by-trial variance.

Similarity and Partial Similarity
Similarity analysis from the STN ROI as a seed reveals three clusters with positive correlations and 16 clusters with negative correlations at a statistical threshold of FDR-adjusted p < 0.05, where the most extensive cluster contains more than 54% of all considered voxel, especially most of the midbrain. To pinpoint more specific regions, we increased the statistical threshold to an FDR-adjusted p-value of p < 0.0001, revealing nine clusters with similar trial-to-trial activity changes as the STN in the nociceptive condition ( Figure 3A) with a minimum cluster extent of 5 voxel (40 mL). The clusters mostly include cortical and subcortical structures already shown to be involved in trigeminal nociception like bilateral Cerebellum, bilateral PAG, bilateral Insula, bilateral Operculum, right Striatum (Caudate, Putamen), and bilateral Thalamus and furthermore the right Cingulate Gyrus (Table 1).
in other areas. The first component of the PCA on the concatenated trial-by-trial beta-images from the GLM explains 9.85% of the overall variance (voxel × trials) and reveals a plausible network within the brainstem during trigeminal nociception. The network includes previous reported hubs of trigeminal nociception in the brainstem and midbrain such as the STN, rostral Pons, Cerebellum, PAG, bilateral Thalamus and bilateral Insula ( Figure 2B). This unsupervised analysis underlines the potential of analyzing trial-by-trial variance.

Similarity and Partial Similarity
Similarity analysis from the STN ROI as a seed reveals three clusters with positive correlations and 16 clusters with negative correlations at a statistical threshold of FDR-adjusted p < 0.05, where the most extensive cluster contains more than 54% of all considered voxel, especially most of the midbrain. To pinpoint more specific regions, we increased the statistical threshold to an FDRadjusted p-value of p < 0.0001, revealing nine clusters with similar trial-to-trial activity changes as the STN in the nociceptive condition ( Figure 3A) with a minimum cluster extent of 5 voxel (40 mL). The clusters mostly include cortical and subcortical structures already shown to be involved in trigeminal nociception like bilateral Cerebellum, bilateral PAG, bilateral Insula, bilateral Operculum, right Striatum (Caudate, Putamen), and bilateral Thalamus and furthermore the right Cingulate Gyrus (Table 1. The control condition with air puffs also showed highly correlated midbrain areas at a threshold of FDR-adjusted p < 0.05. Regarding nociception, we therefore increased the threshold to p < 0.0001 ( Figure 3B, Table 1) and identified seven clusters correlated with the left STN seed region in their trial-to-trial variance. These clusters included the bilateral Cerebellum, bilateral PAG, the left inferior Temporal Gyrus, left Fusiform Gyrus, and the right Planum Polare (Table 1). The control condition with air puffs also showed highly correlated midbrain areas at a threshold of FDR-adjusted p < 0.05. Regarding nociception, we therefore increased the threshold to p < 0.0001 ( Figure 3B, Table 1) and identified seven clusters correlated with the left STN seed region in their trial-to-trial variance. These clusters included the bilateral Cerebellum, bilateral PAG, the left inferior Temporal Gyrus, left Fusiform Gyrus, and the right Planum Polare (Table 1).     Table 2). Most of the regions are well known to process painful trigeminal input.
In the control condition, we observed four clusters with positive correlations and one cluster with negative correlations to the STN ( Figure 4B, Table 2), namely bilateral Cerebellum and the Temporal Pole ( Figure 5).   The included regions are bilateral Cerebellum, PAG, and Cingulate  Gyrus, parts of the temporal lobe (right Planum Temporale, right Middle Temporal Gyrus, left Inferior  Temporal Gyrus), bilateral Striatum (Caudate), and with a negative correlation the right Fusiform Gyrus and both sides of the Thalamus ( Figure 5, Table 2). Most of the regions are well known to process painful trigeminal input.

Discussion
Uncovering cortical and subcortical connectivity in human fMRI data is a challenging field as to date only a few methods exist. While commonly used DCM and PPI both work on raw time courses, they depend on the height of BOLD amplitude, and are further hindered by a long list of assumptions and constraints, ordinary correlation analysis may lead mostly to an overestimation of results. Here we introduce a new way to delineate direct (dynamic) network connectivity complementary to ordinary correlation analyses of beta time series: Partial Similarity.
We conceptually proved that Partial Similarity and its counterpart, (based on ordinary correlations) Ordinary Similarity, reliably reveal dynamic connectivity of brainstem networks known from preclinical studies. As a physiological model, we specifically chose the connectivity of the spinal trigeminal nucleus, which is the first central hub of all incoming peripheral trigemino-nociceptive input into the central nervous system (CNS). Studies on mammals suggest direct anatomical connections to the cerebellum [47][48][49], Thalamus [50], and PAG [50][51][52], which are known to be involved in mechanical and nociceptive input [49,51,53,54]. Additional to these well-known hubs, In the control condition, we observed four clusters with positive correlations and one cluster with negative correlations to the STN ( Figure 4B, Table 2), namely bilateral Cerebellum and the Temporal Pole ( Figure 5).
To gain an overview of the reported results, we sketched the revealed connectivities and their strength from the Similarity and Partial Similarity analyses in Figure 5. While both connectivity analyses show strong connectivities to Cerebellum and PAG during nociception and reveal the direct connectivity from STN to (contralateral) Thalamus, the strong contralateral crossing to higher cortical areas (Temporal lobe, Insula, Fusiform Gyrus) is only revealed using Partial Similarity analysis. In the control condition, the Similarity analyses suggest a network between left STN and bilateral Cerebellum and bilateral PAG. The latter could not be identified in the Partial Similarity analysis.

Discussion
Uncovering cortical and subcortical connectivity in human fMRI data is a challenging field as to date only a few methods exist. While commonly used DCM and PPI both work on raw time courses, they depend on the height of BOLD amplitude, and are further hindered by a long list of assumptions and constraints, ordinary correlation analysis may lead mostly to an overestimation of results. Here we introduce a new way to delineate direct (dynamic) network connectivity complementary to ordinary correlation analyses of beta time series: Partial Similarity.
We conceptually proved that Partial Similarity and its counterpart, (based on ordinary correlations) Ordinary Similarity, reliably reveal dynamic connectivity of brainstem networks known from preclinical studies. As a physiological model, we specifically chose the connectivity of the spinal trigeminal nucleus, which is the first central hub of all incoming peripheral trigemino-nociceptive input into the central nervous system (CNS). Studies on mammals suggest direct anatomical connections to the cerebellum [47][48][49], Thalamus [50], and PAG [50][51][52], which are known to be involved in mechanical and nociceptive input [49,51,53,54]. Additional to these well-known hubs, Similarity further reveals strong connectivity to the ipsilateral Insula and Operculum, which was not observed with Partial Similarity. Partial Similarity on the other hand emphasizes connectivity to contralateral regions like the striatum, Fusiform Gyrus, Temporal lobe, and ipsilateral Cingulate Gyrus. These differences stem from the nature of the two measures, which we demonstrated in our simulations: while Similarity shows overall communication between two regions, Partial Similarity focuses on direct communication between two hubs, diminishing trivial and induced correlations, but uncovering correlations hidden by other sources or a global network signal. Both analyzing techniques contain complementary information (direct communication and overall communication), which should be exploited in future research by locating the origin of these differences. We suggest that imaging studies, which investigate relationships in and between neuronal networks, will strongly profit from the combination of both analytical tools. This is particularly true for research investigating cycling pathologies like migraine and cluster headaches [55], where changes in network dynamics [1] follow the current stage of the patients' disease [2,3,56]. While we concentrate on the first hub of the trigeminal driven network in the CNS, further research will reveal the dynamical connectivity of all hubs in the network and may identify cortical and subcortical structures by combining both analytical tools, which would be missed if only one of these tools were used. Moreover, given the open-ended search feature of the proposed analyses, we aim to identify connectivity and its dynamics and hubs, which is not yet within the focus of current research. Further work will have to include the translation of Partial Similarity to other fields of Representational Similarity applications [21], where trivial correlations currently might hinder a deeper insight into direct relations, for example in the field of complementing fMRI and electroencephalography (EEG) [57]. Partial Similarity and Similarity analyses as established in the current work might thus be an important tool to discover neuronal network dynamics in fMRI and thus grant more detailed insights into the pathophysiological mechanisms underlying neurological diseases, including but not limited to primary headache disorders.

Conclusions
We have shown with both simulated and real task-dependent fMRI data that the Partial Similarity method presented here can represent direct connectivity between two hubs. Specifically, Partial Similarity can unravel brainstem dynamics involved in trigeminal pain. It thereby gains complementary information to Ordinary (also called Representational) Similarity, which reveals a summary of global and direct connectivity. The combination of both methods can be valuable tools in discovering network dynamics in fMRI and an important instrument for future insight into the nature of pathological pain processing and various other neurological diseases.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3425/10/9/603/s1: The Supplementary Materials include detailed descriptions and results of the cases and simulations referred to in the article. It consists of: Document 1: "Cases of Ordinary and Partial Similarity" including: Pseudocode P1: Code used to construct the 4 sample cases for ordinary and partial correlations; Figure S1: Simulated cases where ordinary and partial correlations show the same (Case 1 and 2) and divergent results (Case 3 and 4). Document 2: "Simulations on the effect of a global network signal on Ordinary and Partial Similarity " including: Pseudocode P2: Code used to conduct the simulations on the effect of a global network signal on (non-existing) direct connectivity between seed and target revealed by ordinary and partial correlations.; Pseudocode P3: Code used to conduct the simulations on the effect of a global network signal on several levels of direct connectivity between seed and target revealed by Ordinary and Partial Similarity; Figure S2: Results of the simulation of the Ordinary and Partial Similarity analysis of the simulated beta time series for 101 steps of 10% increment of global network influence while no direct connectivity between seed and target is simulated; Figure S3: Results of the simulation of the Ordinary and Partial Similarity analysis of the simulated beta time series for 11 steps of 10% increment of global network influence while direct connectivity between seed and target is also simulated for 11 steps of 10% increment; Figure S4: Results of the simulation of the Ordinary and Partial Similarity analysis of the simulated beta time series for 11 steps of 10% increment of global network influence. Here the global network signal has an additive effect on the seed and a negative effect on the target voxel.