Inducing a meditative state by artificial perturbations: A mechanistic understanding of brain dynamics underlying meditation

Meditation encompasses a wide range of mental training techniques that allow one to explore the inner self and its relation with the outer world (Hanh, 1991; Kabat-Zinn, 2005). Meditation has roots in different cultures and religions, with its origin dating back to ancient Hindu and Buddhist spiritual traditions in India (Dumoulin, 1988; Sharma, 2015; Wynne, 2009). The different practices of meditation (e.g., attentional, open monitoring, loving kindness) have a variety of positive effects, including present-moment awareness and observation of experience with openness, acceptance, and nonattachment (Dahl, Lutz, & Davidson, 2015; Kabat-Zinn, 2003; Lutz, Slagter, Dunne, & Davidson, 2009; Ricard, Lutz, & Davidson, 2014). They can also promote stress reduction, well-being, social connectedness, self-awareness, and attention and emotional regulation (Basso, McHale, Ende, Oberlin, & Suzuki, 2019; Champion, Economides, & Chandler, 2018; Guendelman, Medeiros, & Rampes, 2017; Ludwig & Kabat-Zinn, 2008; Marchand, 2012; Pascoe, Thompson, Jenkins, & Ski, 2017; Tang, Hölzel, & Posner, 2015; Valentine & Sweet, 1999).

Over the past years, meditation is increasingly being explored as an intervention to enhance well-being and alleviate suffering in the healthy general population, as well as to prevent and recover from disability and disease (Dahl et al., 2015; Ludwig & Kabat-Zinn, 2008). The field of contemplative neuroscience, in particular, uses neuroscience tools—for example, functional magnetic resonance imaging (fMRI)—to study the effects of meditation (Davidson & Begley, 2013; Varela, 1993). Neuroanatomical studies have investigated structural gray matter thickness in meditators relative to controls, finding structural evidence such as the thickening of the prefrontal cortex with experience in long-term meditation practice (Lazar et al., 2005). Moreover, functional brain imaging studies have largely shown that key brain areas relevant to meditation belong to the default mode network (DMN), the salience network, and the control network (Fox et al., 2016; Ganesan et al., 2022; Hasenkamp, Wilson-Mendenhall, Duncan, & Barsalou, 2012; Sezer, Pizzagalli, & Sacchet, 2022). Interestingly, these resting-state networks belong to the triple-network model hypothesis which establishes their interaction as core to cognitive function (Menon, 2011). Moreover, this DMN-salience-control network interplay has been described in models of focused attention meditation developed by Hasenkamp et al. (2012) and Ganesan et al. (2022). Different static and seed-based approaches have analyzed blood oxygen level-dependent (BOLD) signal activity and network connectivity changes of meditation in health (Garrison, Scheinost, Constable, & Brewer, 2014; Ives-Deliperi, Howells, Stein, Meintjes, & Horn, 2013; Kilpatrick et al., 2011; Kyeong, Kim, Kim, Kim, & Kim, 2017; Martínez et al., 2021; Pagnoni, 2012; Vishnubhotla et al., 2021) and brain disorders (Datko et al., 2022; Goldin, Ramel, & Gross, 2009; Li et al., 2022; Lifshitz et al., 2019). Furthermore, theoretical methods have considered the underlying brain dynamics (Bremer et al., 2022; Lim, Teng, Patanaik, Tandi, & Massar, 2018; Martínez et al., 2021; Marusak et al., 2018; Mooneyham et al., 2017), including whole-brain approaches such as turbulent dynamics (Escrichs et al., 2022), dynamical complexity (Escrichs et al., 2019) and structural and effective connectivity (De Filippi et al., 2022). Whole-brain dynamic analysis and computational modeling can shed light on the mechanisms underlying meditation practices (Tang et al., 2015). Nonetheless, they remain largely unexplored.

The awakening framework developed by Deco et al. (2019) has been robust and helpful in elucidating the underlying mechanisms and brain transitions in different brain states towards control regimes (Deco et al., 2019; Escrichs et al., 2023; Mana et al., 2023; Vohryzek et al., 2024). It consists of a model-free approach, Leading Eigenvector Dynamics Analysis (LEiDA) (Cabral et al., 2017), and a model-based approach, consisting of Hopf whole-brain computational models and off-line in silico perturbations. Brain dynamics are studied with reference to the concept of metastability, which corresponds to the ability of a system to maintain equilibrium regardless of slight perturbations (Freyer, Roberts, Ritter, & Breakspear, 2012; Kelso, 2012). The LEiDA framework, in particular, typifies a brain state with its probabilistic metastable substate (PMS) space (Cabral et al., 2017; Deco et al., 2019; Escrichs et al., 2021; Figueroa et al., 2019; Kringelbach & Deco, 2020; Lord et al., 2019). Each PMS corresponds to a discrete repertoire of metastable substates (i.e., dynamic patterns) at critical points between chaos and order (Cabral et al., 2017; Deco, Kringelbach, Jirsa, & Ritter, 2017), in which substate duration and arrangement is a dynamic signature of a particular brain state (e.g., sleep, anesthesia) (Deco & Kringelbach, 2016; Tognoli & Kelso, 2014). By building a whole-brain model composed of a network of coupled local nodes (Botvinik-Nezer, Holzmeister, Camerer, et al., 2020; Deco et al., 2019), the empirical PMS can be simulated which then allows to force a transition to a desired control state with in silico external stimulation. Hence, whole-brain reactivity and its mechanistic underpinnings can be studied by stimulating one brain area at a time (Breakspear, 2017; Deco & Kringelbach, 2014; Deco, Tononi, Boly, & Kringelbach, 2015; Kringelbach & Deco, 2020).

Using the framework outlined above, we investigated the underlying mechanisms of brain activity during meditation. As a first step, we applied LEiDA to define the PMS of meditation and rest conditions. As a second step, we developed a Hopf whole-brain model based on the empirical PMS of each condition at the bifurcation point, where oscillatory and noisy regimes are undifferentiated, and the system is critical. Lastly, we applied off-line in silico external, unilateral, and localized probing to force transitions between conditions. Thus, we could evaluate local brain areas in terms of their sensitivity to stimulation and their mechanistic global effects after perturbation. Overall, our work contributes to the state of the art of contemplative neuroscience, as we try to understand brain mechanisms during meditation. Furthermore, it motivates current and potential future therapies based on meditation (Grossman, Niemann, Schmidt, & Walach, 2004; Ludwig & Kabat-Zinn, 2008). This way, we could benefit actual and future practices for healthy population and brain disorders.

This study selected 20 experienced meditators and 20 healthy controls from a dataset previously described in Escrichs et al. (2019). The meditator group (i.e., expert meditators) consisted of 13 males and 7 females (mean age = 39.8 years (SD = 10.29); education = 13.6 years; and meditation experience = 9,526.9 hours (SD = 8,619.8). Meditators were recruited from Vipassana communities in Barcelona, Catalonia, Spain, and had more than 1,000 hours of meditation experience with a daily practice of over one hour. The healthy control group (i.e., controls) consisted of 13 males and 7 females with no prior meditation experience (mean age = 39.75 years (SD = 10.13); education = 13.8 years). Both groups were well-matched for age, gender, and educational level and reported no history of neurological disorder. All participants provided written informed consent, and the study was approved by the Ethics Committee of the Bellvitge University Hospital in accordance with the Helsinki Declaration.

Magnetic resonance imaging (MRI) scans were performed on a 3 T (Siemens TRIO) scanner using a 32-channel receiver coil. High-resolution T1-weighted images were acquired with 208 contiguous sagittal slices, with the following parameters: repetition time (TR) 1,970 ms, echo time (TE) 2.34 ms, inversion time (IT) 1,050 ms, flip angle 9°, field of view (FOV) 256 mm, and isotropic voxel size 1 mm. Resting-state and meditation functional MRI (fMRI) scans were obtained using a single-shot gradient-echo echo-planar imaging (EPI) sequence, comprising a total of 450 volumes, with TR 2,000 ms, TE 29 ms, FOV 240 mm, in-plane resolution 3 mm, 32 transversal slices, thickness 4 mm, and a flip angle of 80°. Diffusion-weighted imaging (DWI) data were collected using a dual spin-echo diffusion tensor imaging (DTI) sequence, comprising 60 contiguous axial slices, with TE 92 ms, FOV 236 mm, isotropic voxel size 2 mm, no gap, and matrix sizes 118 × 118. Diffusion was measured with 64 optimal noncollinear diffusion directions using a single b value of 1,500 s/mm² interleaved with nine nondiffusion b0 images, with a frequency-selective fat saturation pulse applied to avoid artifacts.

In this study, we evaluated expert meditators scanned during meditation and rest, and controls scanned during rest. During rest, participants were instructed not to think about anything in particular. During meditation, participants had to focus on their natural breathing and return whenever their mind wandered (i.e., Anapanasati in Pali language). Participants were asked to remain still and move as little as possible inside the scanner. In the final contrast analysis, we excluded the first 90 volumes of the meditation fMRI condition to ensure that the initial acquisition stage corresponding to the acclimatization to the scanner, which may have interacted with the ability to enter a stable meditative state, did not influence the physiological measures collected. No consensus is established in meditation research with respect to the amount of time needed for entering a stable meditative state, and the decision was based on past studies considering enough remainer recording needed for a high-quality analysis (Braboszcz, Cahn, Levy, Fernandez, & Delorme, 2017; Davanger, Ellingsen, Holen, & Hugdahl, 2010; Manna et al., 2010). For normalization purposes, we applied the same criteria for each resting-state fMRI condition to have equal time points.

Preprocessing was done using MELODIC version 3.14 (Beckmann & Smith, 2004), which is part of FMRIB’s Software Library (FSL, https://fsl.fmrib.ox.ac.uk/fsl). First, the initial 10 volumes were discarded for magnetization stabilization, followed by motion correction using MCFLIRT (Jenkinson, Bannister, Brady, & Smith, 2002), nonbrain removal using BET (Brain Extraction Tool) (Smith, 2002), spatial smoothing with 5-mm FWHM Gaussian Kernel, rigid-body registration, high-pass filter cutoff = 100.0 s, and single-session ICA with automatic dimensionality estimation. Noise removal was performed using FIX (FMRIB’s ICA-based X-noiseifier) (Griffanti et al., 2014), which independently removes noise components for each participant. FSLeyes in Melodic mode was used to manually classify the single-subject Independent Components (ICs) into “good” for signal, “bad” for noise, and “unknown” for ambiguous components, based on the spatial map, the time series, and the temporal power spectrum (Griffanti et al., 2017; Salimi-Khorshidi et al., 2014). Then, FIX was applied using the default parameters to obtain a cleaned version of the functional data. All participants showed low motion, given the mean framewise displacement (FD) was <2 mm (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012). The mean and standard deviation for controls during rest, meditators during rest and meditators during meditation was 0.0149 ± 0.0052, 0.0113 ± 0.0023, and 0.0151 ± 0.0046, respectively.

We then used FSL tools to extract the time series in native EPI space between 100 cortical (Schaefer et al., 2018) and 16 subcortical brain nodes (Tian, Margulies, Breakspear, & Zalesky, 2020). Specifically, we co-registered the cleaned functional data to the T1-weighted structural image using FLIRT (Jenkinson & Smith, 2001). The T1-weighted image was then co-registered to the standard MNI space using FLIRT (12 DOF) and FNIRT (Andersson, Jenkinson, Smith, et al., 2007; Jenkinson & Smith, 2001). The resulting transformations were concatenated, inversed, and applied to warp the resting-state atlas from MNI space to the cleaned functional data in native space, preserving the labels with a nearest-neighbor interpolation method. Finally, we used custom-made MATLAB scripts to extract the averaged time series using fslmaths and fslmeants.

For each participant, a whole-brain structural connectivity matrix (SC) was computed using the atlas described above following the steps of previous studies (Cao et al., 2013; Gong et al., 2009; Muthuraman et al., 2016). The analysis was conducted using FMRIB’s Diffusion Toolbox (FDT) in FMRIB’s Software Library (FSL). First, the DICOM images were converted into Neuroimaging Informatics Technology Initiative (NIfTI) format using dcm2nii. The b0 image in native space was then co-registered to the T1-weighted image using FLIRT (Jenkinson & Smith, 2001). The co-registered T1 image was then co-registered to the standard space using FLIRT and FNIRT (Andersson et al., 2007; Jenkinson & Smith, 2001). The resulting transformation was inverted and applied to warp the atlas from the MNI space to the native MRI diffusion space using nearest neighbor interpolation. Next, the diffusion-weighted images were processed using the FDT pipeline in FSL. Nonbrain tissues were removed using the Brain Extraction Tool (BET) (Smith, 2002). The eddy current distortions and head motion were corrected using the eddy correct tool (Andersson & Sotiropoulos, 2016), and the gradient matrix was reoriented to correct for motion (Leemans & Jones, 2009). To model crossing fibers, BEDPOSTX was used, and the probability of multifibre orientations was computed to improve the sensitivity of nondominant fibre populations (Behrens, Berg, Jbabdi, Rushworth, & Woolrich, 2007; Behrens et al., 2003). Probabilistic Tractography was then performed for each participant in their native MRI diffusion space using the default settings of PROBTRACKX (Behrens et al., 2003, 2007). Finally, the connectivity probability SC_np between nodes n and p were calculated as the proportion of sampled fibres in all voxels in node n that reach any voxel in node p. Since diffusion tensor imaging (DTI) does not capture fiber directionality, the SC_np matrix was symmetrized by computing its transpose matrix SC_pn and averaging both matrices.

When a pair of nodes are temporarily aligned, the phase coherence is close to one since the difference between their Hilbert transformed signal angle is 0° [cos(0°) = 1]. Contrarily, orthogonal BOLD signals have a phase coherence close to zero [cos(90°) = 0]. The resulting dFC(t) of each participant represented the interregional BOLD synchrony at each timepoint. The size of the matrix was N × N × T, where N is the number of brain areas (116) and T the total time points (350). A total of 14,000 (20 participants * 350 timepoints * 2 conditions) matrices (dimension N × N) were calculated. A dominant connectivity pattern was obtained by reducing the dimensionality of the undirected and symmetric matrices into their leading eigenvectors V₁(t) (dimension N × 1) (Deco et al., 2019). The eigenvectors capture the dominant connectivity pattern at each time point t, representing the contribution of each brain area to the whole structure (Cabral et al., 2017). It is possible to reproduce the dominant connectivity pattern of dFC(t) by calculating the outer product of V₁(t) and its transpose (V1.V1.^T) (Lohmann et al., 2010). The leading eigenvectors of the dFC(t) for each TR and all participants from all states were then clustered with K-means clustering, varying k from 2 to 8. The resulting cloud centroids V_c(t), which represent the dominant connectivity pattern in each cluster, were rendered onto a brain map using the Surf Ice software (https://www.nitrc.org/projects/surfice/).

Lastly, we computed the Probabilistic Metastable Substate Space (PMS) for each condition (i.e., rest and meditation). The PMS includes all the metastable substates extracted in LEiDA. Thus, each condition is characterized by a unique PMS corresponding to the probability of being in each metastable substate. Specifically, the ratio of number of TRs in a cluster divided by the total number of TRs in all clusters was computed separately for resting-state and meditation.

A whole-brain computational model was built for meditation and resting-state. These were generated using the normal form of a supercritical Hopf bifurcation. This Landau-Stuart oscillator, referred to as Hopf, has been used to study transitions from noisy to oscillatory regimes (Deco et al., 2015) (Figure 1B). Emergent brain dynamics were simulated by linking the structural and functional connectivity using effective connectivity (Deco et al., 2015). For each of the 100 cortical and 16 subcortical brain areas (Schaefer et al., 2018; Tian et al., 2020), BOLD activity was emulated and the working point was fitted and optimized with specific parameters of the model (Deco et al., 2019).

The Hilbert transform BOLD signal phase of nodes j and i at time t is expressed by φ(t) and the brackets indicate the average across time. Each time, the empirical and simulated values were compared until the difference was less than 0.001 (Deco et al., 2019).

Transitions between each source state and a target state (i.e., from resting-state toward meditation and vice versa) were evaluated as schematized in Figure 1B. The Hopf model allows shifting the bifurcation parameter of a particular brain area while keeping the rest of the regions at their bifurcation value corresponding to the optimal fit of the model (i.e., a = −0.02), thus representing a local perturbation. In other words, the whole-brain model for a source state was perturbed unilaterally by moving the local bifurcation parameter a of one brain area individually and separately while maintaining the rest fixed to their initial value. As a consequence of a local stimulation, the global brain dynamics of the perturbed model changed and were analyzed by calculating the perturbed PMS using the empirical cluster centroids found during LEiDA analysis. A measure of the fit to the target state was evaluated by computing the KL distance between the perturbed PMS and the target empirical PMS. The perturbation was repeated independently for each of the 116 brain areas. The areas with the lowest KL distance were the optimal to promote transitions. In order to minimize random effects from the Gaussian noise of the model, each simulation was repeated three times (Deco et al., 2019). There were two different protocols studied: synchronization and noise. For the synchronization protocol, the bifurcation parameter a of each brain area was shifted from −0.02 to 0.18 in steps of 0.01, given that for positive values of a there is a stable limit cycle oscillation. As for the noise protocol, the bifurcation parameter a of each brain area was shifted negatively from −0.02 to −0.22 in steps of 0.01 for the local dynamics to have a low activity noisy state. Stimulation intensities are represented by the absolute value of each step (Deco, Kringelbach, et al., 2017).

The statistical analysis was performed using MATLAB R2022a from MathWorks (Natick, MA, USA). To test the results of the PMS obtained by LEiDA, Wilcoxon tests with 1,000 iterations were used to examine the probability of occurrence of all explored clustering conditions (k from 2 to 8). Permutations were compared using the Wilcoxon test with a significance threshold of 0.05. Multiple comparisons were corrected with the False Discovery Rate (FDR) method (Hochberg & Benjamini, 1990) when testing the differences between resting-state and meditation for all metastable substates in each value of k. The reported p values remain significant after FDR correction.

We first analyzed the group of expert meditators during resting-state and meditation. We selected the minimum number of k that provided the maximum significant differences between the probabilities from each condition for each given substate (i.e., k = 5). This selection survived correction by multiple comparisons and was the maximum value acceptable for the model-based analysis in terms of computational costs. We illustrate the probability of occurrence of the PMS of resting-state and meditation in Figure 2A. In addition, we rendered cluster centroid eigenvectors onto the brain maps where the sign of eigenvector elements are positive or negative and indicated by red and blue colors, respectively (Figure 2B). The signs represent the direction that the BOLD signals from different elements follow in relation to their projection onto the leading eigenvector. Hence, classifying elements into communities in which color strength reflects how strongly each area is associated with its community (Alonso Martinez, Deco, Ter Horst, & Cabral, 2020; Cabral et al., 2017). In the first metastable substate, all eigenvector elements had the same sign (negative). On the other hand, substate 2 showed a functional network dominated mainly by somatomotor and salience networks. Additionally, substate 3 displayed a community comprised mostly of the DMN and a few areas from the limbic network, whereas substate 4 showed coordination primarily between the control network and the DMN. Lastly, substate 5 had a network principally led by the visual system and some areas of the somatomotor network. Significant differences between conditions were found in the probability of occurrence of substates 2 and 3. For substate 2, the probability of occurrence was lower during meditation compared to resting-state [0.1364 ± 0.0129 vs. 0.1807 ± 0.0150, p = 0.0157]. The same applied for substate 3, given it was visited less often during meditation compared to rest [0.1196 ± 0.0136 vs. 0.1696 ± 0.0175, p = 0.0161]. Different clustering configurations can be consulted in Supporting Information Figure S1 and Figure S2. As a further analysis, we investigated the relationship between the probability of occurrence of each substate and the hours of experience of expert meditators. We calculated this for both rest and meditation, in separate analysis, in order to look whether the experience of experts meditators affected the probability of any substate during any condition. We found no significant correlations (Supporting Information Figure S3). Lastly, we looked for the correlation between the condition delta FD (meditation - rest) and the delta substate occupancy for each substate finding no significant results (Supporting Information Figure S4). This confirmed that there was no relationship between the head motion during the scan and the probability of occurrence of the substate led by the somatomotor network. Overall, LEiDA allowed us to elucidate distinct patterns of brain dynamics providing a more complete disentangling of brain configuration compared to classical static analysis (Supporting Information Figure S5).

We then analyzed differences between controls during rest and expert meditators during meditation. We found that the minimum number of k that provided a maximal proportion of statistically significant differences between conditions was k = 2 (Figure 2C). The centroid eigenvectors for this clustering configuration are rendered into brain maps in Figure 2D. The metastable substate 1 had all negative eigenvector elements, whereas substate 2 had a coordination led by brain areas from the DMN, and some from the control network and limbic network. All substates revealed significant differences between conditions. In the case of substate 1, the probability of occurrence during meditation was higher compared to rest [0.7099 ± 0.0221 vs. 0.6291 ± 0.0258, p = 0.0246]. Conversely, substate 2 was visited less often during meditation compared to rest [0.2901 ± 0.0221 vs. 0.3709 ± 0.0258, p = 0.0327]. Different clustering configurations can be consulted in Supporting Information Figure S6. Lastly, we sought to determine if there is a relationship between the probability of occurrence of each substate during meditation and the hours of practice of expert meditators. No significant correlation was found (Supporting Information Figure S3). Furthermore, the correlation between the condition delta FD (meditation - rest) and the delta substate occupancy for each substate did not have any significant result (Supporting Information Figure S4). Head motion was not associated with substate occupancy.

We fitted and optimized a whole-brain model for each condition in order to find the parameters that maximally approximated the corresponding empirical PMS. This was identified with the minimal KL distance between the simulated and the empirical PMS. The best fit for the models of meditation and resting-state were found at a global coupling of G = 0.05 with a KL distance of 0.0122, and G = 0.04 with a KL distance of 0.0064, respectively. The similar probabilities of occurrence for all five substates are shown in Figure 3A and Figure 3B. The modeled PMS for k = 4 can be consulted in Supporting Information Figure S7. The evolution of G as a function of the KL distance for each model can be consulted in Supporting Information Figure S8.

The best fit for the models of meditation and resting-state was found at a global coupling of G = 0.04 with a KL distance of 0.0002, and G = 0.04 with a KL distance of 0.0003, respectively. The simulated probabilities of the two substates are shown in Figure 3C and Figure 3D. The evolution of G as a function of the KL distance for each model can be consulted in Supporting Information Figure S8.

In the last step of the analysis, we systematically perturbed the model of each source state and compared the resulting perturbed PMS with the empirical PMS of the target state. For the synchronous protocol, we shifted the bifurcation parameter a independently in each brain area, with positive values. Contrarily, for the noise protocol, we used negative values. We identified the optimal perturbation as the one resulting in the smallest KL distance between the perturbed modeled source PMS and the target empirical PMS.

Results for perturbing the whole-brain model of resting-state in the synchronization and noise protocols at different stimulation intensities are illustrated in Figure 4A. The color-scale indicates the KL distances between the empirical and the perturbed PMS after stimulating each brain area separately. A transition was successful in the synchronization protocol, where KL distances decreased with increasing positive bifurcation point values. The minimum KL distance for the optimal perturbation of each brain area is rendered into the brain maps of Figure 4B. The top 10 areas to perturb, at their particular optimal stimulation intensity in the synchronization protocol, belonged mainly to the somatomotor network, posterior dorsal attentional network, and medial salience network in decreasing order of prevalence (Supporting Information Table S1). We achieved the best fit when perturbing the third segment from the left somatomotor cortex with a = 0.11. Consequently, the perturbed and target PMS were very similar, as reflected in Figure 4C.

On the other hand, transitions from meditation towards resting-state were successful in the noise protocol and not in the synchronization protocol. Here, KL distances decreased with decreasing negative values of the bifurcation parameter a (Figure 4D). The optimal stimulation for each brain area at a particular optimal stimulation intensity for the noise protocol can be seen in Figure 4E, in which the top 10 areas to perturb were found mainly in the DMN, followed by the control and limbic networks as well as the left anterior thalamus (Supporting Information Table S1). The optimal transition was found for the second segment of the left parietal region from the DMN for a = −0.14 (Figure 4F). Perturbational analysis for k = 4 can be consulted in Supporting Information Figure S9 and Supporting Information Table S2.

Transitions from rest to meditation were successful with the synchronization protocol Figure 5A. The optimal perturbation for each brain area in the synchronization protocol is shown in Figure 5B. Here, the top 10 areas to perturb were mainly from the DMN, and some brain areas from the somatomotor, limbic and control networks, and the subcortical areas left amygdala and left hippocampus (Supporting Information Table S3). The optimal transition was found when stimulating the third segment of the left prefrontal cortex from the DMN for a = 0.09 (Figure 5C).

Furthermore, results for stimulating the meditation model are illustrated in Figure 5D. Transitions were achieved in the noisy protocol and not in the synchronization protocol, with decreasing KL distance with more negative stimulation intensities. The best transition for each brain area is shown in Figure 5E, and the top 10 areas to perturb were mainly from the DMN, and the rest from the control network (Supporting Information Table S3). The optimal transition was found when perturbing the second segment of the right precuneus posterior cingulate cortex from the DMN for a = −0.20 (Figure 5F).

We successfully studied brain dynamics using model-free and model-based approaches (Deco, Kringelbach, et al., 2017) to find mechanistic explanations for meditation. LEiDA allowed us to characterize brain signatures of meditation and resting-state by defining their PMS, with an associated probability of occurrence of being in distinct metastable substates (i.e., dynamic patterns) (Cabral et al., 2017). Then, we fitted a Hopf whole-brain computational model to the empirical PMS of each state. Consequently, this allowed us to force transitions from a modeled PMS to an empirical target PMS using in silico perturbations. By varying stimulation intensities and protocols (i.e., noise and synchronization), we were able to determine the mechanistic global effects of exhaustive local perturbations. This way, we revealed how local changes can alter whole-brain dynamics and evaluated the optimal transitions toward a given state. Using a synchronization protocol, we could transition from the resting-state to the meditative state. Furthermore, using a noisy protocol, we could transition from meditation to rest.

In the model-free approach, we used LEiDA to reveal signatures of functional networks underlying meditation and resting-state (Figure 2). In the analysis of expert meditators, we found a significantly lower probability of occurrence in substates 2 and 3 during meditation compared to rest. Firstly, we found that substate 2 was led by the somatomotor and salience networks. The somatomotor network is responsible for processing bodily sensations and motor functions (Biswal, Zerrin Yetkin, Haughton, & Hyde, 1995). According to the review of Ganesan et al. (2022), most published evidence on areas from the somatomotor network showed deactivation during meditation. We suggest that the decrease in this network could align with the disengagement from external sensory interferences in the practice. Furthermore, a reduction of somatomotor cortex activity has been associated with decreased pain-related perception in a meditative state (Nakata, Sakamoto, & Kakigi, 2014; Zeidan et al., 2011). Additionally, the salience network monitors the switching between internal and external processing (Elton & Gao, 2014; Goulden et al., 2014). Reduced functional connectivity within the salience network has been associated with experiential acceptance and present-moment awareness (Sezer et al., 2022). Guidotti et al. (2023) used a trained classifier revealing that this network was able to discriminate meditation styles in an expert group of meditators. On the other hand, substate 3 displayed a community led majorly by the DMN and some areas of the limbic network. The DMN is involved in internal thought, self-related processing, and mind wandering (Alves et al., 2019; Mason et al., 2007). Our result aligns with previous studies reporting a reduction in DMN activity during meditation (Ganesan, Moffat, Van Dam, Lorenzetti, & Zalesky, 2023; Garrison, Zeffiro, Scheinost, Constable, & Brewer, 2015). This reduction has been associated with the phenomenological well-being of meditation (Brandmeyer & Delorme, 2021; Killingsworth & Gilbert, 2010). Furthermore, relevant areas from the limbic network were the orbitofrontal cortex and the temporal pole, involved in emotional processing. The orbitofrontal cortex is engaged in reward-based tasks (Kringelbach & Rolls, 2004), whereas the temporal pole has a role in social processing (Herlin, Navarro, & Dupont, 2021). In turn, due to increased attention to the breath, we suggest the reduction in probability of this substate may decrease emotional and/or reward-related processes.

In the analysis of controls and expert meditators, we found a significantly higher probability of occurrence of substate 1 and a significantly lower probability of occurrence of substate 2 during meditation compared to rest. The first substate, characterized by all negative eigenvector elements, has been associated to synchronization, stability, or noise artifacts (Farinha, Amado, Morgado, & Cabral, 2022; Olsen et al., 2022). On the other hand, substate 2 resembled substate 3 from the previous analysis, and could be also explained by lower mind wandering (Mason et al., 2007) and emotional interference (Ganesan et al., 2022) during meditation.

For the model-based approach, we first built a Hopf whole-brain computational model for each condition. This allowed us to systematically evaluate in silico stimulations to force transitions between resting-state and meditation (Deco, Kringelbach, et al., 2017). In particular, we fitted each model to the empirical data and identified the optimal global coupling value, corresponding to the lowest KL distance between the empirical PMS and the modeled PMS (Figure 3). We then explored local perturbations area by area in order to study mechanistic global effects. We evaluated a synchronization protocol (i.e., toward a system with self-sustained oscillations) and a noise protocol (i.e., toward a system that decays to a stable point and low noisy activity). In the analysis of expert meditators, transitions from resting-state to meditation were successful with a synchronous protocol (Figure 4A–C). Our results show that the most sensitive areas to perturbation were found in the somatomotor network, posterior regions of the dorsal attention network, and medial areas from the salience network. Interestingly, these brain areas overlapped with substate 2 from the LEiDA analysis. In the somatomotor network, a sensitive area was the supplementary motor area, related to inner speech processing (Hertrich, Dietrich, & Ackermann, 2016). Its role here could be supported by the shift from internal distraction during meditation. Moreover, the dorsal attention network is associated with self-directed attention (Devaney et al., 2021). Its increased connectivity has been related to mindfulness and focused attention (Froeliger et al., 2012; Zhang et al., 2021). As a final point, the sensitivity of the salience network is highly relevant given it is hypothesized to play a modulator role in neurocognitive models of focused attention meditation (Ganesan et al., 2022; Hasenkamp et al., 2012). On the other hand, our results show that transitions from meditation to the resting-state were achieved in the noise protocol (Figure 4D–F). Here, the most sensitive areas belonged to the DMN, overlapping substate 3 in the LEiDA analysis. The role of the DMN in forcing a transition could be related to mind wandering and increased DMN activity at rest (Mason et al., 2007).

For the model-based results of controls and expert meditators, transitions from rest to meditation were also effective in the synchronization protocol (Figure 5A–C). Furthermore, highest sensitivity was found in the DMN and subcortical areas (i.e., hippocampus and amygdala). The sensitivity of the DMN for forcing a transition could be associated with its active role during rest, which decreases in meditative states (Mason et al., 2007). On the other hand, the hippocampus and amygdala have been associated with the emotional brain. The hippocampus is involved in emotional behavior, learning, memory, and spatial navigation (Anand & Dhikav, 2012). In addition, the amygdala plays a key role in the assignment of sensory information with emotional value (Šimić et al., 2021). Thus, their sensitivity in the transitions could be explained by the need to reduce interference from emotional information during meditation. On the other hand, our results show that transitions from meditation to the resting-state were achieved with a noisy protocol (Figure 5D–F). The most sensitive areas belonged to the DMN. The importance of the DMN has been revealed throughout the whole analysis, and the perturbational results show its pivotal influence on shifting transitions between rest and meditation. Lastly, the remaining most sensitive areas belonged to the central executive network. This network has a role in executive functions (i.e., attention, cognitive control, and working memory) (Niendam et al., 2012), and its influence on the transition could be related to the shift in cognitive demands from meditation to rest.

In all cases, a transition from rest to meditation was successful in a synchronization protocol, and the opposite transition occurred in a noise protocol. A fundamental brain function is short- and large-scale neural communication achieved with synchronization of oscillations between different brain areas (Schnitzler & Gross, 2005). Synchronization is essential to cognition, and abnormal values are associated with brain disorders (Uhlhaas & Singer, 2006). The relationship between synchronization and a meditative state is not straightforward. EEG literature of expert meditators during meditation showed increased gamma synchronization (Berkovich-Ohana, 2017; Lutz, Greischar, Rawlings, Ricard, & Davidson, 2004) and higher betweenness centrality in the alpha band (van Lutterveld et al., 2017). These could support our results given they were related to a restructuring in information processing and an increase of integrative mechanisms of information exchange for facilitating a meditative state (Engel, Fries, & Singer, 2001; Fell, Axmacher, & Haupt, 2010). In contrast, the transition from meditation to rest in the noisy protocol reflects that neural noise enables the emergence of the dynamic repertoire of the brain at rest (Ghosh, Rho, McIntosh, Kötter, & Jirsa, 2008), promoting less synchronization across the whole brain (Deco et al., 2019). Overall, our model-based results provide novel insights into how local perturbations influence global brain dynamics. We could explore different stimulation protocols and the most sensitive brain areas in terms of their response to perturbations as a model-based biomarker of the underlying mechanisms of meditation.

This study has some limitations. Firstly, we analyzed a relatively small sample size dataset. Therefore, validating our results using a larger dataset would be beneficial to increase the generalizability of our findings (Fox et al., 2016). We also want to acknowledge the fact that no established consensus on the time for the preparatory period in meditation research is yet established. The choice for the time removed due to the build-up period was drawn from past studies (Braboszcz et al., 2017; Davanger et al., 2010; Manna et al., 2010). In addition, our analysis focuses on a specific meditation technique that prevents the translation of our findings to all styles of meditation. Thus, we cannot oversimplify our results to one signature of brain dynamics in meditation given the heterogeneity of the practice (Fox & Cahn, 2019; Young et al., 2021). Moreover, complementary studies would be necessary to evaluate the longitudinal effects of meditation practice on brain structure and function.

In this work, we explored the underlying brain dynamics in meditation. Using the LEiDA model-free approach, we could quantitatively define the brain dynamics during meditation and resting-state. In addition, the model-based approach allowed us to measure the whole-brain reactivity to localized perturbation at resting-state and study transitions toward meditation and vice versa. As a last remark, we want to emphasize the importance of not oversimplifying the Eastern tradition of meditating to a single motivation around attention and emotion regulation (Garcia-Campayo, López Del Hoyo, & Navarro-Gil, 2021). And, if used as a therapeutic intervention, it is essential to recognize that meditation is not a cure-all solution (Wright, 2013). Overall, given the paucity of research in meditation with neuroscience tools, our results shed light on the field of contemplative neuroscience for promoting meditation practices in health and a potential therapy in disease.

We are sincerely grateful to Mind and Life Europe for supporting this project with the Francisco J. Varela Award. We thank participants for their involvement in the study.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00366.

Paulina Clara Dagnino: Conceptualization; Formal analysis; Investigation; Methodology; Project administration; Software; Validation; Visualization; Writing—original draft; Writing—review & editing. Javier A. Galadí: Writing—review & editing. Estela Càmara: Data curation; Resources; Writing—review & editing. Gustavo Deco: Conceptualization; Funding acquisition; Investigation; Methodology; Project administration; Resources; Software; Supervision; Validation; Writing—review & editing. Anira Escrichs: Conceptualization; Data curation; Formal analysis; Funding acquisition; Investigation; Methodology; Project administration; Resources; Software; Supervision; Validation; Writing—original draft; Writing—review & editing.

P.D. was supported by the AGAUR FI-SDUR Grant (no. 2022 FISDU 00229) and by the AGAUR research support grant (ref. 2021 SGR 00917) funded by the Department of Research and Universities of the Generalitat of Catalunya. J.G. was supported by the Spanish Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación (FJC2021-047765-I). A.E. and G.D. were supported by the Grant PID2022-136216NB-I00 funded by MICIU/AEI/10.13039/501100011033 and by “ERDF A way of making Europe,” ERDF, EU. G.D. was also supported by the project NEurological MEchanismS of Injury, and Sleep-like cellular dynamics (NEMESIS) (ref. 101071900) funded by the EU ERC Synergy Horizon Europe. A.E. was also supported by the project eBRAIN-Health - Actionable Multilevel Health Data (id 101058516), funded by the EU Horizon Europe.

Functional magnetic resonance imaging (fMRI):

Noninvasive neuroimaging technique for mapping brain activity in vivo based on relative changes in blood oxygen levels in the capillary beds.

Structural connectivity:

Anatomical links between different brain areas, often measured with diffusion tensor imaging.

Effective connectivity:

Directed influence of one brain area on another, not necessarily showing direct connections but measuring the interactions between different brain areas.

Brain state:

Evolving pattern of self-organized neural activity specific to a condition.

Metastable substate:

Dynamic large-scale functional network patterns that change over time, stochastic subdivisions of brain states.

Diffusion tensor imaging (DTI):

Noninvasive magnetic resonance imaging technique that measures the brain’s white matter tracts in vivo based on the diffusion properties of water through nerve fibers.

Phase coherence matrix:

Cosine difference of phases between each pair of brain areas.

Hopf bifurcation:

In nonlinear dynamics, critical point where a system shifts from a fixed stable point toward a self-sustained oscillatory regime.

Functional connectivity:

Statistical dependencies between different brain areas, often measured with the coherence and correlation between the time series from fMRI.