Introduction

Mind-body interventions–structured practices that harness the interaction between psychological processes and physiological systems–can significantly improve human physical and mental health, yet the neural and molecular effects of increasingly popular mind-body retreat interventions remain poorly understood. In a recent RCT, reconceptualizing pain as the product of plastic brain activity rather than of peripheral tissue injury reduced pain three times more than placebo and six times more than usual treatment1. Placebo effects, another mind-body intervention based on healing-centered rituals, symbols, and behaviors, affect every major organ system2,3 and sometimes surpass routine surgical outcomes4,5,6. Interestingly, open-label placebos (administered without concealment such that the subject is aware of the placebo) conserve their effectiveness against many conditions7,8,9, demonstrating that placebo responses do not require deception, expectations, or conditioning. Meditation, yet another intervention based on self-regulated attentional practices, can produce subjective mystical-type experiences10, mental health improvements11,12, and can alter neural13,14,15, immune16, autonomic17,18 gene expression19,20,21, proteomic22,23,24, and metabolomic25,26 activity.

Each of these interventions operates through partially distinct cognitive mechanisms, raising the possibility that they may complement one another to impact the brain and body synergistically. Reconceptualization operates through conscious, discursive, and volitional alteration of core beliefs, while meditation involves a willful yet non-discursive alteration of consciousness, and open-label placebo involves conscious awareness but operates unconsciously. While research on each technique exists, their combined neural and molecular effect has never been studied. To do so, we conducted an exploratory observational study with functional magnetic resonance imaging (fMRI) and blood plasma-based high-throughput proteomics, metabolomics, exosome-specific miRNA transcriptomics; and neurite growth and real-time metabolism cellular assays (Fig. 1A) on 20 healthy adult participants (14 females, age = 46.35 ± 10.06 (SD) years) (Fig. 1B) sampled before and after a 7-day mind-body retreat (Fig. 1C). While logistical limitations prevented us from including age and gender-matched controls to isolate and mechanistically describe the neural and physiological pathways engaged by each separate mind-body technique–an equally important but separate task–, our study provides evidence for the breadth and depth of physiological effects following a holistic, multifaceted experience commonly described by participants as personally transformational.

Fig. 1: Study design, participants, data collection, and recruitment.

A Outcome measures to capture biological changes associated with brain and body. Created with BioRender. Simpson, S., Jinich, A. (2025). BioRender.com/ryzs1cd. B Participant characteristics, including age, subject, gender, meditation experience level, and top/bottom median scores for BDNF, short-chain fatty acid (SCFA) metabolism, and HDAC proteomic pathways. C Intervention and data collection timeline. (fMRI functional magnetic resonance imaging, MEQ mystical experience questionnaire). D CONSORT flow diagram.

The 7-day retreat combined lectures, meditation, and healing rituals. Daily lectures (25 total hours) emphasized the body’s self-healing abilities, the mind’s capacity to shape lived reality, and the healing power of present-centeredness and mystical-type experiences. All meditations (33 total hours) were guided, delivered with atmospheric music, and taught Kundalini techniques, which combine conscious meta-awareness and conscious breathing exercises with slow, ascending, focused interoceptive attention on purported energetic centers along the midline (e.g., brow, throat, heart) which, according to practitioners, can reprocess embodied trauma and catalyze adaptive mental and physical changes27,28. Guidance also emphasized sustaining a heart-centered state devoid of thinking or judgment and focusing awareness on a void beyond one’s normal sense of space and time—a common theme in some contemplative practices29. Guided healing rituals (5 total hours) brought 6–8 “healers” around one “healee” in which the former were instructed to practice loving-kindness compassion meditation while focusing attention on their heart, hands, and on the latter’s body. A healing mechanism was not presented, but the possibility that healing could occur on either party because of the ritual was mentioned, similarly to how open-label placebos are presented in trials30. All study subjects participated as healers. Of the 20 participants, 11 were “advanced” meditators who had practiced the techniques taught for at least six months, while 9 were “novices” who had not. No pharmacological substances, including any psychedelics, were involved in the retreat.

We characterized the resting and meditation states experienced by participants pre- and post-intervention. Consistent with previous meditation studies[31](https://www.nature.com/articles/s42003-025-09088-3#ref-CR31 “Prakash, R. S. et al. Mindfulness meditation and network neuroscience: review, synthesis, and future directions. Biol. Psychiatry Cogn. Neurosci. Neuroimaging https://doi.org/10.1016/j.bpsc.2024.11.005

(2024).“),32, we observed higher meditation-associated whole-brain functional integration and reduced intra-network connectivity in the default mode and salience networks. Comparing pre- to post-intervention whole plasma, we report evidence of enhanced neuroplasticity and glycolytic metabolism, as well as activation of endogenous opioid and neuromodulatory pathways.

Results

Functional brain activity

To characterize the neural signature of the meditative state, participants underwent structural and blood-oxygenation-level-dependent (BOLD) functional MRI scans during rest (5 min) and meditation (15 min). Mystical experience questionnaire (MEQ) scores reflective of the scanner meditation increased significantly following the intervention (n = 20; pre = 2.37 ± 1.26 (SD); post = 3.02 ± 1.44) (Wilcoxon signed rank test W = 35.0, p = 0.03), indicating a deepening of the meditative state between sessions. One participant’s data was excluded from further analysis due to corrupt T1w data. Participants moved more during meditation than rest, a potential confound revealed by the significant effect of task (meditation, rest) on mean framewise displacement on a two-way (task × time) repeated measures ANOVA (n = 19, F(1,18) = 25.1, p = 0.00009, η²p = 0.58) (Supplementary Tables 1 and 2).

Meditation versus rest

We examined functional connectivity in seven resting state networks (RSNs) (Fig. 2A, B), eight regions of interest (ROIs) (Fig. 2C) (coordinates in Supplementary Tables 3 and 4), and two whole-brain networks (Fig. 2D), comparing meditation with rest (n = 19, pFWE = 0.0018 for 28 network pairs). All resting state networks had higher intra- than inter-network connectivity (Supplementary Fig. 1A). Compared to rest, meditation reduced intra-network connectivity in the salience network pre- and post-intervention (pre: t = −4.32, p = 0.0004, Cohen’s d = −0.87; post: t = −6.71, p = 0.000003, Cohen’s d = −1.76), default mode network post-intervention (DMN) (t = −5.04, p = 0.00009, Cohen’s d = −1.78), and dorsal attention network (DAN) post-intervention (t = −3.71, p = 0.002, Cohen’s d = −1.06) (Fig. 2A). Meditation also showed reduced inter-network connectivity post-intervention between the salience and dorsal attention networks (t = −4.09, p = 0.0007, Cohen’s d = −1.18) and between the somatomotor and auditory networks (t = −5.30, p = 0.00005, Cohen’s d = −1.38).

Fig. 2: Functional magnetic resonance imaging (fMRI) (n = 19 participants).

A Meditation vs. rest RSN functional connectivity pre- and post-intervention. Cell values and color scale represent Cohen’s d effect size values from paired t-tests. Asterisks denote pFWE < 0.05 significance (28 multiple comparisons). ECN (executive control network), DMN (default mode network), DAN (dorsal attention network), SN (salience network), SMN (somatomotor network), VN (visual network), and AN (auditory network). B Meditation vs. rest RSN-region functional connectivity pre- and post-intervention. Cell values and color scale represent Cohen’s d effect size values from paired t-tests. Asterisks denote pFWE < 0.05 (630 multiple comparisons). C Significant (p < 0.05) meditation-induced functional connectivity changes between a priori ROIs pre- and post-intervention. Color scale represents Cohen’s d effect sizes from paired t-tests. D Whole-cortex/brain network measures for 48-region Harvard-Oxford cortical atlas and Power (2011) 264-region brain atlas (n = 19). Asterisks denote * (p < 0.05), ** (p < 0.01), *** (p < 0.001), **** (p < 0.0001). E Anatomical Structure. Coronal and axial views with superior parietal lobule cluster (yellow), where advanced practitioners had greater gray matter volume than novices at baseline.

For individual network-component regions (Fig. 2B and Table 1), both pre- and post-intervention, meditation reduced connectivity between the right posterior cerebellum and medial prefrontal cortex (pre: t = −4.44, p = 0.0003, Cohen’s d = −1.15; post: t = −5.81, p = 0.00002, Cohen’s d = −1.30) and bilateral connectivity between insular cortices, lateral parietal cortices, and anterior prefrontal cortices (all p ≤ 0.0003; Table 1), network hubs associated with sensory, affective, and cognitive functions that shape the prediction-driven conscious state33. Interestingly, the only meditation-driven connectivity increase was post-intervention between the left insula and posterior cingulate cortex (t = 4.50, p = 0.0003, Cohen’s d = 1.47), a finding previously reported during absorptive trance states34. ROIs also revealed lower functional connectivity between key network hubs, including medial prefrontal cortex (mPFC), precuneus, bilateral insular cortices, and bilateral angular gyri (Supplementary Table 5), serving as an internal replication and returning results agreeable with the RSN data-driven approach.

To observe meditation-induced changes at the whole-brain/cortex level, we parcellated the brain into 48 cortical regions using the Harvard-Oxford atlas35,36,37 and into 264 brain regions using the Power atlas38 and calculated whole-brain/cortex network measures (Supplementary Table 6). For the cortical parcellation, a 2 (pre/post) by 2 (rest/meditation) repeated measures ANOVA (Supplementary Table 7) (n = 19) revealed a significant effect of meditation on network modularity (F(1,18) = 15, p = 0.001, η²p = 0.45) and global efficiency (F(1,18) = 48, p = 0.000002, η²p = 0.73), and no significant effect on characteristic path length, with similarly significant whole-brain results. Post-hoc Wilcoxon signed rank tests confirmed that meditation decreased modularity (pre: W = 18.0, p = 0.001, Cohen’s d = −1.07; post: W = 45.0, p = 0.04, Cohen’s d = −0.77) and increased global efficiency (pre: W = 3.0, p = 0.00002, Cohen’s d = 2.04; post: W = 10.0, p = 0.0002, Cohen’s d = 1.26) as compared to rest (Fig. 2D). These results were robust to excluding BOLD runs with mean framewise displacement > 0.3 mm, indicating they were not due to higher meditation-associated head motion.

Meditation thus reduced functional connectivity in the default mode and salience networks and induced a whole-brain functional state less segregated into distinct modules, that allows for more efficient information flow.

Pre versus post

Pre- to post-intervention, functional connectivity during meditation decreased between executive control and salience networks (t = −2.61, p = 0.02, Cohen’s d = −0.62) but did not survive correction for multiple comparisons (pFWE = 0.0018) (Supplementary Fig. 1B). For network component regions, connectivity during meditation decreased significantly pre-to-post between the left posterior intra-parietal sulcus and auditory cortex and between the left posterior intra-parietal sulcus and left lateral parietal cortex (Table 2). No significant changes were found between a priori ROIs (Supplementary Table 8).

Neuroanatomical differences

No significant anatomical changes were observed pre-to-post intervention, but advanced practitioners showed greater gray matter volume in the right superior parietal lobule at baseline (n = 19, p = 0.049) (Fig. 2E), a region linked to spatial awareness and body representation, which has been reported to have greater cortical thickness in experienced Zen meditators39. This difference warrants further investigation with larger samples and comparisons with age-matched non-meditator controls.

Whole plasma

Having characterized the neural state induced by meditation, we turned our attention to the broader metabolic and molecular changes in whole plasma.

Enhanced neuroplasticity

Participants’ anecdotal reports consistently emphasize radical psychological breakthroughs, and previous meditation studies have reported increased BDNF (brain-derived neurotrophic factor) levels22 consistent with enhanced neuroplasticity. To investigate whether the intervention affected circulating plasma factors conducive to neuroplasticity, we treated cultured glutamatergic PC12 neuroendocrine cells with NGF (nerve growth factor) and 1% pre- and post-intervention plasma and quantified neurite outgrowth length. Starting on day 4 post-NGF treatment, post-plasma-treated cells exhibited significantly longer neurites than pre-plasma-treated cells (t = −2.52, p = 0.01, Cohen’s d = 0.59) and continued to do so until the end of the experiment (Fig. 3A, B).

Fig. 3: Neuronal growth (n = 20 participants).

A PC12 neurite growth: Percent change in mean neurite length from baseline (day 2 post-NGF treatment). Asterisks denote * (p < 0.05), ** (p < 0.01), *** (p < 0.001), **** (p < 0.0001). B Phase contrast 20**×** microscopy images of PC12 cells on day 10 post-NGF treated with pre- and post-intervention plasma with longest neurite per cell body traced. C BDNF pathway index pre-to-post fold change. Error bars denote SEM. Fold change levels significantly above zero signal upregulated pathway or protein cluster. D BDNF Pathway Heat Map. Fold change per protein for advanced and novice participants.

To investigate proteomic factors driving this effect, we quantitatively measured 7596 protein targets using the high-throughput SomaScan assay and constructed a BDNF pathway 26-protein pre-to-post foldchange index (Fig. 3C, D), which was significantly upregulated (1-sample t-test: t = 3.21, p = 0.001, Cohen’s d = 0.12). While BDNF itself was not significantly affected, SLITRK1 (SLIT and NTRK-like family member 1), a protein that promotes excitatory synapse development and glutamatergic neurite outgrowth40,41, increased significantly pre-to-post (W = 39.0, p = 0.01, CLES = 0.35). NGFR (nerve growth factor receptor), a TNF receptor that binds to NGF and BDNF and plays an essential role in neural cell differentiation and survival, also increased (W = 58.0, p = 0.08, CLES = 0.35).

Metabolic reprogramming

Previous studies have characterized meditation as a hypometabolic state42 and reported enhanced glycolysis in Tibetan Buddhist monks24. To test the intervention’s effect on real-time metabolism, we treated BE(2)M17 human neuroblastoma cells with 1% plasma from all participants for 60 min and performed Seahorse XF assays. An ATP Rate assay (Fig. 4A) revealed that changes in plasma produced during the intervention induced a compensatory shift from mitochondrial to glycolytic ATP production, with post-plasma-treated cells showing significantly higher basal glycolytic rate than pre-plasma-treated cells (t = 2.95, p = 0.008, Cohen’s d = 0.36) despite no significant differences in mitochondrial (t = −1.04, p = 0.31, Cohen’s d = −0.28) or total ATP production rate (t = 1.46, p = 0.16, Cohen’s d = 0.25). A glycolytic rate assay (Fig. 4B) confirmed that, compared to pre-plasma, post-plasma increased basal glycolytic rate (t = 3.39, p = 0.003, Cohen’s d = 0.63) and lowered basal respiration rate (t = −4.80, p = 0.0001, Cohen’s d = −1.26). Mitochondrial stress parameters (Fig. 4C) did not significantly differ between cells exposed to pre- and post-plasma, while, compared to cells exposed to growth medium only, plasma-exposed cells displayed significantly higher ATP production and glycolytic rates.

Fig. 4: Metabolic effects (n = 20 participants).

A–C Seahorse XF analyzer real time cellular metabolic assays with BE(2)M17 human neuroblastoma cells treated with 1% pre- and post-retreat plasma for 60-min or assay buffer (control). Mean ± 95% CI error bars. A ATP production rate assay with total ATP production rate (mean ± 95% CI) broken into glycolytic and mitochondrial ATP production rates. B Glycolytic rate assay, including basal and compensatory glycolytic rates (top) and basal mitochondrial respiration rate (bottom). Mean ± 95% CI error bars. C Mitochondrial respiration stress assay. Mean ± 95% CI error bars. Asterisks denote * (p < 0.05), ** (p < 0.01), *** (p < 0.001), **** (p < 0.0001). D Glycolysis pathway index pre-to-post fold change. Error bars denote SEM. Index fold change levels significantly above zero signal an upregulated pathway or protein cluster. E Glycolysis pathway heat map. Fold change level per component protein for advanced and novice participants. F, G Oxidative phosphorylation pathway index with pre-to-post fold change (F) and heat map for advanced and novice participants (G). Error bars denote SEM.

As before, to see if the plasma proteome reflected these changes, we pre-selected 19 proteins involved in glycolysis and oxidative phosphorylation and calculated pre-to-post foldchanges and an index average (Fig. 4D–G). The glycolysis index increased significantly pre-to-post (t = 3.37, p = 0.0008, Cohen’s d = 0.23), with 12 upregulated targets led by ENO2 (enolase 2, a neuron-specific converter of 2-phosphoglycerate into phosphoenolpyruvate) (W = 48.0, p = 0.03, CLES = 0.33) and LDHA (lactate dehydrogenase A, converter of pyruvate into lactate in anaerobic glycolysis (W = 38.0, p = 0.01, CLES = 0.22). Oxidative phosphorylation associated proteins trended upwards (t = 1.38, p = 0.17, Cohen’s d = 0.13) but did not reach statistical significance.

Functional cellular signaling

Having investigated specific pathways of interest, we performed an exploratory, hypothesis-free analysis of proteomic and metabolomic results.

Proteomics

Volcano plot analysis (Fig. 5A) revealed 21 significantly altered proteins. Cofilin-2 (COF2) and Enoyl-CoA hydratase were significantly upregulated, which suggests enhanced cellular processes related to cytoskeletal regulation and fatty acid metabolism. IL1-F6 (interleukin-36 alpha), MYPC1 (myosin binding protein C1), LDHA, and FGF-19 (fibroblast growth factor 19) exhibited moderate upregulation.

Fig. 5: Plasma proteome (n = 20 participants).

A Plasma Proteome. Volcano plot of pre/post fold changes (x-axis) and paired t-test –log(p-values) (y-axis) for 7596 proteins. Vertical yellow dashed lines mark a linear fold change of 0.25 in either direction; horizontal yellow dashed line marks p = 0.05 cutoff. Proteins in top right/left quadrants are significantly up/downregulated. B Protein-Protein Interactions. Heat-diffused network clusters with ≤ 20 interactors for proteins with pre/post (i–iii) p < 0.05 and (iv) |foldchange| > 0.5. Top gene enrichment terms for each cluster are: (i) fatty acid metabolism, (ii) proton-transporting two-sector ATPase complex, (iii) complement system, and (iv) nucleosome organization. C Gene Enrichment. Enrichment analysis (Enrichr) with pre/post significantly altered proteins (p < 0.01) and assay targets as background. D Proteome Cluster Analysis. Hierarchical and K-means clustering of |fold changes| > 0.5 proteins. Dendrogram branch lengths indicate degree of similarity between expression profiles. Nodes represent protein targets; colors indicate fold change. E Proteomic A Priori Gene Indices. Mean pre/post fold change for protein clusters or pathways. Index proteins in Supplementary Table 15. Mean ± SEM error bars. F ELISA Assay Biomarkers. Mean ± 95% CI error bars. Asterisks denote * (p < 0.05), ** (p < 0.01), *** (p < 0.001).

Protein-protein interaction networks (Fig. 5B) revealed three significantly altered protein clusters related to mitochondrial energy production (ATP5F1, ATP6V1F); fatty acid metabolism (ECHS1, POP7, and ACAT2); nucleosome organization (HDAC1, RANGAP1, RBL2, C3, and MYOM2).

Pathway enrichment analysis (Fig. 5C) showed upregulated proteins associated with muscle cell apoptosis, amyloid precursor protein catabolism, and mitochondrial proton transport pathways (padjusted < 0.05), suggesting alterations in metabolic pathways and muscle-related cellular processes. Interestingly, pathways related to butanoate, propanoate metabolism, and tryptophan biosynthesis were enriched, suggesting shifts in key metabolic routes. (See Supplementary Figs. 3 and 4 for g:Profiler enrichment and pathway heatmaps.)

We also created proteomic indices to examine inflammation, epigenetic regulation, and short-chain fatty acid (SCFA) metabolism pathways, comparing across timepoint (pre/post) and experience level (novice/advanced) (Fig. 5E) (Supplementary Fig. 5). The epigenetic (HDAC) index demonstrated elevated sirtuins and HDAC expression in the advanced group, reflecting increased epigenetic regulation and possible stress resilience. In contrast, the novice group showed differences in mitochondrial function and chromatin remodeling. SCFA index showed increased expression of acyl-CoA dehydrogenase short-chain (ACADS) and fatty acid synthase (FASN) in the advanced group, suggesting a shift toward enhanced fatty acid oxidation and improved metabolic efficiency.

Inflammation, anti-inflammation, and cellular turnover

To assess whether the intervention elicited inflammatory or anti-inflammatory cascades, we examined a panel of 23 inflammatory and 21 anti-inflammatory proteins (Fig. 5E). We found significant upregulation of inflammatory markers (t = 3.81, p = 0.0001, Cohen’s d = 0.15), driven by increases in S100A8 (calgranulin A) (W = 30.0, p = 0.004, CLES = 0.33) and CCL2 (C-C motif chemokine 2) (W = 47.0, p = 0.03, CLES = 0.33) and trending increases in IL-6, S100A9 (S100 calcium-binding protein A9; calgranulin B), and PTGS2 (prostaglandin G/COX-2). These findings align with the role of S100A8 and S100A9 as alarmins—endogenous molecules released in response to cellular damage or stress known to induce secretion of inflammatory mediators IL-6, IL-8, and CCL243,44.

Interestingly, we also observed a significantly upregulated anti-inflammatory markers index (t = 2.25, p = 0.03, Cohen’s d = 0.09), with positively trending levels of TGF-b1 (transforming growth factor beta-1), NFKBIA (NF-kappa-B inhibitor), STAT6 (signal transducer and activator of transcription 6), CEBPB (CCAAT/enhancer-binding protein beta**)**, IL1, SOCS3 (suppressor of cytokine signaling 3), and TNFAIP3 (TNF alpha-induced protein 3). Concurrent activation of both pathways suggests a dynamic process of immune modulation, possibly reflecting enhanced cellular turnover or repair mechanisms. We also measured plasma nanoparticles and found no significant change in total nanoparticle concentration (t = −0.17, p = 0.87, Cohen’s d = −0.03), but did find a significant decrease in the percentage of particles in the exosome range (20–120 nm diameter) (t = −2.09, p = 0.04, Cohen’s d = −0.65) (Supplementary Fig. 2A, B) consistent with both higher cellular turnover and metabolic suppression.

Endogenous opioids

Placebo effects are known to engage endogenous opioid and endocrine systems. We assessed levels of 15 proteomic targets within the endogenous opioid pathway and found the pathway index to be significantly upregulated (t = 2.14, p = 0.03, Cohen’s d = 0.12), driven by positively trending increases in opioid peptide precursor PENK (proenkephalin-A), opioid peptide PDYN (dynorphin A), GABR2:CD (GABBA B receptor subunit 2: cytoplasmic domain), and CAMK2B (calcium/calmodulin-dependent protein kinase type II), which is known to increase after opioid administration (Fig. 5E). ELISA assays to confirm these findings (Fig. 5F and Supplementary Table 11) revealed significant increases of beta-endorphin (W = 14.0, p = 0.0002, Cohen’s d = 0.42) and dynorphin (W = 37.0, p = 0.009, Cohen’s d = 0.27).

Metabolomics

We performed liquid chromatography-mass spectrometry-based metabolomic analysis on plasma. Partial least squares discriminant analysis (PLS-DA) (Fig. 6A) revealed distinct metabolic profiles between time points, with key metabolites involved in synaptic plasticity, metabolism, RNA modulation, neurotransmitter availability, and inflammation and over 25 metabolites with VIP (variable importance in projection) scores > 2 (Fig. 6B, C).

Fig. 6: Metabolomics (n = 20 participants).

A PLS-DA Scores Plot. Scatter plot shows 2 principal components with greatest variation. Ovals show 95% confidence intervals. Similar observations cluster together. Component 1 (x-axis) contains 7.8% of total variation; component 2 contains 7.1%. Oval cluster spatial separation indicates systematic differences between timepoints. B PLS-DA Variable Importance in Projection (VIP) Scores. VIP scores indicate metabolite’s contribution to timepoint separation, ranked in descending order. Top metabolites are involved in neurotransmitter regulation, lipid signaling, and RNA modification, suggesting the intervention induced broad metabolic and molecular adaptations. C Top PLS-DA Compounds. Correlation coefficients between metabolites and PLS-DA discriminant function. Each metabolite’s correlation coefficient indicates strength and direction of association with timepoint separation (positive = direct, negative = inverse, closer to ± 1 = more influential in driving group separation). Top metabolites highlight changes in neurotransmitter precursors, lipid signaling, and amino acid metabolism. D Pathway Analysis. Metabolite pathways’ impact scores (x-axis) and p-values (y-axis, -log10 transformed). Larger impact values indicate greater contribution to the pathway. Bubble size reflects the pathway’s enrichment score, and bubble color corresponds to significance level (darker indicating lower p-value). E Tryptophan Pathway. Colored boxes represent detected metabolites. Color scale represents magnitude of pre-to-post change. Box plots show pre and post concentration distributions and paired t-test values with raw p-values. F Important Features Heat Map showing top differentially expressed metabolites, including several phenylalanine-related metabolites involved in neurotransmitter synthesis and metabolic regulation. Color scale indicates relative abundance (red = higher; blue = lower). Clustering highlights metabolic pathway shifts, including phenylalanine metabolism and lipid signaling.

MetaboAnalyst pathway enrichment detected 53 impacted pathways (Fig. 6D and Supplementary Table 16). Six showed significant pre-to-post changes (p < 0.05), led by tryptophan metabolism (pFDR = 0.03) (Fig. 6E and Supplementary Fig. 6), with decreases in upstream and downstream metabolites, including L-tryptophan (p = 0.03), tryptamine (p = 0.04), L-kynurenine (p = 0.04), indole-3-acetate (p = 0.001), and 5-methoxyindoleacetate (p = 0.001). The steroid hormone biosynthesis pathway showed lower androstenediol (p = 0.04), testosterone (p = 0.19), cortisol (p = 0.06), cortisone (p = 0.07), 21-deoxycortisol (p = 0.77), corticosterone (p = 0.07), and 11-dehydrocorticosterone (p = 0.02). Lower cortisol levels reported in previous meditation studies were interpreted as signs of improved stress response regulation45.

Among top-ranking metabolites, three were related to phenylalanine metabolism, including CAN-2-Hydroxy-3-phenylpropionic acid, N-acetyl-D-phenylalanine, and phenethylamine, which can act as neurotransmitter storage, contribute to stress adaptation, and enhance dopamine, norepinephrine, and serotonin release. Two RNA-related metabolites, 1-methyladenosine and N6-methyladenosine, were linked to inflammatory responses, with roles in RNA metabolism, methylation, and cellular signaling. 3-indolebutyric acid, associated with tryptophan metabolism, may influence neurotransmitter synthesis and synaptic plasticity. Several dipeptides were identified as enhancers of gut microbial activity, promoting SCFA production, which is crucial for regulating inflammation and immune function. S1P (d18:1), a signaling lipid, may reflect gut barrier integrity and inflammatory pathway alterations. Figure 6F shows a distinct clustering of metabolites by relative abundance, highlighting these key contributors.

Exosome-specific transcriptomics

We analyzed differentially expressed exosome-specific extracellular microRNAs, non-coding RNAs, and RNAs mapping to protein-coding mRNAs. Data from 6 participants was excluded at preprocessing (n = 16). At least 18 non-coding exRNAs (p < 0.1, log2FC > ±0.58) exhibited distinct expression profiles on pre and post timepoints (Fig. 7A, B). Principal components (PC) explain 46.3% of the total timepoint variation (Fig. 7C), indicating significant shifts in non-coding exRNA expression during the intervention. At least 5.99% of the variance was attributed to experience levels, with partial separation between novice and advanced meditators (Fig. 7D). Correlations between principal components and key variables such as timepoint and experience reveal that PC1 and PC2 strongly correlate with timepoint (R² = 0.35 and 0.28, respectively), while PC9 moderately correlates with experience (R² = 0.31), indicating that both influence exRNA expression (Fig. 7E).

Fig. 7: Exosome transcriptomics (n = 16 participants).

A Differential non-coding exRNAs (miRNAs, ncRNAs) on volcano plot of -log10p versus log2FC (pre/post difference). Blue/red dots denote exRNAs prevalent in pre/post-intervention exosomes, respectively; green dots denote exRNAs without significant difference. Dashed orange lines indicate p < 0.1 and log2FC > ±0.58 thresholds. Full list of exRNAs on Supplementary Table 17. B Differential expression (log2-transformed TMM-normalized counts) of 18 non-coding exRNAs. Colors indicate counts. Sample name: (M male, F female)_participantID_(pre; post)_(E, advanced; N, novice). C, D Principal component analysis of top normalized exRNA counts. Encircled are clusters corresponding to C pre/post, and D novice/advanced. E Eigenvector plot showing the correlation of principal components to variables’ metadata and test significances (** denotes p < 0.01; * denotes p < 0.05). Values correspond to Pearson R2 values and are colored by significance. F Differential protein-coding exRNAs (excluding ncRNAs) on volcano scatter plot of -log10p versus log2FC (pre/post difference). Blue/red dots denote exRNAs prevalent in pre/post exosomes, respectively; green dots denote exRNAs without significant difference. Dashed orange lines indicate p < 0.1, log2FC > ±0.58 thresholds. G TMM-normalized counts of log2-transformed protein-coding exRNAs with highest variance across samples. Color scale indicates counts. Sample name: (M male, F female)_participantID_(pre; post)_(E, advanced; N, novice). H Reactome pathway analysis of protein-coding exRNAs enriched in pre-vs-post intervention exosomes. Color indicates p-value; circle size indicates exRNA count per pathway.

We also identified exRNAs mapped to at least 66 annotated protein-coding mRNAs upregulated or downregulated post-intervention (Fig. 7F, G). Some of the upregulated exRNAs could be mapped to ras/rab interactor 1 (RIBC1), synapsin 3 (SYN3), and glutamate ionotropic receptor kainite type subunit 3 (GRIK3), genes linked to synaptic function and neurotransmission. The enrichment of exRNAs specific to solute carrier family 27 member 1 (SLC27A1) and proprotein convertase subtilisin/kexin type 9 (PCSK9) could enhance neural signaling and metabolic regulation. Downstream functional analysis of exRNAs in pre-vs-post exosome fractions against the Reactome database46 (Fig. 7H) predicted pathways related to neurotransmission (Serotonin/Dopamine Neurotransmitter Release Cycle), further highlighting the enhancement of synaptic activity post-intervention. Metabolic pathways, including Transport of Vitamins and Nucleosides, could reflect broader metabolic changes.

Machine learning and feature-MEQ correlation

To identify features differentiating timepoints and experience levels across metabolomic, proteomic, fMRI, and RNA outcomes, we evaluated the data using Random Forest and XGBoost models. After normalization and dimensionality reduction, performance was assessed with tenfold cross-validation using F1 score, precision, recall, and area under the receiver operating characteristic (AUROC). Both models demonstrated robust classification across conditions, with AUC values ranging from 0.70 to 0.93. SHAP analyses revealed key contributors spanning metabolites, fMRI connectivity features, proteins, and non-coding RNAs (Fig. 8A, B).

Fig. 8: Machine learning and MEQ-features correlations (n = 20 participants).

A ROC curves and SHAP plots for XGBoost and Random Forest models predicting pre/post classification: AUC = 0.86 (XGBoost) and 0.90 (Random Forest) indicate good classification performance. SHAP plots display the top contributing features ranked by impact on model output. B ROC curves and SHAP plots for XGBoost and Random Forest models predicting novice/advanced classification: AUC = 0.70 (XGBoost) and 0.93 (Random Forest). Abbreviations: R (right), L (left), AN (auditory network), SN (salience network), DMN (default mode network), SMN (somatomotor network), ECN (executive control network), VN (visual network), DAN (dorsal attention network), dlPFC (dorsolateral prefrontal cortex), AI (anterior insula), AG (angular gyrus), GAPDH (glyceraldehyde-3-phosphate dehydrogenase), ATP5PB (ATP synthase peripheral stalk membrane subunit B), FGF (fibroblast growth factor), NTRK2 (neurotrophin receptor tyrosine kinase 2). C Heatmaps of Spearman R correlations between Mystical Experience Questionnaire (MEQ) scores and the top machine learning features per model and time point, with * denoting FDR-adjusted statistical significance (pFDR < 0.05).

Focusing first on