Acupuncture stimulation and paradigm
Each fMRI scan run was 31.5 minutes in duration, during which the subject experienced one of two different stimulations: verum acupuncture (VA) or sham point acupuncture (SPA), which controlled for the location of stimulation. Each fMRI scan run began with a 90sec rest block period followed by 15 120sec periods of alternating stimulation (60sec) and rest (60sec) ().
Figure 1 Experimental protocol. (A) Stimulation location for both the verum acupoint (ST-36) and sham point (SH-1). (B) Cardiac-gated fMRI was performed using coronal slices covering the brainstem, as well as other cortical and subcortical structures in this field-of-view. (more ...)
For VA, electro-acupuncture stimulation was used with alternating 3-second periods of 2Hz and 15Hz electro-stimulation (known as “dense-disperse”) at acupoint ST-36 on the lower leg (). For SPA, electro-stimulation was performed with the same parameters as for VA at a sham acupoint, ~8cm above the proximal edge of the patella, on the midline of the thigh (). The reference electrode for VA was 1cm distal to ST-36, while for SPA, it was 1cm proximal to the sham point. Acupuncture needles were inserted before the start of each fMRI run in either the right or left leg (pseudo-randomized). Different legs were chosen for VA and SPA to avoid the possibility that long duration electro-stimulation altered segmental neural reactivity. Due to the 3T magnetic field, we used pure-silver, non-ferrous acupuncture needles (0.25mm diameter, 40mm length, Maeda Corporation, Japan), inserted to a depth of approximately 2–3cm. Electrical current was delivered with a modified current-constant H.A.N.S. (Han’s Acupoint Nerve Stimulator) LH202 (Neuroscience Research Center, Peking University, Beijing, China) for both VA and SPA. Current intensity was set by verbal subject response to be “moderately strong, but below the pain threshold” for both VA and SPA equally. Current intensity did not differ significantly (paired t-test, p>0.1) between VA (2.10 ± 0.96mA) and SPA (1.94 ± 0.89mA). Relatively strong stimulation and alternating frequencies were used to limit the potential for current intensity dropping below the sensory threshold over the lengthy scan run.
At the end of each fMRI scan run, subjects rated the intensity of sensations they felt during the run using a computerized i/o device. Subjects were presented with a 10 point visual analog scale (VAS) and were asked to rate the intensity of sensations commonly associated with the experience of deqi¥
, and listed on the MGH Acupuncture Sensation Scale (MASS) (Kong et al., 2007a
). The 10 point VAS used specific guidewords: “none” corresponding to a 0, and “unbearable” corresponding to a value of 10. Intermediate guidewords also included “mild,” which was linked to a value of 2, “moderate” which was linked to a value of 5, and “strong,” linked to a value of 8. The bar scale was presented and button box mediated responses were acquired (in increments of 0.5) with a laptop and Labview Software (Labview 7.1, DAQCard 6024E, National Instruments, Austin, Texas). A sensation was counted as present if the intensity was greater than 0.5. Subjects were also asked to assess the intensity of throbbing and sharp pain (this latter sensation is not considered to be characteristic of deqi
). Additionally, subjects were asked to assess the extent of “spreading” that may have occurred for any of the listed sensations. A modified version of this procedure has been successfully used by our group in the past to assess psychophysical response in conjunction with neuroimaging (Dhond et al., 2008
; Hui et al., 2005
; Hui et al., 2007
; Napadow et al., 2007
; Napadow et al., 2005b
). In order to quantify the overall, combined amplitude of deqi
experienced, we used the MASS-Index
, which is calculated by an exponentially decreasing weighted sum of all sensations (see (Kong et al., 2007a
) for details). We have developed this index as an attempt to balance the breadth and depth of sensations as well as the number of different sensations chosen by the subject. The MASS-index was compared between stimulation groups using a paired t-test, significant at p<0.05 (SPSS 10.0.7, Chicago, Illinois). Furthermore, frequency counts of specific sensations were compared between stimulation type with a Pearson Chi-squared test, significant at p<0.05.
Functional MRI analysis
Structural MRI data were co-registered to the MNI152 brain template using the Automated Brainstem Co-registration (ABC) method (Napadow et al., 2006
) in order to maximize sensitivity in important brainstem and subcortical brain regions. Functional MRI data were first motion corrected with AFNI (NIH) through an iterated, linearized, weighted least-squares method with Fourier interpolation (Cox, 1996
). Data runs were excluded if gross translational motion exceeded 3mm on any axis. Any residual motion was then removed by performing probabilistic independent component analysis (ICA) and removing components related to motion artifact (e.g. positive/negative fMRI response on opposing edges of the brain, IC timeseries spikes consistent with motion correction timeseries spikes) (Beckmann and Smith, 2004
We performed cardiac-gated fMRI in order to reduce cardiogenic motion artifact, an approach taken by several groups at our Center (Guimaraes et al., 1998
; Mainero et al., 2007
; Napadow et al., 2008
; Napadow et al., 2005a
; Zhang et al., 2006
). However, cardiac gating produces data with variable TR, and hence variable T1-weighting. If the temporal variability of T1-weighting is significant it will overwhelm any BOLD fMRI T2*-weighting. Thus, any MR signal fluctuation due to T1-variability must be accounted for before any model-based processing is applied. Data were corrected for T1 variability (due to variable TR) using the method detailed by Guimaraes et al. (Guimaraes et al., 1998
). Their study clearly demonstrated improved SNR (i.e. mitigated T1-variability) in brainstem regions (the inferior colliculus for auditory stimuli) compared to non-gated fMRI. As reported by Guimaraes et al., T1 was estimated by minimizing the variance in fitting an exponential T1 decay curve to fMRI signal intensity within each voxel. The signal intensity in any given voxel can be expressed as:
where the known variables include: Si,n
=signal in ith
time point, and tn
=TR at nth
time point. This expression is fit to find: Ai,n
=max signal in ith
voxel without T1 weighting, and T1i
=longitudinal relaxation time in ith
voxel. The algorithm then uses the estimated Ai,n
to correct Si,n
in each voxel to the average tn
across all time points (Tave
This procedure uses the fMRI data itself to perform the correction, thereby avoiding spatial registration errors which would be incurred if T1 were estimated in a separate MRI scan. While the amount of correction at each time point varies with the difference between the TR and Tave, and theoretically the TR could vary with variable heart rate between ON and OFF blocks, we found no difference between ON blocks and OFF blocks for both VA (ON:63.6±7.8bpm, OFF: 64.0±7.9, paired T-test: p=0.20) and SPA (ON:64.2±10.6, OFF: 64.2±10.9, paired T-test: p=0.92). Thus, the difference in mean R-R between ON and OFF blocks for VA was on the order of 5.9mS,
In addition, error in the estimated T1
may have been propagated to an error in BOLD signal correction. In order to estimate this propagated error we calculated the variability in T1
(standard deviation of T1
, or σ(T1
)) within several significant brainstem clusters. Using our collected data we assume an average TR to which all time points are corrected to be Tave
= 1900ms. Also using our data, we take the average correction time to be Δt = 81.5ms. Then, using our calculated mean T1
for each region, and assuming a T1
′ = T1
), as calculated for each region, we can compute the propagated percent error in the correction of signal intensity (%CorErr) as:
This error represents the effect of using T1′ instead of T1 in correcting a TR at Tave + Δt to Tave. We found the following mean T1, standard deviations, and propagated correction error for each of 3 representative ROIs: locus ceruleus (T1: 1152.6ms, σ(T1): 226.3ms, %CorErr: 8.92%), PAG (T1: 1261.9ms, σ(T1): 131.1ms, %CorErr: 4.06%), substantia nigra (T1: 1108.1ms, σ(T1): 101.5ms, %CorErr: 5.45%). Thus, the propagated error to the correction of a typical time point in these brainstem regions is roughly only 4–8% of the original correction. Our actual T1 correction for the ROI’s above (comparing corrected with uncorrected BOLD data) corresponded to an average percent signal change of 2% or less: locus ceruleus (1.8 ± 0.2%, μ±σ), PAG (2.0 ± 0.3%), substantia nigra (1.6 ± 0.3%). Thus we do not expect that a 4–8% error in this correction (resulting in a % signal change error of 0.08 – 0.16%), arising from a roughly 10–20% error in T1, would have a large impact on percent signal change calculations in our regions of interest.
It should also be noted that other viable options for physiological noise correction exist, and include cardiac-gating with T1-correction using a T1-map calculated in scans different from the original BOLD data (Malinen et al., 2006
), as well as non cardiac-gated retrospective algorithms using concurrently collected ECG data (Glover et al., 2000
Only minimal spatial smoothing, on the order of a voxel, (FWHM=3mm) was performed on the fMRI data, as some brainstem nuclei are on the order of image voxel size. The data were also high-pass filtered in the temporal domain (fhigh
= 0.0083 Hz) to remove baseline signal drifts. Corfield et al. found that high-pass filtering with this cutoff frequency yielded similar results to more intensive artifact mitigation procedures such as regressing out cardiac and respiratory phase timings in a general linear model (Corfield et al., 1999
Statistical parametric mapping at the single subject level was completed via a generalized linear model (GLM) by using FMRI Expert Analysis Tool (FEAT, FSL). Time-series statistical analysis was carried out using FILM (FMRIB’s Improved Linear Model) with local autocorrelation correction. The hemodynamic response function utilized in the GLM analysis was defined by the block design paradigm convolved with a prescribed gamma function (standard deviation = 3s, mean lag = 6s). The block paradigm was defined to be “1” during the 1-minute ON blocks and “0” during the 1-minute OFF blocks. As we wanted our analysis to be sensitive to time-variant signal increase and signal decrease response patterns, we did not use a single ON-OFF block design that encompassed the full run duration. Instead, each of the 15 ON blocks was modeled with a separate block step function regressor convolved with a canonical gamma HRF model. This approach allowed for analysis of the time varying response of the BOLD signal for each voxel, something not typically done in traditional fMRI block design analysis which has fewer blocks and much shorter scan run durations.
Main effect and time-variant group response was calculated at the second level by passing up each subject’s parameter estimate (and its variance) for each of the 15 regressors defined above. For the brainstem specific analysis, ABC-adjusted standard space parameter estimates (and their variance) from individual subjects were imported to a group-level analysis with FLAME (FMRIB’s Local Analysis of Mixed Effects, using a Markov Chain Monte Carlo (MCMC)-based mixed effects analysis). The same was done in MNI-space to infer activation/deactivation in cortical and subcortical structures rostral to the brainstem. As stimulation site was pseudo-randomized to be right or left, we performed two MNI-space analyses, one which flipped left-sided stimulation data across the mid-sagittal plane in order to infer activity in brain structures likely to respond contralateral to the stimulus side (thalamus, SI and SII, posterolateral parietal), and one which kept the data in its original right-left configuration (all other structures). For both VA and SPA, contrast parameters at the group level were specified which derived (1) the “main effect” across all 15 blocks from all subjects and (2) linear time-variant (decrease or increase) response. The latter was calculated by creating a contrast where the group response for each of the 15 blocks was weighted by a linearly decreasing factor, from 7 to −7 (zero mean). We chose not to set a single EV with linear scaled block-weighting at the first level because we wanted our analysis to differentiate between significant linear decreases manifest as strong activation in the beginning of the run trailing off to no activation at the end, versus a significant linear decrease where there is no activation in the beginning of the run but strong deactivation at the end. By estimating time-variant response as a group level contrast (collapsed over all subjects in the analysis for each block), and by plotting the time-variant GLM coefficients, the above two scenarios can be differentiated.
Both VA and SPA parameter estimates were also placed in a mixed-effects group level analysis in order to directly compare the two forms of stimulation with a VA - SPA contrast. For brainstem data, where activation clusters were expected to be small, resultant statistical parametric maps for the main effect and linear time-variant effect were cluster corrected (Gaussian Random Field theory) for multiple comparisons at p<0.01. For cortical and subcortical (but supra-brainstem) data, the main effect and linear time-variant group maps were threshold at a cluster corrected p<0.001. Average time courses () from distinct brain regions were created by first resampling single-subject fMRI data to an average TR of 2 seconds.
Figure 5 Verum acupuncture (VA) response in the substantia nigra [A] demonstrated linearly decreasing bimodal response over time - initial activation and later deactivation [B]. [C] There was a fundamentally different fMRI response to VA in this brain region during (more ...)
Classical habituation was defined in our linearly time-variant results by using the following criteria: (1) robust activation (positive percent signal change greater than 1 standard deviation from nil) for the first 3 blocks, (2) linear decrease in time, and (3) minimal response (percent signal change within 1 standard deviation of nil) for the last 3 blocks (see ).
Table 2 Linear time-variant response to acupuncture (VA and SPA). All listed structures demonstrated linearly decreasing (and not increasing) response. The percent signal change mean and standard deviation for the first 3 blocks and the last 3 blocks in the scan (more ...)
In order to test if the results of our psychophysical analysis correlated with inter-subject variability in fMRI response to the stimuli, we performed a mixed-effects ANCOVA, with the demeaned MASS-index from each subject as the covariate of interest. This analysis was performed separately for both VA and SPA, using group analysis methods identical to those above. For the MNI-analyses, data were corrected for multiple comparisons at a cluster corrected p<0.05 (reflecting the relatively lower degrees of freedom in this ANCOVA).
Anatomical localization for forebrain neuroimaging data was aided by a stereotaxic (Talairach coordinate system) atlas (Mai et al., 2004
). Given that the Talairach or MNI coordinate system is not specific for the brainstem, we adopted a previously validated approach (DaSilva et al., 2002
) under supervision of one of the investigators (NM), who is an experienced neuroanatomist. Briefly, localization was performed using a landmark-based topographical methodology wherein the human brainstem has been subdivided into 28 distinct regions of interest or parcellation units (PUs). These PUs included the rostral dorsal medulla (site of DMNX, NTS, and NAmb) and rostral dorsal pons (site of PBN, and locus ceruleus, LC), which were of principal interest in the present study. Activity within the PUs was assessed for consistency with known brainstem nuclei using an atlas (Paxinos and Huang, 1995
). In lieu of an MNI coordinate, we reference the location of brainstem activation clusters in an axial plane relative to a medullary landmark - the obex (located roughly at MNI z=−57mm on our group averaged anatomical volume, see ,). Furthermore, based on our estimates from postmortem human material and those reconstructed from anatomical atlases (Paxinos and Huang, 1995
), we can calculate the values for three representative structures - one small (locus ceruleus: approximately 27mm3
), one medium-sized (inferior colliculus: approximately 180mm3
), and one large (substantia nigra: approximately 570 mm3
). The rest of the structures are in-between these estimates. Given our voxel size, partial-volume effects would certainly be present for the smaller brainstem nuclei, and so we have used the phrase “consistent with” when attributing results to specific brainstem nuclei, reflecting this uncertainty.
Figure 2 [A] Group maps of the main effect for brainstem activation/deactivation elicited by verum (VA) and sham point (SPA) acupuncture. While VA induced both activation and deactivation in several brainstem nuclei, brainstem fMRI response at the same threshold (more ...)
Figure 3 Group map of linear time-variant response for verum acupuncture (VA) in brainstem, cortical, and subcortical brain structures. Red/yellow clusters represent fMRI response which linearly decreases from block to block, while blue/cyan clusters (none found) (more ...)