Participants
Seven neurotypical participants (mean age: 24.9 ± 2.1 years; four females, three males) were recruited via word-of-mouth advertisement. All participants were right-handed, with no self-reported injuries or sensory impairments to their hands, and no self-reported history of neurological disorders. Participants received €30 for their participation. The experiment was approved by the local ethics review committee of Maastricht University and participants gave written informed consent about their participation in the experiment.
Study procedure
In this within-subject design, participants attended two MRI sessions (separated by 14.3 ± 7.4 days; Range: 7–24 days). No specific instructions were given for the time between the two sessions; it was assumed that participants continued their regular activities. Each 2-hour-session consisted of one structural MR acquisition and four functional digit mapping measurements. Eight additional functional MR measurements obtained in each session are not considered in the present analysis. These other measurements included mapping of the non-dominant hand, both hands at the same time, and active finger tapping sequences. Since the focus of the present analysis is on retest reliability, both sessions included the same measurements.
Mapping procedure
During the functional measurements, participants were asked to keep their eyes open and to fixate on a cross back-projected to the center of a screen at the end of the scanner bore. At the same time passive vibrotactile stimulation was applied through five separate modules of a mini piezotactile stimulator (mPTS; Dancer Design, Merseyside, United Kingdom; Fig. 1B) attached to the most distal phalanx of the five digits of the right hand. A stimulation frequency of 25 Hz was delivered through a metal probe (diameter: 6 mm) positioned centrally at the top of the module, moving approximately 0.5 mm up and down. This frequency corresponds to the flutter range, optimally activating the rapidly adapting type 1 afferent fibers [51]. Order and timing of the stimulation was defined by two different stimulation designs, the traveling wave (TW) and the blocked design (BD), as specified below.
Stimulation device and stimulation designs. A Traveling wave (TW) and blocked design (BD) consist of two stimulation directions: Forward (Fw; upper row) with stimulation moving from D1 to D5 and Backward (Bw; lower row) with stimulation moving from D5 to D1. The TW consists of 4 s stimulations of each digit, repeated 15 times in both stimulation directions. The BD consists of 12 s stimulations of each digit and a 12 s non-stimulation rest period (gray) after every fourth digit stimulation. TW was analyzed using cross correlation with two predictors per digit shifted by one TR (2 s). BD was analyzed using a standard GLM analysis with one predictor per digit. B Mini piezotactile stimulators (mPTS) with a round, MRI-compatible, metal probe (left panel), delivering a 25 Hz vibrotactile stimulation, were attached to the most distal phalanx of all five digits of the participants right hand (right panel). C The order of the stimulation designs within the two sessions was counterbalanced across the participants by assigning them to either Version A or B. In each version, both sessions began with the Fw runs of the two designs, starting with either BD or TW (counterbalanced across sessions). Then, the Bw runs were acquired in the opposite order. During the first session, multiple other measures not relevant for the current experiment (MoM) were acquired between the Fw and Bw runs. In the second session, MoM were acquired after the Bw runs. Lastly, the T1-weighted scans (T1w) were always preceded and followed by functional scans with reversed phase encoding direction (RP), required for EPI distortion correction of the functional scans of interest. All sessions ended with the acquisitions of MoM
Traveling wave
During the measurement runs with TW stimulation (Fig. 1A), the five digits of the right hand were repeatedly stimulated in successive anatomical order: thumb (D1), index finger (D2), middle finger (D3), ring finger (D4), little finger (D5). Two different directions of stimulation were acquired in separate measurement runs: forward (D1-D2-D3-D4-D5-D1-D2-D3….) and backward (D5-D4-D3-D2-D1-D5-D4-D3…). For one digit, each stimulation period lasted 4 s, resulting in a TW-cycle of 20 s for all five digits. Each forward or backward TW fMRI measurement run consisted of 15 cycles, resulting in a total acquisition time of 5 min and 20 s, including 10 s non-stimulation periods at the beginning and the end. The combined acquisition time of both TW stimulation directions was 10 min and 40 s, i.e. each digit was stimulated for 120 s across both stimulation runs. The forward and backward stimulation directions were acquired in a counterbalanced fashion (Fig. 1C).
Blocked design
The BD (Fig. 1A) was based on the approach used by RS and colleagues [52, 53, 56]. Digits were stimulated for 12 s in anatomical succession (from D1 to D5), but, contrasting the TW approach, with non-stimulation rest periods included after every fourth digit stimulation. This resulted in a long BD cycle in which each digit was stimulated four times and five rest periods were incorporated after every fourth digit stimulation (D1, D2, D3, D4, rest, D5, D1, D2, D3, rest, D4, D5, D1, D2, rest, D3, D4, D5, D1, rest, D2, D3, D4, D5, rest), i.e. once before and once after the stimulation of each of the five digits. Based on the 12 s stimulation periods for each digit stimulation and the non-stimulation periods, one BD-cycle lasted 5 min. The two stimulation directions were acquired in a counterbalanced fashion (Fig. 1C).
The long stimulation of 12 s, as well as the non-stimulation rest periods, are both elements taken from the classical block design [46, 58]. As the stimulation of the single digits occurs in successive anatomical order, the sparse number of rest periods also follows a specific”wandering” pattern of changing position in an orderly manner. Being placed at least once before and after each digit stimulation allows a more balanced general linear model for the analysis of the associated BOLD-activation of each digit.
The BD runs of the first two participants were matched in acquisition time to the TW run of 5 min and 20 s consisting of one BD cycle (5 min) plus a non-stimulation period of 10 s at the start and end, resulting in a total acquisition time of 10 min and 40 s for combined forward and backward stimulation. The total stimulation time per digit was 96 s total compared to 120 s total stimulation time in the TW design.
In the remaining five participants the BD paradigm was lengthened to compensate for the unequal total digit stimulation time. Five 12 s stimulations, one for each digit, plus a non-stimulation rest period after the fourth digit stimulation were added at the end of the runs, resulting in the total duration of stimulation of 120 s for each digit. This increased the acquisition time of the BD to 6 min and 40 s per run (8 s baseline, 5 min BD cycle plus 1 min and 12 s stimulation and a 20 s non-stimulation period after the last digit stimulation) and the total acquisition time to 13 min and 20 s for both stimulation directions.
Attentional task
Since attention to touch increases somatosensory cortical activation [31] and counteracts potential habituation effects, short time segments without stimulation were inserted within the period of the vibrotactile stimulation [3, 52,53,54, 57]. In the TW design, seven 100 ms long interrupts were included in every 4 s digit stimulation spaced by 500 ms. Since the interrupts were rather short and difficult to perceive, participants were asked to count the total number of TW stimulation cycles to assure that their attention was focused on the stimulation.
In the BD, several short non-stimulation periods of 150 ms were included within the 12 s stimulation of each digit. The number of interrupts and their onset within the stimulation period were constant for each digit stimulation, but differed across stimulations of different digits. The stimulation period of D1 and D5 included five interrupts, D2 and D3 six interrupts, and D4 four interrupts. Participants had to count the total number of interrupts in the vibrotactile stimulation during each run of the BD [52,53,54, 57] and verbally report this number at the end of each run via the intercom used for communication.
Data acquisition
MRI data were obtained at a Siemens 3 T Prisma Fit system with a 64-channel head coil (Siemens Healthcare, Erlangen, Germany). Anatomical T1-weighted data were acquired with the 3D whole brain coverage Alzheimer’s Disease Neuroimaging Initiative (ADNI) Magnetization Prepared Rapid Gradient Echo (MPRage) sequence (TR = 2300 ms; TE = 2.98 ms; flip angle = 9°; Bandwidth = 240 Hz/Px; FoV = 256 × 256 mm; number of slices = 192; spatial resolution = 1 × 1 × 1 mm3). Functional data were acquired using T2*-weighted Gradient Echo Echo Planar Imaging (GE-EPI) (TR = 2000 ms; TE = 30 ms; flip angle = 77°; Bandwidth = 1786 Hz/Px; Multiband Acceleration Factor = 2; GRAPPA Factor = 2; FoV = 200 × 200 mm; number of slices = 64; spatial resolution = 2 × 2 × 2 mm3). An additional functional measurement of two volumes with reversed phase encoding direction (posterior to anterior) was recorded before and after the structural measurement (Fig. 1C), to be used for EPI distortion correction during preprocessing [10]. For two participants these reversed phase encoding data were, due to technical difficulties, not available, so no distortion correction could be performed. Careful inspection of the T2* weighted images showed no constraints for the further analysis of the data of these two subjects, since the digit area of the somatosensory cortex being located mid brain is usually only minimally impacted by the distortion.
To minimize head motion as well as to assure participants’ comfort, foam padding was used inside the head coil. To further reduce head motion during data acquisition medical tape was attached from one side of the head coil across the participants’ forehead to the other side of the head coil so that even small head movement resulted in a slight pull on the participant’s skin, providing detailed tactile feedback, in addition to the minimal fixation of the head [37].
Data analysis
MRI preprocessing and the GLM-analysis of the BD was carried out using BrainVoyager (Version 22.0, Brain Innovation, Maastricht, The Netherlands). The cross-correlational analysis of the TW design as well as the feasibility and reliability measures were performed using the NeuroElf toolbox (Version 1.1) for Matlab (Version 2020a, MathWorks).
Data preprocessing
Structural MR images were intensity inhomogeneity corrected, non-linearly transformed into MNI space (ICBM-MNI 152) and for each participant averaged across the two sessions. Functional MRI data were slice-time corrected (cubic spline interpolation), as well as motion corrected and aligned (detection: trilinear; correction: sinc interpolation) to the run and volume that was temporally closest to the anatomical scan during each session. To correct for low-frequency noise, the functional data were high-pass filtered with a cut-off frequency of 0.01 Hz. EPI distortion correction was obtained with the BrainVoyager plugin COPE 1.1 and to preserve the 2 mm isotropic spatial resolution, no extra spatial smoothing was performed.
Coregistration of functional MRI data to within-session structural MRI scans was achieved with boundary-based registration and transformation of functional data into MNI space was performed with sinc interpolation; in the same step, the spatial resolution of the functional data was interpolated to 1 mm isotropic resolution.
Definition of the region of interest
The region of interest (ROI) was defined following the procedure described by Valente and colleagues [63], taking inter-individual anatomical differences into account without the need to manually draw ROIs for each participant.
Individual MNI-normalized cortical surface meshes of the left-hemispheric white–gray matter boundary were created and aligned to a standard sphere using cortex based alignment (CBA) in BrainVoyager [22, 26]. An earlier study has shown a significant and meaningful improvement of the macro-anatomical correspondence of the primary somatosensory hand region between participants after CBA as compared to a volumetric normalization approach [24]. For CBA, a dynamic group average, representing the average curvature information of the specific sample, is created by an iterative alignment of the individual curvature maps from each participant’s hemisphere. This averaged and cortex-based aligned cortical mesh of the seven participants was used to draw a ROI manually, covering the anterior wall of the postcentral gyrus opposite of the hand knob area [69] stretching from the fundus of the central sulcus towards the crown of the postcentral gyrus (see Fig. 2). This ROI is expected to cover the hand area of somatosensory BA 3b as well as parts of BA 3a and BA 1, while reducing the risk of including the large blood vessels over BA 1 [23]. This ROI was subsequently back projected to each participant’s cortical surface mesh using the information from the CBA procedure and in a second step sampled from the surface space to the volume space by expanding it by 2 mm both towards the white matter and the cerebrospinal fluid to assure coverage of the entire gray matter ribbon in that region.
Average and individual regions of interest. The region of interest (ROI) was defined on the averaged surface mesh (top left) as the posterior bank of the central sulcus opposite to the motor hand area. Subsequently, it was backprojected to each participant’s individual surface mesh, covering a similar area in all participants
First level analyses
Cross-correlational analysis of the traveling wave design
The TW was analyzed based on the cross-correlation approach outlined by Kolasinski and colleagues [34, 35] and originally used in retinotopy [19]. A reference predictor was created using a boxcar function including a 4 s “on” period, reflecting the digit stimulation and a 16 s “off” period, reflecting the stimulation of the other digits, which was repeated 15 times to reflect the TW stimulation cycles. This reference predictor was then iteratively shifted 10 times by 2 s, i.e. the TR of the functional measurement, resulting in 10 predictors. To account for the hemodynamic delay each predictor was convolved with a two-gamma canonical response function (Onset = 0,Time to response peak = 6 s; Response dispersion = 1; Response Undershoot ratio: 6; Time to undershoot peak = 16 s; Undershoot dispersion = 1).
Next, the time courses of all voxels within the ROI were correlated with all 10 time-shifted predictors. To allow for statistical tests of the resulting correlation coefficients and correct for their non-normal sampling, Fisher z-transformation was applied [59]. Each of the 10 predictors was then assigned to one of the five fingers, i.e. two predictors to each finger according to the boxcar predictor’s on-period. This approach resulted in two different HRF latencies per digit, i.e., 6 and 8 s (Fig. 1A) [31, 32, 35].
Within each voxel, the Fisher z-transformed correlation values of the two predictors assigned to the same digit were averaged, as were the correlation values from the forward and backward run. This resulted in five correlation values within each voxel, each representing the correlation to one of the five digits. Values in these maps that exceeded a false discovery rate (FDR) corrected threshold (based on all voxels within the ROI) of q(FDR) < 0.05 were labeled as active. Individual FDR thresholds are reported in Table 1.
A two-way repeated-measures analysis of variance (rm-ANOVA) was performed to investigate systematic differences in FDR thresholds, including the factors SESSION and DIGIT. No significant difference in FDR thresholds between the two sessions could be observed (F(1, 6) = 0.2, p = 0.671). The five digits differed significantly in their FDR thresholds (F(2.1, 12.7) = 21.8, p < 0.001), with D3 and D5 exhibiting higher FDR thresholds.
As part of the adaptation to the clinical usage we did not apply the “winner map” approach in which each voxel is exclusively assigned to the digit with the highest statistical value [31, 32, 35]. Utilizing multiple statistical values per voxel enables the analysis of activation overlap between digit activations, facilitating a more thorough investigation of its occurrence and retest reliability. Moreover, this approach permits identifying and excluding voxels that may be influenced by blood vessel activation (see 2.5.3.3. Functional Activation Associated with Veins).
General linear model analysis of the blocked design
For the analysis of the BD, a general linear model (GLM) analysis was applied within the ROI in BrainVoyager. Each digit stimulation was modeled by a boxcar predictor, which was convolved with the same two-gamma canonical response model as used for the TW design (Onset = 0; Time to response peak = 6 s; Response dispersion = 1; Response Undershoot ratio: 6; Time to undershoot peak = 16 s; Undershoot dispersion = 1) (Fig. 1A). A fixed effects analysis including both the forward and backward run was applied in each participant. Each single digit predictor was contrasted against the mean of the predictors of the four remaining digits. This contrast, delimiting the overlap of activation between digits, was introduced to counterbalance the distinguishably larger volumes of activation in the BD design, probably caused by the longer single digit stimulation periods in the BD (12 s) compared to TW (4 s). These larger areas of activation are presumably not reflecting a higher sensitivity to the area of the underlying neuronal digit representation but probably caused by an involvement of larger areas of the vascularity which is not expected to align with the discrete digit layout in BA 3b.
Voxels with statistical values that exceeded an FDR corrected threshold of q(FDR) > 0.05 were labeled as active in the somatosensory ROI. Individual FDR thresholds are reported in Table 2.
A two-way rm-ANOVA with the factors SESSION and DIGIT was applied on the FDR thresholds of the BD. No significant differences between the FDR thresholds of the two sessions could be detected (F(1, 6) = 0.2, p = 0.637). There was a significant difference between the FDR thresholds of different digits (F(2.5, 14.9) = 6.0, p = 0.009).
Functional activation associated with veins
To control for the influence of draining veins on the BOLD signal, single voxels that showed significant activity for three or more digits [57] were excluded from further analysis, both for TW and BD (Tables 1 and 2. This was done separately per session and design; thus, a single voxel that was excluded in one session or design was not necessarily excluded in the other session or design. The rationale behind this is that in a normal clinical setting, data from only one design and one session will be available and therefore should be sufficient to exclude all potential draining vein artifacts.
Visualization of digit maps
To visually inspect the somatosensory activation maps that were elicited by the passive vibrotactile stimulation of the participants’ digits, the data of each session and design were projected to the flattened cortical surface representation of each participant’s left hemisphere. The functional data were projected onto the flattened surface using trilinear interpolation. Data were sampled from a 4 mm ribbon around the reconstructed white–gray matter boundary covering approximately 1 mm into the white and 3 mm into the gray matter direction, using only the maximum value in this range. This was done solely for the purpose of visualization. All analyses were performed in volume space, except for the calculation of the geodesic distance measure, which were obtained from the vertices of the folded surface map.
Identification of single digit activation clusters
The activation of each digit in each design and session was defined as its largest cluster of significant BOLD-activation in volume space that contained the peak voxel (i.e., the voxel with the highest statistical value) within the ROI. In case the largest activation cluster did not contain the peak voxel, either the largest cluster or the cluster containing the peak voxel was chosen, depending on which cluster was closest to the clusters chosen for the neighboring digits.
Parameterization of digit activation clusters
Activation clusters were parameterized for comparison across the two designs (validity) and across the two sessions (reliability) using ‘center of gravity’ and ‘D1-D5 distance’ as parameters for the location of the activation as well as ‘volume of activation’ and ‘activation overlap between neighboring digits’ as parameters associated with the volume of activation.
Center of gravity
The center of gravity (CoG) was defined as the weighted center of activation of that digit’s activation cluster. The coordinates in volume space of the CoG were determined based on the cluster’s average coordinates weighted by each voxel’s statistical value [20]:
$${CoG}_{j}=\frac{{\sum }_{i=1}^{{n}_{j}}{t}_{i}{x}_{i}}{{\sum }_{i=1}^{{n}_{j}}{t}_{i}}$$
with x being the coordinate of voxel i in cluster j with size n and t being the statistical value of that voxel.
D1—D5 distance
The coordinates of the CoGs of D1 and D5 of each session and each design were projected onto the folded surface mesh of each participant and the geodesic distance, i.e., the shortest path along the surface between the vertex representing the D1 CoG and the vertex representing the D5 CoG was obtained using Dijkstra’s algorithm [15].
Even though the distance between D1 and D5 is a measure to describe the size of the cortical digit representation it is included as a location-based parameter, as it critically depends on the localization of the single digits’ activation.
Volume of activation
The volume of activation of the single digit representation was defined as the volume of the digit’s activation cluster within the anatomically predefined SI ROI, calculated in mm3.
Overlap of activation
The overlap of the activation clusters of two anatomically neighboring digits (D1 + D2, D2 + D3, D3 + D4, and D4 + D5) was calculated using the Dice Coefficient (DC) [14]:
$$Dice_{A, B} = \frac{{2 \times \left| {A \cap B} \right|}}{\left| A \right| + \left| B \right|}$$
where |A ∩ B| represents the activation volume shared between the two clusters A and B and |A| and |B| represent the activation volumes of the clusters A and B, respectively. It describes the volume of activation that is overlapping for the two digits relative to the additive volume activated for the two digits. The DC ranges between 0 and 1, with higher values indicating more overlap between the activation of the two digits.
Validity: statistical testing
To quantify the validity of mapping the somatosensory digit area using the two adapted mapping procedures, we statistically compared the four extracted map parameters of the first session across the two designs. We expect that both mapping designs produce similar results, similar to those reported in existing literature, and that they likely reflect the main aspects of the underlying somatosensory digit map.
The three coordinate values of the CoG were compared with three two factorial repeated measures (rm)-ANOVAs including the factors DESIGN (TW and BD) and DIGIT (D1 – D5).
The D1-D5 distance was compared across designs using a paired sample t-test.
The volumes of activation were compared with a two factorial rm-ANOVA including the factors DESIGN (TW and BD) and DIGIT (D1 – D5).
The overlap of activation of neighboring digits were compared with a two factorial rm-ANOVA with the factors DESIGN (TW and BD) and DIGITPAIR (D1 + D2, D2 + D3, D3 + D4, D4 + D5).
For all rm-ANOVAs the normality assumption was checked using the Shapiro–Wilk test (Table 3) and homoscedasticity was acknowledged by considering the Greenhouse–Geisser corrected statistical values.
Retest reliability: statistical testing
Stability of activation between sessions
To investigate the reliability of the activation pattern of both designs, the DC was utilized to compare the activation patterns across the two sessions by describing two different types of spatial correspondence. Firstly, the spatial correspondence of same single digit activations across the two sessions, with high DC values indicating high retest-reliability of that single digit cluster, both in terms of location and volume, across the two sessions. Secondly, the spatial correspondence of the overlap between neighboring digit activations across sessions. Here, high DC values indicate high retest-reliability of the overlap of neighboring digit activations across the two sessions both in terms of location and volume.
Correlational analysis of map parameters
To investigate the retest reliability of the map parameters, they were correlated across the two sessions separately for each mapping design using Pearson correlation. CoG data were correlated separately for each axis (X, Y, and Z), digit, and design, resulting in 30 individual correlation values (3 axes × 5 digits × 2 designs). The D1-D5 distance was correlated separately per design, resulting in two correlation values. The volume of activation was correlated separately for each digit and design, resulting in 10 individual correlation values (5 digits × 2 designs). The overlap of activation was correlated separately for each digit-pair and design, resulting in 8 individual correlation values (4 digit-pairs × 2 designs). Thus, in total, 50 correlation values were calculated and tested whether they were significantly larger than zero (one-tailed test). To correct for the multiplicity problem, FDR correction using a linear step-up procedure was used [2], values that exceeded a threshold of q(FDR) < 0.05 were considered significantly larger than zero. Additionally, for each correlation value, the slope and intercept of the best fitting line was calculated using the least squares method. The best fitting line of a highly reliable design is expected to be close to the line of equality (i.e., 1x + 0).
The individual stability of the CoG locations between the first and second session was described through Euclidean distances. The Euclidean distance was preferred over the surface-based geodesic distance as no transformation to the gray-white matter boundary is necessary, potentially projecting the very small inter-session differences onto the same vertex. Another advantage is the option to relate the inter-session distances between the CoGs to the voxel dimensions of the functional measurement.
Add Comment