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 (Figure 1a). 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 (Figure 1a). 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.

BOLD functional imaging was performed using a gradient echo T2*-weighted pulse sequence (TE = 30msec, matrix = 64 × 64, FOV = 200mm, flip angle = 90°) using a GRAPPA (Griswold et al., 2002) acceleration factor of 2 with our multi-channel coil. Volume acquisition was gated to the subjects’ every other ECG R-wave; thus, TR varied with each acquisition (average TR was approximately 2sec). Cluster acquisition of 17 coronal slices was completed in sequential order within 926ms. This was the maximum number of slices we could acquire (with our coil and scanner setup) within a time period which was approximately a single cardiac R-R interval. Each slice was 3mm thick with 0.6mm gap (voxel size 3.13 × 3.13 × 3.6mm) and was slightly tilted to be parallel to the brainstem longitudinal axis (Figure 1b). This orientation was chosen in order to minimize through-slice brainstem motion. Motion correction was incorporated into the fMRI scanning sequence with Prospective Acquisition Correction (PACE, Siemens Medical Systems, Erlangen, Germany) – a critical addition due to increased opportunity for head motion during our long duration scan run. More details of our cardiac-gated fMRI procedures have been previously published (Napadow et al., 2008).

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: S_i,n=signal in i^th voxel, n^th time point, and t_n=TR at n^th time point. This expression is fit to find: A_i,n=max signal in i^th voxel without T1 weighting, and T1_i=longitudinal relaxation time in i^th voxel. The algorithm then uses the estimated A_i,n and T1_i to correct S_i,n in each voxel to the average t_n across all time points (T_ave).

This procedure uses the fMRI data itself to perform the correction, thereby avoiding spatial registration errors which would be incurred if T₁ were estimated in a separate MRI scan. While the amount of correction at each time point varies with the difference between the TR and T_ave, 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 T₁ may have been propagated to an error in BOLD signal correction. In order to estimate this propagated error we calculated the variability in T₁ (standard deviation of T₁, or σ(T₁)) within several significant brainstem clusters. Using our collected data we assume an average TR to which all time points are corrected to be T_ave = 1900ms. Also using our data, we take the average correction time to be Δt = 81.5ms. Then, using our calculated mean T₁ for each region, and assuming a T₁′ = T₁+σ(T₁), 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 T₁′ instead of T₁ in correcting a TR at T_ave + Δt to T_ave. We found the following mean T₁, standard deviations, and propagated correction error for each of 3 representative ROIs: locus ceruleus (T₁: 1152.6ms, σ(T₁): 226.3ms, %CorErr: 8.92%), PAG (T₁: 1261.9ms, σ(T₁): 131.1ms, %CorErr: 4.06%), substantia nigra (T₁: 1108.1ms, σ(T₁): 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 (f_high = 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 (Figure 5) from distinct brain regions were created by first resampling single-subject fMRI data to an average TR of 2 seconds.

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).

The hippocampus also demonstrated linearly time-variant activation/deactivation to VA. This limbic structure is crucial for memory consolidation and, as a component of the default-mode network (Raichle et al., 2001), is deactivated by acupuncture (Hui et al., 2000; Hui et al., 2005; Napadow et al., 2005b) as well as other externally directed stimuli (Buckner et al., 2008). The time-variant response found in this study suggests that activity in this area changes over time during acupuncture stimulation and is consistent with our previous finding that connectivity between this area and the default-mode network is enhanced not just during but also in a resting state after acupuncture stimulation (Dhond et al., 2008).

Studies in animal models have shown that acupuncture modulates DA neurotransmission. For example, DA down-regulation has been noted in several animal (Han et al., 1999; Wang et al., 1999). Our results in the human clearly show that VA modulates activity in the VTA and SN, a source of dopaminergic tone to higher brain structures (Figures 2, ​,3).3). VA induced activation in the SN, and a slightly different SN region was found to have activation to early VA stimuli, and deactivation to late (>25min) stimuli (Figure 5). In comparison, acupuncture-induced SN deactivation in previous studies (Hui et al., 2005) may have been due to differences in brain response to electrical versus manual needle stimulation. However, we did find SN deactivation in latter blocks, which is also consistent with DA downregulation in animal models following longer duration (40 minutes) EA at the same acupoint on the leg, ST-36 (Wang et al., 1999). In general, difficulties in comparing our human neuroimaging results with those using invasive methods in animal models are significant. While similarities exist, differences could be easily attributed to species differences and awake versus anesthetized states; both of which lead to variable influences of cognitive, top-down modulation on lower brainstem regions.

Acupuncture modulation of serotonin (5-HT) has also been investigated with animal models and this system is important for both disorders of pain and affect. Electroacupuncture appears to increase 5-HT synthesis and utilization in the rat (Han et al., 1979), specifically in the DRN and NRM (Kwon et al., 2000). Our results in the human demonstrated NRM deactivation for VA.

Noradrenaline (NA) modulation by acupuncture was supported by LC deactivation in response to VA. The LC is the source of noradrenergic input to higher forebrain structures and is associated with arousal and the hypothalamic-pituitary-axis stress response, relaying afferent signaling to the hypothalamus and amygdala (Ziegler et al., 1999). Thus, deactivation in the human LC is consistent with decreased stress response during VA.

We would like to thank the National Center for Complementary & Alternative Medicine, NIH for funding support: K01-AT002166 (to VN), P01-AT002048 (VN, KKSH), K01-AT004481 (RD), F05-AT003770 (KP) and K24-AT004095. We also acknowledge the NCRR (P41-RR14075, GCRC M01-RR01066), and the Mental Illness and Neuroscience Discovery (MIND) Institute. The content is solely the responsibility of the authors and does not necessarily represent the official views of our funding agencies.