Twenty-one long term practitioners of MM and 23 meditation-naïve control participants underwent fMRI scanning as part of a previous study (Berkovich-Ohana et al., 2016a), including several participants from a pilot study that was run with the same settings. All participants were Caucasian, right-handed, having University level education, and with no history of neurological or psychiatric disorders. Four long term meditators and four control participants were excluded from the current study due to missing data, leaving 17 meditators (MM, mean ± SD age 41.9 ± 9.3 y.o., 6 females; see Table 1) and 19 control participants (C, mean ± SD age 40.3 ± 10.98 y.o., 7 females; see Table 1). We have considered the size of our sample as sufficiently reliable, based on previous studies, which showed that for a liberal threshold of 0.05, about 12 subjects were required to achieve 80% power at the single voxel level for typical activations (Desmond and Glover, 2002; see also Friston, 2013 for a critical view on large sample sizes and over-powered studies). In addition, we also calculated the effect size for each instance of significant between-groups effect (explained further on).

The MM group consisted of long-term practitioners (mean ± SD 8442 ± 8430.42 h of formal practice, range 936–29, 293 h of formal practice over mean ± SD 15.2 ± 5.22 years; see Table 1). Hours of formal practice calculations included both accumulated hours in formal retreats, as well as accumulated daily practice throughout the years. All the meditators practice mindfulness meditation according to the Satipatthana and Theravada Vipassana traditions and were mostly recruited via the Israeli Insight Society TOVANA. All participants provided written informed consent for their participation. The experimental procedures were approved by the Tel Aviv Sourasky Medical Center - Ichilov Hospital ethics committee. This study was aligned with the Declaration of Helsinki and performed in accordance with the relevant guidelines.

Before entering the scanner, participants completed Hood’s Mysticism Scale (Hood, 1975). MRI images were acquired on a 3 T Trio Magnetom Siemens scanner, at the Weizmann Institute of Science, Rehovot, Israel. Participants underwent a T1-weighted 3D Magnetization Prepared - RApid Gradient Echo (MP-RAGE) sequence, which was used for anatomical segmentation and template normalization (1x1x1mm3 resolution, TR = 2,300 ms, TE = 2.98 ms, TI = 900 ms, and FA = 9°), followed by resting-state fMRI (T2*-weighted, obtained with gradient echo planar imaging sequence), which was acquired continuously for 9 min (TR = 3,000 ms, TE = 30 ms, FA = 90°, inplane resolution—3x3mm², slice thickness—3 mm, 46 axial slices; see also Berkovich-Ohana et al., 2016b; Berkovich-Ohana et al., 2016a). This scanning sequence is sensitive, but not exclusively, to the BOLD signal response (Chen et al., 2020). The participant’s head was placed on a foam cushion for stabilization, and MR compatible earphones (MR confon, Magdeburg, Germany) that substantially reduce external noise were placed on the ears. Participants were instructed to rest with their eyes closed, without falling asleep. In order to rule out cases of sleep or meditation practice during the resting-state, each participant was interviewed while in the scanner immediately after the conclusion of the scanning session via an intercom, using a semi structured interview. The participants were asked to describe their experience and thought content—a procedure which lasted around 5 min (Berkovich-Ohana et al., 2016a). In addition, each participant performed a visual recognition memory task (Hasson et al., 2003), the results of which have been previously published (Berkovich-Ohana et al., 2016b; Berkovich-Ohana et al., 2016a).

The anatomical MRI and resting state fMRI scans of the participants were acquired in DICOM format and converted into NIFTI files in the BIDS format using dicm2nii conversion tools¹ in MATLAB. The fMRI data in this manuscript were preprocessed using fMRIPrep 20.2.1 with default settings (Esteban et al., 2019; Esteban et al., 2022a), which is based on Nipype 1.5.1 (Gorgolewski et al., 2011; Esteban et al., 2022b). Anatomical data preprocessing included the following steps: correction for intensity non-uniformity, skull stripping, brain tissue segmentation, brain surface reconstruction, and normalization to standard Montreal Neurological Institute (MNI) space. Preprocessing of the resting state fMRI T2* scans included correction of head motion, slice-time correction, co-registration to structural image and non-linear spatial normalization to MNI space. Many internal operations of fMRIPrep use Nilearn 0.6 (Abraham et al., 2014), mostly within the fMRI processing workflow (see the Supplementary methods section in the Supplementary materials for additional details).

The data used here was only minimally preprocessed, without special care to remove non-neuronal components, and therefore no assumptions are made regarding the origins of the fluctuations in the fMRI signal (Chen et al., 2020). The aim of the current study is to examine potential differences in functional connectivity dynamics between veteran practitioners of mindfulness meditation and meditation-naïve individuals, without separating between neuronal and non-neuronal sources.

The fMRI signals obtained at the voxel level were parcellated into N = 360 cortical parcels using the Cole-Anticevic Brain Network Parcellation (CAB-NP), allowing this work to benefit from the higher precision of the multimodal Glasser parcellation scheme (Glasser et al., 2016), while also associating each parcel with an established resting-state functional network (Ji et al., 2019). Since our data lacks the T2-weighted image required for the Human Connectome Project (HCP) pipelines, the relevant images were first converted from T1-weighted space into CIFTI format using Ciftify toolbox version 2.0.3 (Dickie et al., 2019), which uses a surface-based pipeline similar to that of the HCP (Glasser et al., 2013). The resulting parcellated data were stored N × T matrices in MATLAB-compatible files for the next stages of the analysis, where N = 360 is the number of cortical parcels and T = 180 is the number of volumes in each scan.

To detect recurrent states of functional connectivity (FC), first the mean time-series of each of the N = 360 parcels was transformed into an analytic signal (represented, at each instant of time, by a phase θ n t and an amplitude A n t , where n ∈ 1 … 360 and t ∈ 2 , … , T − 1 ) using the Hilbert transform. The first and last volumes were excluded due to possible signal distortions induced by the Hilbert transform (Vohryzek et al., 2020). The dynamic phase coherence (dPC) was computed by estimating the phase alignment between pairs of brain regions n and p at each timepoint t, using the cosine similarity equation:

The dynamic phase coherence dPC is a tridimensional tensor with size N × N × T ′ obtained for each scan, where N is the number of brain parcels and T ′ = 178 is the number of time points in each fMRI scan. The dPC values range between −1 and 1, where 1 corresponds to the cases where the signals in brain regions n and p , at timepoint t have the same phase, d P C n p t = cos 0 ° = 1 , whereas −1 corresponds to a phase difference of half a cycle (or 180°).

The phase relationships at each instant of time can be captured with reduced dimensionality and increased sensitivity using the Leading Eigenvector Dynamics Analysis (LEiDA) algorithm (Cabral et al., 2017). LEiDA consists of extracting the leading eigenvector, V 1 n , of the dPC matrices at each timepoint t. The leading eigenvector V 1 is a vector with size N x 1 that captures the dominant pattern of phase relationships between brain regions. Each of the N elements of the vector V 1 n represents how the signal in each brain region n is aligned in phase with respect to the other brain regions. Elements with the same sign (positive or negative) in the leading eigenvector indicate regions that are phase-aligned and can be considered part of the same community. The magnitude of the elements V 1 n indicates the extent to which the respective brain region n belongs to its community. Given that the relative sign in eigenvectors is arbitrary, a convention is applied such that the largest community has a negative sign, while the smallest community has a positive sign (Eraifej et al., 2023; Farinha et al., 2022; Hancock et al., 2022).

Recurrent states of phase coherence were detected by applying a k-means clustering algorithm on all 6,408 leading eigenvectors obtained from the dPC matrices of all 36 participants. The k-means algorithm serves to detect recurrent patterns in the data, but requires determining a priori the number of clusters k. In the resting state literature, the number of functional networks is unclear, and the current consensus is that there is some hierarchy, which means that increasing the level of granularity might reveal more specific functional networks. The aim of this work is to attempt the detection of connectivity patterns, the occurrence of which differs between meditators and controls. To perform such an exploratory analysis, the k-means algorithm was run with k ranging from 2 to 20 and using the cosine distance as a measure of similarity between eigenvectors. Focusing on connectivity patterns that are not exclusive to a specific partitioning model, we opted not to rely on the metrics that are usually used to determine the optimal k, such as the Elbow method or Silhouette coefficient. For each value of k, this algorithm partitions the eigenvectors into k clusters. Averaging all the eigenvectors assigned to a given cluster returns the cluster centroid, V c , with size Nx1, which is considered to represent a state of phase-locking, or PL state, in brain signals. This term is distinct from the phase locking value (PLV), which is used in other studies (e.g., Ponce-Alvarez et al., 2015). In essence, the centroids of the various PL states capture modes of phase synchronization across brain regions, where elements with the same sign represent in-phase locking, and anti-phase locking is represented by elements with opposite signs (Hancock et al., 2022). It has been consistently found that the smallest subset of regions with the same sign in these cluster centroids exhibit a significant spatial overlap with known intrinsic connectivity networks and therefore this smaller subset of brain regions is considered as representative of the functional subsystems activated in each PL state (Eraifej et al., 2023; Larabi et al., 2020; Lord et al., 2019).

The clustering assigns to each timepoint a single PL state, making it possible to characterize each fMRI scan by its sequence of states over time. Using the state time course, we calculated the probability of occurrence of each state, which is simply the number of timepoints assigned to a given PL state divided by the total number of timepoints in each scanning session. For each partition model (i.e., for different numbers of clusters k) the probabilities of each PL state were calculated for each subject.

Differences in probabilities of occurrence between the groups were statistically assessed using a permutation-based non-parametric two sample hypothesis test, which estimates the null distribution by permutations of group labels, instead of relying on the test-type standard distributions. For each of 10,000 permutations, a t-test was applied to compare populations, and a p-value is returned. Correction for multiple comparisons was performed for each partitioning model separately using the conservative Bonferroni correction, i.e., by dividing the threshold for significance by the number of states in the model.

As an additional means of validating the robustness of our results, the effect size of each significant between-groups difference was estimated (Farinha et al., 2022). This was done using the biased form of Hedges’ g statistic, also known as Cohen’s d, which measures the effect size for differences between means (Hedges, 1981), and calculated as

where s ∗ is the pooled standard deviation, which is computed as

Having found PL states that exhibit significant between-groups differences in terms of their probability of occurrence in their respective partitioning models, these PL states were grouped into coupling modes (CMs). Each CM represents a recurring mode of cortical phase synchronization that can be found in several partitioning models, and therefore can be considered independent of the granularity of the model. This grouping was performed on the basis of visual inspection and validated by computing correlations for each pair of relevant PL states between their respective centroids, which stand for the connectivity pattern they represent.

A core mechanism of mindfulness meditation is changes to self-awareness (Tang et al., 2015), mainly studied via alterations to the DMN and largely thought to support the sense of self (Northoff et al., 2006). Importantly, the same cohorts reported in this study were previously shown to exhibit significant differences in various nodes of the DMN, both on the functional level (Berkovich-Ohana et al., 2016a), as well as the structural level (Berkovich-Ohana et al., 2020). In order to assess any potential link between functional connectivity dynamics and changes in self-awareness, all participants answered Hood’s Mysticism Scale (Hood, 1975), which is considered to be a measure of ego-dissolution, which is an alteration to the sense of self experienced during deep meditative states (Millière et al., 2018). This questionnaire is composed of eight subscales that relate to different qualities of mystical experiences: Ego quality refers to experiencing the loss of the sense of self while maintaining consciousness (e.g., “I have had an experience in which something greater than myself seemed to absorb me”); Unifying quality refers to the perception of unity among separated objects (e.g., “I have had an experience in which I realized the oneness of myself with all things”); Inner Subjective quality refers to perceiving an inner subjectivity even to material objects (e.g., “I have had an experience in which all things seemed to be conscious”); Temporal/Spatial quality refers to the modification of the temporal and spatial parameters of the experience (e.g., “I have had an experience which was both timeless and spaceless”); Noetic quality refers to treating the experience itself as a valid source of information (e.g., “I have had an experience in which a new view of reality was revealed to me”); Ineffability refers to the inability to verbally express the nature of the experience (e.g., “I have had an experience that cannot be expressed in words”); Positive Affect refers to the experience of joy or blissful happiness (e.g., “I have experienced profound joy”); and Religious quality, which refers to the sacredness of the experience (e.g., “I have had an experience which I knew to be sacred”). Correlations were calculated between the probabilities of occurrence for each PL state and participants’ scores on the various subscales of Hood’s Mysticism Scale, along with the mean score across the different subscales, both for the whole sample and for each group separately.

Considering that expertise in meditation develops across prolonged timescales, spanning years of practice and training, correlations were also calculated between the probabilities of occurrence for the different PL states and participants’ assessed meditation expertise (by years and by hours) only within the MM group. Important factors that may contribute to, or confound, findings of interest, relate to the age of the participants (age-related changes in the neural profiles of practitioners have been shown to differ from those among meditation-naive individuals; Luders, 2014), and the amount of time of training that they experienced throughout the years (e.g., Berkovich-Ohana et al., 2020). Therefore, the probabilities of occurrence of the various PL states were also correlated with participants’ age, aiming to separate between changes in functional connectivity dynamics that relate to age from those that relate to the accumulation of experience in the practice of meditation. Similarly to the correlations with the scores of Hood’s Mysticism Scale, the correlations between the probabilities of occurrence and age were computed both for the whole sample and for each group separately, while the correlations with meditation experience were computed only for the meditators group. All correlations were corrected for multiple comparisons for each partitioning model separately, using Holm-Bonferroni corrected p-values.²

The participants were explicitly requested not to use meditation techniques during the resting-state scans. In order to rule out the possibility that they did, as well as to collect incidents of drowsiness, all participants were interviewed immediately after the scan about their experience during it (Berkovich-Ohana et al., 2016a). None of the participants reported drowsiness, and more importantly, none of the MM practitioners reported entering a meditative state, such as anchoring on the breath or body sensations, which are the typical instructions in MM practice (Gunaratana, 2002).

In line with previous applications of LEiDA (Alonso Martínez et al., 2020; Cabral et al., 2017; Farinha et al., 2022; Lord et al., 2019; Vohryzek et al., 2020, 2022), our analysis revealed a repertoire of phase-locking (PL) states in the resting-state fMRI signals that reoccurred over time and across participants. These PL states, obtained by clustering the instantaneous leading eigenvectors into a different number of clusters in each partitioning model (ranging from k = 2 to k = 20), can be characterized as dividing all cortical parcels to two communities according to their instantaneous phases. The smaller community is considered to be in anti-phase with the rest of the brain, and thus only the cortical parcels that belong to that community are rendered in the representation of the PL states (Figure 1B), with each rendered parcel colored according to the resting-state network to which it was exclusively assigned to in Ji et al. (2019), as detailed in Figure 1A. Each row shows the PL states obtained by increasing the number of clusters from 2 to 20, representing a separate partitioning model. Within each row, PL states are sorted according to the likelihood of occurring differently (i.e., rejecting the null hypothesis) when comparing their probability of occurrence between meditators and controls, with this likelihood increasing from left to right. Aligning with previous applications of LEiDA, each partitioning model contains one PL state referred to as the global mode, in which all parcels are relatively aligned in phase. These PL states are rendered in Figure 1B without the representation of any cortical parcels, as none are in anti-phase with the rest of the brain. These PL states have the highest probability of occurrence in their respective partitioning model (see Supplementary Figure S1), aligning the current study with other applications of the method. Having used a purely cortical and finer grained parcellation scheme, the patterns found here closely overlap with those obtained in fMRI data from the Human Connectome Project (Hancock et al., 2022; Vohryzek et al., 2020).

The probability of occurrence of each of the phase-locking (PL) states (Figure 1B) was compared between the long-term meditators (MM) and the control (C) participants. The clusters returning the most significant differences between groups are marked by highlighting the resulting p-values above the corresponding brain network images in either red or green, corresponding to p < 0.05 and corrected p < 0.05/k, respectively. Since the clusters are not independent for different k partitions, we grouped the clusters that showed significant differences into different coupling modes (CM), according to the correlations between them (Figure 2; see Supplementary Figure S2 for a visual comparison of the centroids of the individual PL states). Since the scores and probabilities are not normally distributed, Spearman’s rank correlation coefficient ρ was used in the entire study. The first coupling mode consists of PL states that were detected more often in meditators than in control participants (Table 2 CM 1; Figure 3A). These PL states mostly include regions that belong to primary and secondary visual networks, alongside several dorsal brain regions assigned to the dorsal attention network, DMN, and FPN (Figure 3B). The PL states in the second CM were also detected more often in meditators compared to controls (Table 2, CM 2; Figure 3C), highlighting a small group of regions from the secondary visual network (Figure 3D).

In contrast, the last CM of PL states was detected more often in control participants compared to meditators (Table 2, CM 3; Figure 3E). While the PL states in this CM did not survive the correction for multiple comparisons, they are quite consistent across partition models, involving coordination between latero-frontal and parietal components of the FPN and specific regions from the cingulo-opercular and default networks (Figure 3F). An additional pattern was found to differ between groups when clustering into k = 19 clusters (Table 2, Partition-specific pattern; Supplementary Figure S3), implicating brain areas distributed among ventral, dorsal, frontal, and parietal areas. However, given that comparing a higher number of states increases the chances of false positives and the difference did not survive correction, we did not consider this result in the subsequent analysis.

While the probability of occurrence of the global mode in all partitioning models did not significantly differ between the groups, meditators exhibited a slight preference of these PL states regardless of k (Supplementary Figure S4).

The probabilities of occurrence of certain PL patterns were found to be correlated with either age or scores in the Hood’s Mysticism Scale (Table 3), but not with meditation expertise (among meditators). For the groups’ comparison in the Hood’s Mysticism Scale scores see Supplementary Table S1. Among all participants, irrespective of meditative experience or partitioning model, synchronization between certain parts of the anterior cingulate cortex (ACC), medial prefrontal cortex, and orbital or polar frontal cortex, was consistently detected more often in older participants (Table 3, Whole Sample; Figure 4A). The probability of occurrence of the global mode was negatively correlated with age among meditators (Supplementary Table S2), while a trend for positive correlations was found between age and probability of occurrence of the global state among the control participants in all partitioning models. However, these correlations achieved significance only in a subset of the partitioning models for the meditators group (Table 3, Meditators; Figure 4C), but not in any partitioning model for the control group. Moreover, no significant correlations were found for the control group between the probability of occurrence of any PL state, in any partitioning model, and their scores in any of the subscales of Hood’s Mysticism Scale.

Two of the PL patterns exhibited significant correlations between their respective probability of occurrence and scores of the Mysticism scale. Both patterns encompass regions that are both spatially distributed in the brain and are associated with a variety of documented functional networks. Across the entire sample, scores in the Unifying quality aspect of the scale were found to be positively correlated with the probability of occurrence of a PL pattern composed of posterior regions that mostly belong to the visual, frontoparietal, dorsal attention, and default networks (Table 3, Whole Sample; Figure 4B). Meditators exhibited a negative correlation between their scores for ‘Noetic Quality’ and the probability of occurrence for a state in which synchronization occurs between orbitofrontal regions from the DMN, alongside regions from the somatomotor and secondary visual networks (Table 3, Meditators; Figure 4D).

The current study bears several limitations. First, as with any cross-sectional study, it is difficult to disentangle the effects of meditation and those of other lifestyle parameters. Second, the use of the extensive Glasser atlas in the present study meant that each observation (i.e., each time point for each participant) contains many features, which renders the study underpowered and warrants a larger participant pool than available. This also led to the decision to limit the focus of the analyses on the cortical parcels defined by Glasser et al. (2016), without including the sub-cortical expansion performed by Ji et al. (2019). Future studies could produce more robust results by using a larger sample size and should implement a longitudinal design to better distinguish between the effects of mindfulness meditation practice from those of other life experiences. In addition, examining FC dynamics during the practice of various forms of meditation could provide more detailed information regarding the effects of these practices. Another possibility would be to measure the performance of experienced mindfulness practitioners in certain tasks, with the aim of evaluating any potential link between the impact of this practice on cognitive functions. Lastly, future studies could apply temporal filtering to the data and evaluate the effects of long-term mindfulness practice on specific frequency components of the fMRI signal.