Psilocybin has rapidly emerged as the leading therapeutic candidate amongst the new wave of psychedelic medicine, particularly within the context of major depressive disorder [1, 2]. Despite its pre-regulatory status, the drug is widely accessible in organic-form due to the recent legislative initiatives of several countries and US-states, giving rise to self-medication practices and the emergence of retreat-style guided psilocybin experiences [3] to fill this need. Although the majority have been marketed towards personal-growth and wellness, one of the most frequently cited motivations for undertaking self-guided or retreat-style guided psilocybin experiences is self-medication for pre-existing depressive symptoms, amongst other physical and mental complaints [4,5,6]. Despite their flourishing popularity amongst those seeking alternative routes to improved depressive symptoms [7] very few empirical studies have been conducted amongst individuals participating in retreat-style psychedelic experiences.

Notwithstanding the magnitude and longevity of depressive symptom improvements observed in psychedelic clinical trials, up to 30% of participants in research studies fail to demonstrate adequate response [2], suggesting that person-level or environmental factors may determine likelihood of deriving a clinically-meaningful benefit. This may be explained, in part, by the heterogenous nature of depressive symptoms, which span several interconnected yet highly divergent constructs, such as mood, appetite, cognition, and sleep. Perhaps the most scarcely described of these factors amongst psychedelic research is sleep, and currently very little is known about whether psilocybin alters sleep, particularly beyond acute states of drug metabolism window. Only one small study [8] has examined the effects of non-therapeutic psilocybin administration on sleep-EEG measures amongst 19 healthy humans, reporting prolonged REM latency following a single daytime administration session. The same study also reported a reduction in delta/slow-wave power during SWS in the first sleep cycle, suggesting alterations in homeostatic sleep drive. However, these sleep-EEG recordings were observed the night immediately after psilocybin administration, and therefore the long-term effects of psilocybin on sleep remain entirely unknown. This is pertinent, as the active metabolites and transient alterations in receptor trafficking [9] which are present during the acute psilocybin administration period are likely to influence sleep in ways that are different to any neuroplastic changes which may be present several weeks after the drug has been metabolized. Similar findings of prolonged REM latency have emerged from animal models [10], following administration of psilocin (the pharmacologically active metabolite of psilocybin) in mice.

However, in animal studies enhanced oscillations in a narrow window of the delta (0.5-4 Hz) frequency band, around 4 Hz were observed, which contrasts with the reduction in 0.5-4 Hz power observed in humans[8]. To our knowledge, no human-subjects study has investigated whether improvements in self-report sleep quality are observed following Administration Of Psilocybin-based Therapeutics (ADOPT), or whether these changes are sustained beyond the period of acute administration. Using naturalistic survey data obtained from participants who took part in facilitated psilocybin sessions within retreat settings, or guided psilocybin experiences, we aimed to examine the impact of psilocybin on sleep and other depressive symptoms at 2-weeks and 4-weeks post-intervention, and the possible impact of baseline sleep quality on psilocybin-mediated antidepressant response. Indeed, given sleep’s impact on many facets of emotional wellbeing, including mood and anxiety [11], poor sleep quality has considerable potential to contribute to the physiological ‘set’ of an individual’s psychedelic experiences, which may contribute to their overall long-term benefits. We explored two key questions; 1) Does guided psilocybin administration lead to sustained improvements [up to 1 month] in sleep and depressive symptoms. 2) Do baseline sleep disturbances, and or changes in sleep symptoms predict improvement in depressive symptoms following psilocybin administration.

Data were obtained from participants (n = 886) who identified an intent to participate in psychedelic use in the near future (for detailed methods see [12]). Participants were required to be 18 years or older, have sufficient English comprehension, and must have indicated that they were about to participate in a retreat, ceremony or other guided experience involving the use of a classic psychedelics containing 5-HT2A receptor agonist (e.g., Psilocybin, DMT, Mescaline, or LSD). Data were collected prospectively from the participants using online surveys completed before, during, and for up to four weeks after their psychedelic experience. We restricted our analyses to respondents who reported psilocybin use, given that it has the greatest clinical relevance to depression (n = 653). Ethical approval was granted by the joint research compliance office and the Imperial College Research Ethics Committee (ICREC reference 181C4346). Data originated from 111 centers, some of which contributed substantially to the proportion of the sample (max = 297 participants), whereas others provided only a single or several participants.

We measured depressive symptom severity using the Quick Inventory for Depressive Symptoms (QIDS [13]), a well-validated 16 item scale which covers the key diagnostic domains of depressive symptomology. The QIDS includes a sleep sub-scale, comprised of four items assessing sleep onset insomnia (Item 1 ‘Difficulty falling asleep’), sleep maintenance insomnia (Item 2 ‘Sleep during the night’), early morning insomnia (Item 3 ‘Waking up too early’), and hypersomnia (Item 4 ‘Sleeping too much’). In keeping with standard scoring instructions for the QIDS, we used the maximum value of the individual item scores (each scored 0 [no disturbance] to 3 [max disturbance]) of these four items to define the severity of sleep disturbance symptoms (max score 27). Although the QIDS was designed for the assessment of global symptom severity, previous studies have used the QIDS-sleep sub-scale in isolation to quantify sleep symptom improvements in other rapid acting antidepressants, as well as the relationship between depression related sleep-disturbance and functional outcomes [14,15,16,17]. We used the QIDS total score to measure depressive symptom severity (QIDS-depression). Naturally, as the score from the sleep sub-scale contributes to the total QIDS total score, a significant relationship between QIDS-sleep and QIDS-total is implicit. Therefore, we excluded the sleep scale from the total QIDS score (i.e., QIDS-depression = QIDS-total – QIDS-sleep). We defined depressive symptom remission at both two-week and four-weeks as a QIDS-depression reduction > 50% relative to baseline, and QIDS-total score < 6 at the respective timepoint.

Sleep and depressive symptoms share a highly complex and bi-directional relationship. Therefore, to better aid in parsing these pathways, we used a nested sequence of a priori and confirmatory post hoc models that we further outline below. We first examined univariate change in sleep and depressive symptoms using random-effects models. Next, we used structural equation models to examine contemporaneous and longitudinal interactions between sleep and depressive symptoms at baseline, two-weeks and four-weeks. Finally, we probed symptom-level effects and predictors of treatment response using a novel network analysis approach, adapted previously for use in sleep and mental-health symptom networks [18].

Outcomes were change in QIDS-depression total score and QIDS-sleep total score at two weeks and four weeks. Fixed effects in both models included timepoint (baseline, two-weeks, four-weeks), and the respective outcome measure at baseline to account for potential confounding of initial symptom severity. We performed model diagnostics using visual assessments of the normality and heteroskedasticity of residuals (see supplement for further information on diagnostics). We interpreted unadjusted models as our primary outcome, to optimise model parsimony and facilitate hypothesis generation, and to aid in comparison with future studies, given this early exploratory stage. However, we conducted pre-planned sensitivity analyses to conduct a secondary corrected model by incrementally adjusting for available covariates (Age, ethnicity, marital status, employment, income, country of retreat, number of sessions, psilocybin dose, retreat-center) exploring their impact on model fit and interpretation, which we report in full in the supplement. To limit overfitting and ensure appropriate variance component estimates, we based our model selection on both the Restricted Maximum Likelihood (REML) and Bayesian Information Criterion (BIC) values. REML was used to provide unbiased estimates of variance components, while BIC helped identify the most parsimonious model by balancing fit and complexity.

Given that large variations in sample size can give rise to spurious fixed effects in nested random-effects models, we also performed a Bayesian-simulation using Montecarlo Markov Chains to generate robust estimates of expected variation in our two outcome variables (sleep and depressive symptoms) across the range of center sample sizes, to distinguish any meaningful center-based effects from natural variance. Bootstrapping was performed on the uncorrected, and best-fitting corrected models to obtain estimates of the stability of the 95% Bootstrapped Confidence Intervals (BCI) of the fixed-effect coefficients. Missing outcome data at each two weeks (53%) and four-week (55%) were determined to be Missing Not at Random (MNAR), as greater observed values for: age, number of sessions, duration of retreat, and outcomes at baseline were significantly associated with increased missingness. Therefore, we performed imputation sensitivity analyses using Baseline Observation Carried Forward (BOCF) to simulate how a ‘worse-case’ scenario might bias our findings, in the event that that those with missing values demonstrated no improvement. We also used a deep-learning imputation approach to impute values based on auxiliary variables predictive of outcome (see supplement for methods).

Structural equation modeling was used to examine the dynamic relationships between QIDS-depression (represented as qids in formulae notation – without sleep items) and sleep quality [QIDS-sleep] across three time points. A cross-lagged panel framework was used to model the variance and covariance (Σ) where βij​ and γij​ represent the coefficients for the QIDS-depression and QIDS-sleep variables, capturing the lagged events across different time points. Error terms ε-qids​ and ε-sleep​ account for unobserved factors influencing depressive symptoms (qids [without sleep]) and sleep quality (sleep) at each time point. This approach allowed us to examine how each variable influences the other across time and to probe both the temporal order effects and potentially causal dynamics between sleep and qids variables measured at multiple time points. Path coefficients were estimated using Maximum Likelihood Estimation and reported as standardized λ values. Lagged analyses controlled for the baseline value of the lagged variable (i.e., Sleep0), as well as the value of the covariate (i.e., Depression2) at the same time point (i.e., Depression3 ~ Depression2 + Sleep0 + Sleep2). Diagnostics and methods are reported in the supplement.

To investigate symptom-specific relationships between individual QIDS items and overall depressive-symptom improvement, we adapted a Network Interventional Analysis [18] and estimated, a Mixed Graphical Model [19] at baseline where each node is represented by one of the 16 QIDS items. We then estimated their conditional relationship with remission likelihood by introducing a remission node, which represented whether or not participants experienced a remission at two weeks (network 1) or four weeks (network 2).

Network analyses were used to generate mixed graphical models [20] via K-degree nodewise regression, which estimates the conditional dependencies between all nodes in the network. LASSO regularization was used to calculate the adjusted edge weights of nodewise regressions. In accordance with previous recommendations [21, 22]. We generated 1000 new resampled networks using the ‘resample’ function of the MGM package, and verified the accuracy and stability of the estimated network models [23] by assessing the probability (%) that a given edge (i.e., link between symptoms) was present in the network after resampling, as well as the strength of the included links, by assessing the range of the adjusted edge weights observed during resampling. To facilitate comparisons between edges, we performed post-hoc edge difference tests ( Bootnet::differenceTest() in R), to compare the magnitude of differences between observed edges.

In our network analyses, all variables (e.g., the QIDS items and the response variables [remission]) are included as individual nodes and visualised in a network. Nodes are linked by edges that represent conditional-dependence among them, i.e., the unique association between two variables after conditioning on all other variables in the network.. Therefore, the relationships between the individual symptoms and remission we report implicitly controls for the confounds of the remaining symptoms. We used the magnitude and path of adjusted edge weights to interpret the magnitude of these relationships, whereby positive (Greater likelihood of remission) and negative (Smaller likelihood of remission) relationships were represented by blue (positive) and red (negative) edges.

In order to quantify remission cases, we confined the sample used for the network analysis to those who demonstrated clinically significant depressive symptoms at baseline (n = 392), as defined as a score > 5 on the QIDS total score, a common cutoff for ‘mild depressive symptoms’[13]. For transparency, we also report networks derived from the whole sample in the supplement. N.B: The multilevel models and SEMs were first performed using the entire sample, followed by a sub-group analysis with clinically significant symptoms. We further probed the extent to which any relationships between sleep symptoms and depressive symptom improvement may have been influenced by the presence of sleep items in the outcome variable, after reintroduction of sleep item. We conducted primary analyses uncorrected but performed sensitivity analyses by sequentially introducing additional nodes for study center, psilocybin dose, and the number of sessions to examine their impact on the network relationships demonstrated in our primary analyses (see supplement).