Using large-scale brain modelling and dynamic sensitivity analysis of resting-state data in responders versus non-responders, the study identified specific brain regions whose perturbation predicts a transition from depressive to healthy brain dynamics following psilocybin therapy. These regions co-localise with in vivo 5-HT2A and 5-HT1A receptor density maps, providing causal mechanistic evidence that serotonergic targets mediate psilocybin's therapeutic effect in treatment-resistant depression.
Psilocybin therapy for depression has started to show promise, yet the underlying causal mechanisms are not currently known. Here, we leveraged the differential outcome in responders and non-responders to psilocybin (10 and 25 mg, 7 days apart) therapy for depression—to gain new insights into regions and networks implicated in the restoration of healthy brain dynamics. We used large-scale brain modelling to fit the spatiotemporal brain dynamics at rest in both responders and non-responders before treatment. Dynamic sensitivity analysis of systematic perturbation of these models enabled us to identify specific brain regions implicated in a transition from a depressive brain state to a healthy one. Binarizing the sample into treatment responders (>50% reduction in depressive symptoms) versus non-responders enabled us to identify a subset of regions implicated in this change. Interestingly, these regions correlate with in vivo density maps of serotonin receptors 5-hydroxytryptamine 2a and 5-hydroxytryptamine 1a, which psilocin, the active metabolite of psilocybin, has an appreciable affinity for, and where it acts as a full-to-partial agonist. Serotonergic transmission has long been associated with depression, and our findings provide causal mechanistic evidence for the role of brain regions in the recovery from depression via psilocybin.
Papers cited by this study that are also in Blossom
Kringelbach, M. L., Cruzat, J., Cabral, J. et al. · PNAS (2020)
Vohryzek and colleagues frame their study around a broader challenge in systems neuroscience: although depression is associated with altered brain dynamics, it remains difficult to identify which dynamic features are most relevant to symptoms and treatment response. The paper notes that major depressive disorder is often linked to rigid, self-referential thinking and to abnormalities in resting-state organisation, including default mode network overactivity and reduced control network function. At the same time, prior work has suggested that psychedelic therapy may increase the flexibility and repertoire of brain states, but the mechanisms underlying clinical response remain uncertain. The authors therefore aim to use empirical fMRI data and large-scale computational brain modelling to identify brain dynamics predictive of response to psilocybin in treatment-resistant depression. More specifically, they ask which spatio-temporal brain substates distinguish responders from non-responders, whether a model fitted to pre-treatment data can reproduce these dynamics, and which brain regions appear most important for driving a transition towards a healthier post-treatment state. They also hypothesise that regions involved in successful transition will relate to the distribution of 5-HT2a and 5-HT1a serotonin receptors. The paper presents an exploratory modelling analysis built on data from a clinical psilocybin study and uses this as a test case for dynamic sensitivity analysis, with the broader implication that such models might eventually support mechanistic and personalised psychiatry.
The researchers analysed previously published fMRI data from patients with treatment-resistant major depression treated with psilocybin at Imperial College London. Of the original 19 patients, 15 were included in the present analysis after excluding scans with excessive movement and other artefacts. Participants underwent MRI before treatment and 1 day after the second psilocybin dosing session. The treatment consisted of two oral doses, 10 mg and 25 mg, given 7 days apart. Clinical response was defined at 5 weeks post-treatment using the Quick Inventory of Depressive Symptomatology (QIDS), with response defined as more than 50% symptom reduction from baseline; 6 of the 15 patients met this criterion. To characterise brain dynamics, the authors used Leading Eigenvector Dynamics Analysis (LEiDA) to identify probabilistic metastable substates (PMSs) from the fMRI time series. They computed phase coherence across 90 brain regions and clustered the dominant phase-locking patterns using k-means across cluster solutions from k=2 to k=10. A three-cluster solution was selected as the best balance of clustering quality metrics and statistical separation between groups. The probability of occurrence of each PMS was then calculated for each recording session. The modelling component used a large-scale brain network model in which each of the 90 regions was represented by a Hopf bifurcation oscillator, coupled according to a structural connectivity template derived from diffusion tensor imaging in 16 healthy young adults. The model parameters were fitted separately for responder and non-responder pre-treatment data by matching simulated and empirical PMS probabilities. The authors used the global coupling parameter G to optimise model fit, and also derived intrinsic frequencies from the empirical fMRI spectra. To probe causal sensitivity, the researchers performed dynamic sensitivity analysis by systematically perturbing regional dynamics in silico. This was done by varying the bifurcation parameter a for homological pairs of regions, effectively moving local dynamics between more noise-driven and more oscillatory regimes. They then assessed how well each perturbation moved the simulated system towards either a healthy target state, defined by the responder post-treatment PMS distribution, or a depressive target state, defined by the non-responder post-treatment PMS distribution. The main fit metric was the symmetrised Kullback-Leibler divergence between empirical and simulated PMS probabilities. Statistical testing used 1,000 non-parametric permutations for experimental analyses, and Pearson correlation for receptor-map associations. Finally, the authors compared the regional perturbation effects with previously published PET-based maps of serotonin receptor density, focusing on 5-HT2a, 5-HT1a, 5-HT2b, 5-HT4 and the serotonin transporter (5-HTT).
Using LEiDA, the authors identified three recurrent PMSs and compared their probabilities before and after psilocybin treatment in responders and non-responders. The occurrence of substate 3 differed significantly in responders when comparing pre- and post-treatment data (P = 0.0258), and also differed in the post-treatment comparison alone (P = 0.0141). In contrast, global brain connectivity, metastability and related summary measures did not show significant differences, suggesting that the relevant treatment effects were better captured by spatio-temporal state dynamics than by static or global measures. Functional connectivity dynamics showed additional group differences. The authors report significant differences in temporal similarity of spatial patterns between pre- and post-treatment responders (P = 0.0163) and also in the comparisons involving non-responders and post-treatment responders (P = 0.0183 and P = 0.0273). These findings were taken to indicate that temporal reorganisation of large-scale brain activity was a key feature of response. The large-scale brain model was fitted successfully to the empirical PMS distributions of both pre-treatment groups. The optimal global coupling values were G = 0.185 for responders and G = 0.165 for non-responders, with corresponding Kullback-Leibler divergences of 0.0064 and 0.0187. This suggests that the model reproduced the observed spatio-temporal dynamics reasonably well, and more closely for the responder group. Dynamic sensitivity analysis showed that perturbations in the noise-driven regime (a < 0) worsened the match to the healthy target state, whereas more oscillatory perturbations (a > 0) initially improved the match before becoming less effective at stronger values. When the target was switched to the depressive state, the best fit was obtained with no perturbation, and increasing oscillatory drive or noise did not improve the fit. The optimal perturbation intensity for moving towards the healthy state was a = 0.07. At that stimulation intensity, the regions most associated with promoting transition towards the healthy state in responders, but not non-responders, were the temporal superior pole, rolandic operculum, fusiform gyrus, supplementary motor area, inferior parietal gyrus, angular gyrus, supramarginal gyrus, opercular inferior frontal gyrus, orbital middle frontal gyrus and parahippocampal gyrus. The authors interpret these as regions implicated in shifting away from a depressed brain configuration towards the post-treatment responder pattern. The regional perturbation effects were positively correlated with serotonin receptor maps. Differences in perturbation effects between responders and non-responders correlated with 5-HT2a receptor density (Spearman’s ρ = 0.227, P = 0.032) and 5-HT1a receptor density (ρ = 0.284, P = 0.007). Correlations with 5-HT2b, 5-HT4 and 5-HTT were not significant or were only weakly suggestive in the case of 5-HT2b (ρ = 0.064, P = 0.055).
The authors argue that their combined empirical and in silico approach identified a subset of brain regions that may be important for recovery from treatment-resistant depression after psilocybin therapy. They interpret the regional perturbation findings as evidence that these areas support a transition from a depressed brain state to the healthier configuration seen in treatment responders. The overlap with 5-HT2a and 5-HT1a receptor density is presented as biologically plausible, because psilocin has affinity for these receptors and 5-HT2a agonism may activate signalling pathways relevant to antidepressant effects. They also relate their results to the wider literature on depression, noting that the regions implicated in the transition may be connected to default mode network dysfunction and serotonergic projections from raphe nuclei. In their view, the findings are consistent with the idea that psychedelic therapy can modulate large-scale brain dynamics, potentially allowing a more flexible state that supports clinical improvement. A further point in the Discussion is methodological. The authors present probabilistic metastable substates as a way of describing complex, non-stationary brain dynamics in a relatively simple but interpretable framework. They suggest this approach could be extended to other neuroimaging modalities such as EEG and MEG, and potentially to more fine-grained spatial data. They also position their in silico perturbation framework as a useful way to explore how interventions might change brain-wide dynamics without relying solely on direct experimental manipulation. The authors acknowledge several limitations. The use of regional parcellation and k-means clustering introduces methodological choices that may affect the substates identified. The clinical dataset was an uncontrolled open-label feasibility pilot study with a small sample size and no placebo group, so they describe the work as exploratory and call for replication. They also note that the healthy target state was defined using the 1-day post-treatment scan, whereas response classification was based on the 5-week outcome, which may not fully align. Finally, the brain models were built from group-average functional and structural data rather than individual-specific models, and the authors state that future work should move towards personalised modelling for clinical application.
We carried out the analysis on previously published data set of patients with treatment-resistant depression undergoing treatment with psilocybin at Imperial College London.In brief, we investigated 15 patients (without excessive movement and other artefacts from the original 19 patients) who were diagnosed with treatment-resistant major depression. The MRI scanning sessions were completed pretreatment with psilocybin and 1-day post-treatment with the treatment consisting of two oral doses of psilocybin (10 and 25 mg, 7 days apart). The patients were split into responders and non-responders to the treatment based on the Quick Inventory of Depressive Symptomatology (QIDS) at 5 weeks post-treatment with 6 out of the 15 patients meeting criteria for response.
In this study, white matter (structural) connectivity of 90 automated anatomical atlas brain areas from a previously obtained data set was used for the large-scale brain network model. In brief, the group consisted of 16 healthy young adults (5 females, mean SD age: 24.7 ± 2.54). Diffusion tensor imaging was applied following the methodology described in Cabral et al.Undirected structural connectivity (SC) C np was obtained were n and p are brain areas and the connectivity weights are defined as the proportion of sampled fibres in all voxels in region n that reach any voxel in region p. Finally, the individual structural connectomes were averaged across the 16 subjects to obtain a group-based template.
Firstly, we calculated the instantaneous phased relationship between individual brain regions by expressing the demeaned regional fMRI signal x(t) as an analytical signal, i.e. in terms of its time-varying phase θ(t) and amplitude A(t) as x(t) = A(t) * cos(θ(t)).We excluded the first and last three timepoints to account for the boundary artefacts introduced by the Hilbert transform. Hence, for every timepoint t and pair of brain regions n and m, we obtain the phase coherence matrix dPC as follows: ( By decomposing the signal in this way, we can look at when the brain regions n and m are aligned with similar angles, cos (0) = 1, orthogonal to each other, cos (π/2) = 0 , and anti-aligned, cos (π) = -1. As the phase coherence is a measure of undirected connectivity, the phase coherence matrix dPC is symmetric and all the meaningful information is captured in the upper triangular matrix. For further analysis, we used only the 1xN leading eigenvector V 1 (t) of the dPC matrix as described in the Leading Eigenvector Dynamics Analysis (LEiDA).In detail, at every timepoint t of the dPC(t), we performed the eigendecomposition taking the first (most dominant) eigenvector to describe the dPC(t) pattern. The dPC(t) is decomposed as dPC(t) = V(t)D(t)V -1 (t), where D is the diagonal matrix carrying the real-valued eigenvalues and V 1 (t) and V -1 1 (t) are the left and right corresponding orthogonal eigenvectors, respectively. The dominant connectivity pattern can be simply reconstructed by the following matrix multiplication: The next step was introduced to find recurrent phaselocking substates in the dominant connectivity patterns varying in time. To that end, we clustered all the leading eigenvectors obtained from all the fMRI scans obtained from both responders and non-responders, thus achieving common phase-locking substates for both groups. For the clustering, we used the unsupervised k-means algorithm, with varying cluster numbers k from 2 to 10 clusters. The algorithm was run in an n-dimensional space where n = 90 (number of automated anatomical atlas brain regions) with condition × subject × timepoints number of observations. For each run, we randomly initialized the algorithm 20 times to ensure stability in the clustering. Again, by computing the matrix multiplication of the 1xN cluster centroids V cα as V cα (t)V T cα (t), we obtain the dominant connectivity pattern of each cluster. In the current analysis, we considered the cluster solution k = 3 as an optimal choice between the quality measures-Dunns, Davies-Bouldin and silhouette score, and Davies (Supplementary Fig.), and the maximizing of the statistical significance between patient groups (P-values). After calculating the phase-locking states, we defined the probability of occurrence of the individual substates by simply dividing their occurrence in each recording session by the total number of timepoints recorded (same for all recordings).
In order to simulate the ultra-slow fluctuations in fMRI signal detected during rest, we used the Landau-Stuart oscillator canonical model, describing the transition from a noisy to an oscillatory dynamics.The so-called supercritical Hopf bifurcation model was used locally at every brain region (node) to emulate the local dynamics.To achieve a large-scale brain-level description, the individual Hopf models were coupled in a SC network, describing the large-scale white matter map of the human brain.The emerging and complex interactions in the large-scale brain network of coupled Hopf models have been shown to describe many aspects known from experimental recordings in MEGand fMRI.Formally, the normal form of the supercritical Hopf bifurcation model for a single uncoupled region of interest (n) in Cartesian coordinates is described by the following set of coupled equations: with βη n (t) being the Gaussian noise with SD of β = 0.02. The bifurcation parameter a positions the system at the supercritical bifurcation point when a = 0, noise activity governed by βη n (t) in regime when a < 0 and stable limit cycle with oscillatory behaviour of frequency defined by f n = ω n /2π when a > 0. The values of the intrinsic frequency ω were calculated from the experimental fMRI signals in the 0.04 -0.07 Hz band by taking the peak frequency of the Gaussiansmoothed power spectrum of each brain area. To describe the coupled large-scale brain computational model, we introduced the coupling term (modelled as the common difference coupling, i.e. describing the linear term of a general coupling function) between the individual nodes weighted by the corresponding values of the SC matrix. To be noted, we do not consider the next non-linear coupling term following Taylor expansion of the full coupling, in case the linear coupling is non-existent.Equations () and () can be hence expanded as follows: where C np is the SC weight between node n and p, and G is the global coupling weight with equal contribution between all the nodal pairs. The SC matrix was rescaled to have the mean value 〈C〉 = 0.2 in order to be consistent with previous literature's range of parameters.The simulated signal is described by the x n equation for every node n. The variables G and a are the control parameters used for the model fitting to the experimental data and the stimulation protocol, respectively.
In order to validate the simulated signal, different realizations of the experimental data can be used.The most standard approach is comparison of the simulated data with grand-averaged static functional connectivity as computed by the Pearson correlationor metastability defined as the SD of the Kuramoto order parameter (Supplementary Fig.-Metastability). To account for the temporally varying nature of the blood oxygen level-dependent signal, recent literature has focused on the comparison between the simulated and empirical functional connectivity dynamics (FCD) spectrums (quantified by the Kolmogorov-Smirnov distance), i.e. the distributions of the cosine distance between the consecutive timepoints as described by the leading eigenvector(Supplementary Fig.-Functional Connectivity Dynamics). As alluded to in the previous section, the fMRI signals organize into spatially meaningful phase-locking states. Here, we compare the simulated data to the probabilities of occurrence of the phase-locking states found in the experimental recordings.We used the symmetrized Kullback-Leibler divergence (KL divergence) of the simulated and empirical probabilities of occurrence as follows: KL(P emp , P sim ) = 0.5 with P emp and P sim being the empirical and simulated probabilities of occurrence of the same phase-locking states, respectively.
We use non-parametric permutation tests for the experimental analysis with 1000 permutations. For the receptor analysis, we used standard Pearson's correlation. Both tests were implemented in MATLAB.
In summary, a quantitative characterization of the spatiotemporal dynamics recorded with fMRI was obtained using LEiDA, resulting in the definition of PMSs (Fig.), whose probability of occurrence was compared across conditions (i.e. within-subject design-therefore, before versus after treatment). We then constructed two large-scale brain models representative of the pre-treatment brains to psilocybin therapy. This was done by fitting their PMS descriptions to those obtained from the experimental data (Fig.). Finally, a dynamic sensitivity analysis was implemented to both responder and non-responder pre-treatment models to identify the brain regions that permit a transition to the healthy PMS [described by responders' (as predicted by the 5-week end-point) 1-day post-treatment brains; Fig.and]. As described in the 'Materials and methods' section, we computed the PMS pre-and post-treatment with psilocybin (where 'post' = 1 day post psilocybin dosing session two), for both responders and non-responders (determined 5 weeks hence). Here, we focused on a three-substate solution-the lowest k-level with statistically significant differences between the two groups as well as optimal quality measures across clustering solutions (Supplementary Fig.). When contrasting responders versus non-responders, the occurrence of substate 3 was significantly different pre-versus post-treatment (P = 0.0258, signed rank-sum test), as well as in the post-treatment data alone (P = 0.0141, rank-sum test; Fig.). Furthermore, we also computed the global brain connectivity, metastability and FCD measures (Supplementary Fig.). These results clearly indicated the necessity of considering both spatial and temporal dimensions to differentiate between conditions as global brain connectivity, synchrony and metastability show non-significant results. Conversely, the FCD measure showed significant differences in the temporal similarities of spatial patterns between pre-and post-treatment responders (P = 0.0163, signed rank-sum permutation test) and pre-and posttreatment non-responders with post-treatment responders, respectively (P = 0.0183 and P = 0.0273, rank-sum permutation test), further supporting the use of spatio-temporal measures to capture the alterations in large-scale brain dynamics across conditions. To obtain large-scale brain computational models representative of the two groups of patients (responders and non-responders before treatment), we first defined a generalized brain network model, where each of 90 cortical and subcortical brain regions (defined using automated anatomical labelling) was described by a Stuart-Landau oscillator (see the 'Materials and methods' section), and regions were coupled according to realistic SC obtained from diffusion MRI. To adjust the model to each group of patients, first, the intrinsic frequency of each brain region was set to the peak frequency in fMRI signals averaged across patients in the same group (Supplementary Fig.). Subsequently, the global coupling parameter, G, was tuned to optimize each model to its appropriate working point. This was achieved by minimizing the divergence between the experimental and simulated PMS spaces (see Fig.). In Supplementary Fig., we report optimization curves for other observables such as the static FC, metastability and FCD. For the responders and non-responders before treatment, we found G = 0.185 (KL divergence = 0.0064) and G = 0.165 (KL divergence = 0.0187), respectively, to minimize the difference. Figureshows, on the left, the experimental results for both groups before treatment; in the middle, the optimal simulated fits for both groups; and on the right, the experimental results after treatment (with the results of responders after treatment serving as the target PMS for rebalancing). These findings clearly demonstrate the large-scale model's ability to fit the spatio-temporal dynamics, as described by the PMS space, of the studied groups of responders and nonresponders before the treatment. Subsequently, we considered dynamic sensitivity analysis to determine the optimal perturbation strategies to rebalance the PMS distribution to the healthy state (as defined by the PMS space of responders after 1 day after treatment). Figureillustrates the dynamic sensitivity analysis, whereby the bifurcation parameter a is used to change the nodal dynamics in terms of its response to added noise, ranging from a more noise-driven regime (the more a is negative) to an oscillatory regime (with larger amplitude, the more a is positive). We focused on homological nodal perturbation of the large-scale brain model, meaning that bilateral regions were perturbed equally, resulting in 45 pairs of regions perturbed at gradually varying values of a. Figureshows the dynamic sensitivity analysis of driving a transition to the healthy state for models of both responders and nonresponders before treatment. Again, an average of the KL divergence between either the perturbed pre-treatment responder or non-responder models and the healthy PMS space was shown. In the noise-driven regime (a < 0), a deterioration of the fit was observed for both groups, while in the oscillatory regime (a > 0), an initial improvement across all 45 runs was depicted, before subsequent deterioration away from the optimal fit for both groups. Conversely, when replacing the target healthy state by the depressive state (i.e. by comparing with the average PMS in nonresponders after treatment), we found that the KL divergence was minimal without perturbation (i.e. keeping a = 0), showing a worse fit for both groups when brain areas became more oscillatory and no effect of the noisy perturbation (Fig.). This analysis demonstrates the optimal response for the amelioration of brain dynamics towards the optimal brain dynamics, as described by the PMS space of posttreatment responders, to be in the oscillation-driven regime (a > 0) for both responders and non-responders. To evaluate which regions permitted transition to a healthy state, we first defined the optimal perturbation strength as the minimum of the averaged KL divergence (across the 45 runs) of the responder group to the treatment. This stimulation intensity was found at a = 0.07. Then, we inspected the difference between the responders and nonresponders at that given value of a to assess what nodal perturbations were permitting the transition to the healthy state in responders but not in non-responders (Fig.). Figureshows the rank-ordered regional differences in KL divergence between perturbations of the responder and nonresponder models before treatment at a stimulation intensity
FigureEvaluation of dynamic sensitivity analysis. (A) Perturbation to induce a transition to a healthy state. Each homological pair of brain regions was perturbed by varying the bifurcation parameter a, which modulates the intrinsic oscillatory behaviour of the dynamical units. The more a is positive, the larger the amplitude of intrinsic oscillations, whereas for negative a, the units decay to a fixed point equilibrium and the local dynamics is dominated by noise. The performance of the perturbations is evaluated by computing the KL divergence between the simulated PMS and the empirical PMS from patients who recovered after treatment with psilocybin. Optimal intensity of a = 0.07 was achieved for the responder group (highlighted rectangles). (B) Perturbation to induce a transition to a depressive state. A transition to the depressive state showed worse or no effect at varying values of the bifurcation parameter a. This is expected since the models were adjusted to patients in the depressive state before treatment. of a = 0.07. We highlighted the regions with the largest KL divergence working in responders but not non-responders to promote a transition to the healthy state. These regions are the temporal superior pole, rolandic operculum, fusiform gyrus, supplementary motor area, parietal inferior gyrus, angular gyrus, supramarginal gyrus, frontal inferior gyrus (opercular), frontal middle gyrus (orbital) and parahippocampal gyrus. Figureshows the cortical rendering of these differences. In summary, these regions are implicated in transition away from 'depressed brain' pathology and towards the 'healthy brain' configurations of treatment responders.
Given the unique and known neuropharmacology of psychedelics whereby psilocybin-the active component of magic mushrooms-binds with high affinity to the serotonergic receptors (mainly the 5-HT 2a ) and in effect increases neural activity, we assessed whether the regions working in responders but not non-responders overlapped with the 5-HT density maps derived from PET imaging data previously obtained by an independent research group.Figureand B show the correlation between the 5-HT 2a and 5-HT 1a receptor density maps and regional differences in KL divergence between perturbations of the responder and non-responder models before treatment at a stimulation intensity of a = 0.07 (same as in Fig.; Spearman's ρ = 0.227, P = 0.032, and Spearman's ρ = 0.284, P = 0.007, respectively). Figureshows non-significant correlations to other 5-HT components-namely the 5-HT 2b (Spearman's ρ = 0.064, P = 0.055) and 5-HT 4 receptors (Spearman's ρ = 0.055, P = 0.607) and the 5-HT transporter (5-HTT; Spearman's ρ = -0.172, P = 0.106). This analysis shows that the ability to force transition towards the optimal spatio-temporal dynamics of a given region correlates with density of 5HT 2a and 5HT 1a neuroreceptors of that region.
In this work, we employed a large-scale brain modelling approach to evaluate potential brain change causes of response to psilocybin therapy for treatment-resistant depression. Using a novel combination of empirical data and in silico modelling, systematic perturbations to brain regions modelled in silico revealed a subset of regions implicated in transition away from 'depressed brain' pathology and towards the 'healthy brain' configurations of treatment responders. Notably, these regions matched those with the highest density of 5HT 2a and 5HT 1a neuroreceptors. This relationship is plausible given that psilocin (psilocybin's active metabolite) is known to have an appreciable-to-high affinity for the 5-HT 1a and 5-HT 2a receptors, respectively, where it acts as an agonist, in the case of the 5-HT 2a , potentially stimulating plasticity-related signalling cascades relevant to an antidepressant action.Another complementary perspective of the findings, which show the correlation between 5-HT 2a and 5-HT 1a serotonergic receptors and regions with difference in the intervention between the groups, is that of the overactivation of DMN in depression both in terms of its connectivity profile and temporal characterization.It has been shown that serotonergic raphe nuclei have direct anatomical projections to the regions of DMN and 5HT 2 receptors expressed in DMN regions,and it is thus possible that the psychedelic action enables modulation of the DMN via serotonergic pathways in responders but not non-responders. A summary of complex spatio-temporal dynamics, in terms of brain substates and their transitions, has drawn a lot of attention in systems neuroscience due to its utility to evaluate the impact of pharmacological and electromagnetic interventions for treating brain and behavioural disorders. Brain substates have been characterized in different ways, by minimal energyas attractor landscapesand, more heuristically, through sliding window analysis and unsupervised clustering.However, it has been challenging to find a model that is sufficiently simple and yet accurate to account for temporally and spatially complex and non-stationary data sets. Here, PMSs are built on a description of the data in terms of a probabilistic 'cloud' in substate space and as such can be extended to different neuroimaging modalities with higher temporal resolution, such as EEG and MEG, or potentially to more fine-grained spatial resolutions.Cutting-edge non-invasive brain stimulation techniques, such as transcranial magnetic stimulation and direct electrical stimulation, and new neuropsychopharmacological drugs for the treatment of psychiatric disorders have heralded a new era of localized and system-wide brain perturbations as medical interventions. For example, transcranial magnetic stimulation has been considered for the treatment of many psychiatric disorders, such as depression, schizophrenia and addiction,and classic psychedelic (drug) therapy, which, in part, targets a specific neuroreceptor (i.e. principally the 5-HT 1a receptor), is showing efficacy in the treatment of a broad range of conditions such as depressive, anxiety and addiction disorders.However, it seems highly likely that the mechanistic action of these interventions lies potentially well downstream of their initial action, and this action may not be straightforward.For example, how direct electrical stimulation-induced signal propagates within neuronal microcircuits remains unclear and often paradoxicaland motivates theoretical neuroscience studies and in silico perturbation protocols. In this perspective, our study attempted to understand how the psychedelic-induced perturbation impacts regional activity that alters the brain-wide dynamics, thus potentially bridging the initial (as done experimentally) and downstream actions (as done in silico) of such intervention. Beyond in silico perturbations, the exhaustive stimulation protocol can also be used as a dynamic sensitivity analysis tool from the complex systems perspective. Traditionally, statistical differences in measures summarizing spatiotemporal dynamics are obtained using signal detection theory. Such approaches can be enhanced by considering large-scale brain models and their structural differences between conditions, for example as described by the global coupling (G) parameter. Moreover, rather than describing and assessing expressions of spatio-temporal dynamics, an exhaustive protocol allows a shift of focus onto transitions to a target state, and this can be used to identify differences between groups, such as treatment responders versus nonresponders, as we have done here. Forcing transitions in large-scale brain networks has also been investigated through the prism of control network theory. In such scenarios, control strategies are deployed to navigate complex systems from a source (initial) state to a target (final) state.This approach has obtained a lot of attention due to its wide-ranging engineering applicability in technological, social and cyberphysical systems across various experimental scenarios.However, the conceptual understanding of controlling neuronal signals from source to target might be problematic as the brain operates in self-sustained and non-equilibrium state, and the notion of well-defined pathway between them might be ill posed.On the contrary, the approach considered in this work describes spatio-temporal dynamics in terms of PMSs and, through systematic perturbation, rebalances the spatiotemporal dynamics between two PMS spaces. Through this approach, the brain is rebalanced to its healthy working point, without specific instructions of what the relevant pathway might be. To obtain a PMS approximation of the brain substate of interest, several methodological choices are made, which inevitably introduce several caveats. First, a regional parcellation must be chosen, which might introduce artificial spatial boundaries especially when dealing with dynamics. Secondly, the choice of clustering algorithm defines the type of substates that can be obtained. Here, we use the unsupervised learning algorithm k-means clustering, which has been shown to adequately represent functionally meaningful brain substates.However, alternative algorithms could be used for this purpose (e.g. k-medoids). Related to the experimental data, the design is an uncontrolled open-label feasibility pilot study, and as such has no placebo group and suffers from small sample size. In this context, it is relevant to consider the current study as exploratory. Hence, future replication studies are warranted to ensure robustness of the findings. Moreover, the healthy state is defined here in terms of the 1-day post-treatment scan but the responders/ non-responders' assessment is done 5 weeks after. Lastly, the large-scale brain models constructed are based on group approximations of the functional brain information and SC group template. For clinical relevance, further research will be needed to create individual-based large-scale brain models that might allow for future in silico-assisted personalized psychiatry.
Create a free account to open full-text PDFs.
Lord, L. D., Expert, P., Atasoy, S. et al. · NeuroImage (2019)
Daws, R. E., Timmermann, C., Giribaldi, B. et al. · Nature Medicine (2022)
Carhart-Harris, R. L., Goodwin, G. M. · Neuropsychopharmacology (2017)
Atasoy, S., Leor, R., Kaelen, M. et al. · Scientific Reports (2017)
Singleton, S. P., Luppi, A. I., Carhart-Harris, R. L. et al. · Nature Communications (2022)
Carhart-Harris, R. L. · Current Opinion in Psychiatry (2019)
King, C., Nichols, D. E. · Nature Reviews Neuroscience (2013)
Deco, G., Cruzat, J., Cabral, J. et al. · Current Biology (2018)
Carhart-Harris, R. L., Roseman, L., Bolstridge, M. et al. · Scientific Reports (2017)
Papers in Blossom that reference this study
Brendstrup-Brix, K., Ozenne, B., Fisher, P. M. et al. · Journal of Psychopharmacology (2026)
Urban, M. M., Zillich, L., Rieser, N. M. et al. · Translational Psychiatry (2026)
Socoró-Garrigosa, M., Sanz Perl, Y., Kringelbach, M. L. et al. · Annals of the New York Academy of Sciences (2025)