Twenty-four participants, experienced in the Tibetan Buddhist meditation practices of NDA and/or FA, participated in the experiments (16 male, eight female; median age, 50 years, range 32–62 years, all but one right-handed; Caucasian, 19, Asian four, Hispanic one). The data from two of the participants were discarded because they fell asleep during the experiment. Of the remaining 22 participants, 14 participated in both meditations, four participated only in NDA, and four only in FA. Participants were selected, on the basis of recommendations from their teachers, from various Buddhist meditation centers in the US. Their experience ranged from 8 to 33 years of meditation; this corresponds roughly to 4000–37,000 h of practice, (mean 13,727 ± 8920 STD), based on participant's reports of their daily practice, and assuming an average of 10 h per day for retreat practice (8–14 h/per day in most Tibetan Buddhist traditions). While we strived for homogeneity, and most of our practitioners (17), had experience only with the Tibetan Buddhism, seven had some experience with other contemplative traditions as well, such as Mindfulness (1), Vedanta (2), Taoism (3), and Zen (1). All participants had normal or corrected-to-normal vision, provided written informed consent, and were paid for their participation in the study. Experiments were conducted in accordance with the safety guidelines for fMRI research, as approved by the University Committee on Activities Involving Human Subjects at New York University.

Instructions to participants for “fixation” were: “Keep your gaze on the fixation point and allow your mind to wander naturally, and avoid engaging in any type of meditation.” For viewing of movie clips the instruction was the same, except without fixating. FA meditation is common to all schools of Buddhism and, in one form or another, to all other contemplative traditions. Thus, the instructions were standard: “Keep a steady focus on the object of meditation, the fixation point, to the exclusion of other mental content, and if thoughts arise do not follow them but remain focused on the fixation point. If your mind wanders off, bring it back gently to the object of focus.” For NDA meditation the instructions were: “Rest in reflexive NDA, equally aware of inside and outside of your body, allowing experiences to arise and subside of their own accord.” NDA meditators were also told to avoid common pitfalls in this style of meditation, such as, fixating attention either exclusively on the external or exclusively on the internal experience or, trying to keep one's mind empty of thoughts and feelings through holding effortful focus.

We did not have a reliable method of assessing the quality of meditation experience in the scanner or the effectiveness with which subjects were able to enter meditative states, beyond the most rudimentary, i.e., whether subjects fell asleep or not (eye-tracker), or whether they were able to meditate at all (self-report after each run). However, the possibility that their meditation was of a poor quality would not present a confound for our results, but rather would have diminished any differences in brain activity between conditions.

Stimuli were presented via an LCD projector and custom optics onto a rear-projection screen in the bore of the MRI scanner. Participants were supine and viewed the projected stimuli through an angled mirror. An infrared, video eye-tracker (ASL R-LRO-6, Bedford, MA, or EyeLink 1000, SR Research, Mississauga, ON, Canada) was used to visually monitor whether fixation was maintained and if the eyes were kept open, but the eye-position measurements were not recorded.

Each MRI scanning session consisted of four scans. Each of the first three were 480 s in duration, during which participants fixated a red dot located in the center of the screen on a neutral gray background. The participants either fixated (fixation), or engaging in NDA, or FA meditations (see Discussion for an explanation of differences between NDA, FA, and fixation). The order of these first three scans was counterbalanced across participants. During the fourth scan, participants viewed movie clips (10 s each), alternating with blank periods of neutral gray (8 s each). The movie clips showed faces, urban buildings, natural outdoor scenes, and various objects, and were presented in randomly shuffled order (Hasson et al., 2003). The experiment began and ended with neutral gray (20 s), and lasted for a total of 608 s per scan. For the meditation scans, participants were given time to enter into meditation before fMRI acquisition commenced. This way we avoided confounding the neural correlates of NDA or FA with the neural correlates of techniques to get to these states.

MR images were acquired at the NYU Center for Brain Imaging, using a Siemens Allegra 3T head-only scanner equipped with a head coil (NM-011, Nova Medical, Wakefield, MA). Blood oxygenation level-dependent (BOLD) contrast was obtained using a T2^*-sensitive echo planar imaging (EPI) pulse sequence (TR = 2000 ms, TE = 30 ms, 36 axial slices, 3 × 3 × 3 mm, 64 × 64 matrix in a 192 × 192 mm FOV). In addition to functional images, a high-resolution, T1-weighted anatomical volume of images was acquired for each participant, using an MPRAGE sequence (TR = 2500 ms, TE = 3.93 ms, TI = 900 ms, flip angle = 8^°, 176 sagittal slices, 1 × 1 × 1 mm, 256 × 256 matrix in a 256 mm FOV).

fMRI data were processed using standard tools available in the Brain-Voyager QX 2.2 software package (Brain Innovation, Maastricht, The Netherlands) (Goebel et al., 2006), with NeuroElf (http://neuroelf.net/, JW), with custom software written in Matlab (MathWorks, Natick, MA), and with Excel (Microsoft, Seattle, WA). The first four volumes of each fMRI scan were discarded, to allow the magnetization to reach steady state, and to allow participants time to recover from the onset of scanner noise. Functional MRI time-series were pre-processed with slice-time correction (to compensate for the slightly different times during which each slice was acquired), 3D motion correction (which involved re-sampling with windowed sinc interpolation), and temporal high-pass filtering (with a cutoff frequency of 0.006 Hz, removing trends with periods greater or equal to 150 s). The fMRI time-series were then corrected for serial correlation (Goebel et al., 2006), and temporally smoothed using a Gaussian kernel with a temporal extent of 8 s (four TRs, full width at half maximum). Anatomical images were corrected for spatial inhomogeneity in intensity, and transformed, using trilinear interpolation, to the Talairach coordinate system (Talairach and Tornoux, 1988). Functional images were aligned with the high-resolution anatomical volume using trilinear interpolation. Other than the interpolation during these transformations, no additional spatial smoothing was performed.

Artifacts due to non-neural signals such as head-motion, breathing, and heart rate were minimized as follows. The head motion variables were projected out (using linear regression) from the fMRI time-series. A pair of global components were extracted and removed: (1) the mean fMRI response time-series for all cortical gray matter voxels, and (2) the mean fMRI response time-series for all white matter and subcortical voxels. The cortical surface was reconstructed from the high-resolution anatomical images, separately for each participant, to identify the gray and white matter voxels. These global components were removed by projecting them out of the response time-series for each individual voxel (Golland et al., 2007; Fox et al., 2009).

Fifteen regions of interest (ROIs) were defined separately in each individual participant (Figure 1). Eight of these ROIs were part of the extrinsic system, exhibiting more activity during external stimulation, while seven were a part of the intrinsic system, exhibiting less activity during external stimulation (Fox et al., 2005; Damoiseaux et al., 2006; Fransson, 2006; Goldberg et al., 2006; Nir et al., 2006; Fox and Raichle, 2007; Golland et al., 2007, 2008; Castellanos et al., 2008; Vincent et al., 2008; Vanhaudenhuyse et al., 2010a). The ROIs were defined based on a combination of anatomical and functional criteria using the movie clip scan. Once the ROIs were defined, the main experimental measures (e.g., level of correlation/anti-correlation between these ROIs during NDA, FA, and fixation) were assessed using statistically independent fMRI data acquired in separate scans.

Brain regions were identified as responding either more strongly to the movie clips than neutral gray (all movie clips > neutral gray) or less strongly to the movie clips than neutral gray (all movie clips < neutral gray), following previous reports of segregating the extrinsic and intrinsic networks by examining brain activity evoked by watching movies (Hasson et al., 2003, 2010; Golland et al., 2007). In addition, area Lateral Occipital (LO) was defined by identifying brain regions that responded more strongly to movie clips of objects than all other movie clips (Hasson et al., 2003). We used Brain Voyager QX to run a standard statistical parametric mapping (SPM) analysis (Friston et al., 1994). In short, we constructed a general linear model (GLM) for the underlying neural responses to each movie clip. The regression matrix contained one row for each time point, in which neural activity was modeled as either “on” = 1 or “off” = 0, and five columns, one for each type of movie clip (visual motion, faces, objects, and houses) as well as a column of ones for the mean image intensity. The periods of neutral gray between movie clips served as the baseline. The columns of this regression matrix were convolved with a canonical hemodynamic impulse response function (HIRF) to create a model of the expected fMRI response time-series (Boynton et al., 1996; Glover, 1999). We used linear regression to estimate response amplitudes (beta values) for each voxel and each condition, solving an equation of the form y = Ax + ε, where the vector y was the measured fMRI time-series during one scan, the vector x contained the estimated response amplitudes, A was the regression matrix described above, and ε was the residual variance, orthogonal to the model A. Finally, voxels exhibiting statistically significant response differences were identified (Figure 1A; p < 0.05; t-test, not corrected for multiple comparisons across voxels, but with a minimum cluster volume of 162 mm³).

A set of anatomical landmarks were identified by hand, one for each ROI in each participant's brain, and the functional activation in the vicinity of these landmarks was used to define the ROIs. Specifically, we identified the voxel within 8 mm of each landmark, separately for each individual participant, which exhibited the strongest functional activation (Figure 1B). This voxel was selected as the ROI center for that ROI in that participant. Each ROI was then created as a cube (12 mm on each side) surrounding the ROI center. Finally, each ROI was restricted to include only voxels exhibiting statistically significant (although not corrected for multiple comparisons) response differences (Figure 1C). Hence, the ROI centers and sizes differed across participants (Table 1).

Eight extrinsic ROIs (four in each hemisphere) were defined (Table 1): “Vis,” an early visual area located in and around calcarine sulcus; LO, area in LO and inferior temporal sulcus; V3+, in and around middle occipital gyrus; SPL, superior parietal lobule in and around the superior intra-parietal sulcus. Each ROI was defined separately for each hemisphere.

Seven intrinsic ROIs were defined: (Table 1): PCC, in and around the rostral portion of posterior cingulate gyrus; right and left LP, lateral parietal area in the inferior parietal lobule defined separately in each hemisphere; Rsfg, right superior frontal gyrus (BA8/6); mPFC, perigenu area of mPFC; cPrec, in and around ventral portion of subparietal sulcus and posterior to precuneal sulcus; dPrec, in and around dorsal caudal portion of subparietal sulcus. There were two peaks of de-activation in the precuneus, which is why we defined two separate precuneus ROIs (cPrec and dPrec) for each participant. This division reflected the previously observed functional and anatomical differentiation in precuneus (Immordino-Yang et al., 2009; Margulies et al., 2009). LP was defined separately in each hemisphere. rSFG was identified only in the right hemisphere. The other intrinsic ROIs were medial, and in close anatomical proximity in the two hemispheres; consequently, we picked a single ROI center and defined a single ROI around that center location that included voxels from both hemispheres.

The main analysis computed inter-area correlations (also called “functional connectivity”) between the time-series of brain activity in different ROIs (Friston et al., 1996). We extracted the time-series of each voxel in each ROI and averaged across voxels to give an average time-series for each ROI, each participant, and each condition (fixation, NDA, FA). Each such ROI time-series was then z-scored by subtracting its mean and dividing by its variance. We obtained the grand average time-series for the extrinsic and intrinsic systems by averaging the z-scored time-series of all ROIs belonging to each of the two systems. We then computed the correlation between those two grand average time-series as an overall measure of the level of anti-correlation between the two systems for each of the three conditions (Figures 3 and 4A).

In addition, we computed the inter-area correlation separately for each pair of ROIs, and averaged the Fisher-transformed correlation values across pairs, and across participants. This was done separately for each pair of extrinsic ROIs, each pair of intrinsic ROIs, and for each extrinsic-intrinsic pair (Figure 4B).

Statistical analyses tested for differences in inter-area correlations between conditions (e.g., NDA vs. fixation) and across participants. The correlations were computed separately for each individual participant and for each condition (fixation, NDA, FA), using individually defined ROIs. Paired t-tests were then used to assess differences in the correlation values between conditions (e.g., NDA vs. fixation). Because correlations are not normally distributed, the t-tests were computed on Fisher-transformed correlation values.

In addition, to rule out the influence of temporal autocorrelation on the time-series, which might have inflated the correlation statistic, particularly as a smoothing kernel of four TRs had been applied, we also used non-parametric randomization tests to assess the statistical significance of the inter-area correlation differences (Hasson et al., 2008; Lerner et al., 2011; Wang et al., 2010). Specifically, we took the discrete Fourier transform of the time-series in each ROI, randomly permuted its Fourier phase without changing the amplitude, inverted the Fourier transform, and recomputed the correlation between the resulting time-series. This was done for each condition (fixation, NDA, FA), and we computed the differences in the resulting correlation values between conditions (e.g., NDA vs. fixation). This procedure was repeated 1000 times to yield a null distribution for the correlation-differences. Finally, we determined a p-value as the fraction of the null distribution (regardless of sign so as to be a two-tailed test) that was as large or larger than the correlation-difference observed without randomization (Figure 2).

We computed the standard deviation of the fMRI response time-series in each ROI to determine if there was any difference in the amplitude of modulation of the fMRI responses between conditions. If there was any such a difference, it might reflect differences in the neural activity between conditions or it might reflect confounding factors such as differences in eye movements, head movements, pulse, or respiration. Specifically, we extracted the time-series of each voxel in each ROI and averaged across voxels to give an average time-series for each ROI, each participant, and each condition (fixation, NDA, FA). Each such time-series was divided by its mean to convert from arbitrary units to percent modulation. We averaged the time-series across voxels for each ROI and computed the standard deviation of the resulting mean time-series. Finally, these standard deviations were then averaged across ROIs and across participants, separately for each condition (fixation, NDA, FA), and separately for the intrinsic ROIs and the extrinsic ROIs (Figure 5).

We used non-parametric randomization tests to assess the statistical significance of any differences in the standard deviations. Specifically, we randomly permuted the condition labels (fixation, NDA, FA) for each ROI separately, and we recomputed the standard deviations. This procedure was repeated 1000 times to yield null distributions for the standard deviations. Finally, we determined a p-value as the fraction of the null distribution (regardless of sign so as to be a two-tailed test) that was as large or larger than the observed value without randomization.

The inter-area correlation between the extrinsic and intrinsic systems depended on whether participants fixated, or engaged in NDA or FA meditations (Figures 3 and 4). We computed the grand-averaged time-series for all voxels in the extrinsic and intrinsic ROIs (see Figure 1 and Methods for details), and then computed the correlation between those two grand average time-series as an overall measure of the level of anti-correlation between the two systems (see Figure 3 for an individual participant example and Figure 4A for the correlations from all participants). Brain activity was significantly more anti-correlated during FA than fixation, and significantly less anti-correlated during NDA than fixation (p < 0.001, paired t-tests; p < 0.005, phase-randomization test, see Figure 2 and Methods for details).

These differences in anti-correlation were specific to interactions between the extrinsic and intrinsic systems, as there was no evidence for differences in correlations between pairs of intrinsic ROIs nor between pairs of extrinsic ROIs (Figure 4B and Table 2). We computed the inter-area correlation separately for each pair of ROIs. This was done separately for each pair of extrinsic ROIs, each pair of intrinsic ROIs, and for each extrinsic-intrinsic pair. We then averaged the correlation values across pairs of ROIs, and across participants (see Methods for details). Averaged pairwise correlations for extrinsic-intrinsic pairs were negative and differed significantly between conditions. For pairs of extrinsic system ROIs, on the other hand, the correlations were positive and did not significantly differ between conditions. The correlations were also positive for pairs of intrinsic ROIs; there was a marginally significant difference between FA and NDA but not between the other pairs of conditions.

There was no evidence for any differences between conditions in the modulation of brain activity, neither in the extrinsic system ROIs nor in the intrinsic ROIs (Figure 5). Because correlation depends on both covariance and variance, the observed differences in the correlations (see above) might have been due to changes in covariance or variance or both. We computed the standard deviations of fMRI response time-series for each ROI in each participant and each condition, then averaged them across ROIs separately for extrinsic ROIs and intrinsic ROIs, and separately for each of the three conditions (fixation, FA, NDA) (see Methods for details). Standard deviations were higher in the intrinsic network than in the extrinsic one, across all three conditions (p < 0.005, randomization test). However, the standard deviations were not significantly different across the three conditions (fixation vs. FA: p = 0.9880 extrinsic, p = 0.4560 intrinsic, randomization test; fixation vs. NDA: p = 0.4540 extrinsic, p = 0.2920 intrinsic, randomization test; FA vs. NDA: p = 0.7440 extrinsic, p = 0.3580 intrinsic, randomization test). Hence, the differences in the correlations between conditions reflected differences in the interactions between extrinsic and intrinsic system ROIs, not in the modulation of activity within these brain areas.

We wanted to ensure that the observed differences in anti-correlation between extrinsic and intrinsic systems across meditation styles were not due to differences in the variance of the global component that had been regressed out. If the global component had large amplitude of modulation during one condition than for another, it might have artificially generated stronger anti-correlations. We found, however, that there was no difference in the variance of the global components across conditions (p >0.25, two-tailed t-test and randomization test). Furthermore, the variance of the global component was largest (though not statistically significantly larger) in the NDA condition where anti-correlation between extrinsic and intrinsic was weakest, suggesting that global component variability did not contribute to the anti-correlations.

A follow-up study would test naïve participants, provided with a minimal training, to determine if our results are the effect of long-term learning/practice, or if similar results could be produces in a shorter time as well.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.