In summary, a quantitative characterization of the spatio-temporal dynamics recorded with fMRI was obtained using leading eigenvector dynamics analysis (LEiDA), resulting in the definition of Probabilistic Metastable Substates (PMS), whose probability of occurrence was compared across conditions (i.e., within-subjects design – therefore, before versus after treatment). We then constructed two whole-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. 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 endpoint) one-day post treatment brains).

As described in the 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 (SI Fig. 2). 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. 2, A). Furthermore, we also computed the Global Brain Connectivity (GBC), metastability and Functional Connectivity Dynamics (FCD) measures (see SI Fig. 1). These results clearly indicated the necessity of considering both spatial and temporal dimensions to differentiate between conditions as GBC, 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 post-treatment 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 whole-brain dynamics across conditions.

To obtain whole-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 (51)) was described by a Stuart-Landau oscillator (see methods), and regions were coupled according to realistic structural connectivity 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 (see SI Fig. 3). 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. 2B. In SI Fig. 4, we report optimisation 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 minimise the difference. Figure 2C shows 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).

Subsequently, we considered a 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 one-day after treatment). Figure 3 illustrates 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 whole-brain model, meaning that bilateral regions were perturbed equally, resulting in 45 pairs of regions perturbed at gradually varying values of *a*. Figure 3, Left, shows the dynamic sensitivity analysis of driving a transition to the healthy state for models of both responders and non-responders before treatment. Again, an average of the KL divergence between either the perturbed pre-treatment responders 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 non-responders 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. 3, Right).

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 non-responders 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. 4).

Figure 4.A shows the rank ordered regional differences in KL divergence between perturbations of the responder and non-responder models before treatment at a stimulation intensity 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 the Parahippocampal gyrus. Figure 4, Right, shows the cortical rendering of these differences.

## Correlation with Serotonin Receptor Maps

Given the unique neuropharmacology of the psychedelic-induced state through serotonergic receptors, 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 (39). Figure 5A shows correlations between the *5-HT**2a* and *5-HT**1a* receptor density maps and the KL divergence differences for the two groups at optimal *a* = 0.07 (Spearman \(\rho\) = 0.227, p = 0.032 and Spearman \(\rho\) = 0.284, p = 0.007 respectively). Figure 5B, shows non-significant correlations to other 5-HT components – namely the *5-HT**2b* (Spearman \(\rho\) = 0.064, p = 0.055) and *5-HT**4* receptors (Spearman \(\rho\) = 0.055, p = 0.607) plus the 5-HT transporter (*5-HTT*) (Spearman \(\rho\) = − 0.172, p = 0.106).