Using individualised whole‑brain computational modelling fitted to fMRI and diffusion MRI, the authors simulated administration of LSD and psilocybin and applied in silico perturbations to compare states of consciousness and assess treatment effects in disorders of consciousness. Simulated psychedelics shifted patients' brain dynamics closer to criticality—especially in the minimally conscious state—with responses in UWS linked to structural connectivity and in MCS to baseline functional connectivity, providing a computational rationale for psychedelic therapies and personalised prediction.
Disorders of consciousness (DoC), including unresponsive wakefulness syndrome (UWS) and minimally conscious state (MCS), have limited treatment options and are characterized by low complexity of brain activity. Recent research suggests that psychedelic drugs, which enhance the complexity of brain activity, could offer promising therapies. Here, individualized whole‐brain computational models are developed for patients with DoC, optimized with empirical functional magnetic resonance imaging data and diffusion‐weighted imaging data, upon which the administration of lysergic acid diethylamide (LSD) and psilocybin is simulated. An in silico perturbation protocol is applied to assess brain dynamics, first distinguishing between different states of consciousness, including DoC, anesthesia, and the psychedelic state. Then, brain dynamics are assessed before and after a simulation of psychedelic drugs on patients with DoC. Findings indicated that the simulation of LSD and psilocybin shifted the brain activity of patients with DoC closer to criticality (the point at a phase transition between order and chaos), with a greater effect in patients in the MCS. In patients with UWS, the treatment response correlated with structural connectivity, while in patients in the MCS, it aligned with baseline functional connectivity. These results offer a computational foundation for using psychedelics in DoC treatment and highlight the potential future role of computational modeling in drug discovery and personalized medicine.
L. and colleagues situate their study within the view that conscious experience depends on spatiotemporally complex brain dynamics supported by structural connectivity (white matter, from diffusion-weighted imaging) and functional dynamics (BOLD fMRI). Healthy brain dynamics are argued to operate near a critical point — a transition between order and disorder — which affords high sensitivity to perturbation, long-range correlations and a balance of integration and differentiation thought to underpin rich conscious experience. Disorders of consciousness (DoC) are described as states of pathologically reduced complexity and stability (sub‑critical), whereas psychedelic states show increased complexity and sensitivity to perturbation. Because administering psychedelic drugs to DoC patients raises practical and ethical challenges, the authors propose using whole‑brain computational models — digital twins — to explore mechanistic effects and perform a "phase zero" virtual clinical trial simulating drug administration. The study sets out to (1) create individualized whole‑brain models for patients with DoC (based on each patient's DWI and fMRI), (2) validate a model‑based perturbational biomarker, the perturbative integration latency index (PILI), for discriminating states of consciousness, and (3) virtually simulate administration of LSD and psilocybin to each patient model to evaluate whether psychedelics shift DoC brain dynamics toward criticality. The authors also aim to identify baseline structural or functional biomarkers that predict the simulated treatment effect, thereby assessing personalised treatment potential with in silico methods.
L. and colleagues combined six previously published datasets to build group‑level and individual whole‑brain models. Key samples cited include a DoC dataset with 35 healthy controls (CNT), 20 unresponsive wakefulness syndrome (UWS) patients and 26 minimally conscious state (MCS) patients; anesthetic datasets (propofol N = 13, dexmedetomidine N = 11); and psychedelic datasets (LSD N = 12, psilocybin N = 15, plus an auxiliary psilocybin dataset). For DoC subjects, diffusion MRI (DWI) and resting‑state fMRI from the same individuals were used; models for healthy drug datasets used an average structural connectome derived from the DoC healthy controls because individual SC was not available for those datasets. Ethical approvals and consent procedures are reported for all source studies. Preprocessing: All fMRI data were reprocessed with a common pipeline (FSL/MELODIC). Steps included discarding initial volumes, motion correction, ICA denoising with manual component inspection, band‑pass filtering (0.01–0.08 Hz; modelling analyses used 0.04–0.07 Hz for signal extraction), registration to T1, and extraction of regional BOLD time series from the 90 AAL atlas regions (cortical and subcortical, cerebellum excluded). Pairwise Pearson correlations produced 90×90 functional connectivity (FC) matrices (Fisher z‑transformed for group averaging). Structural connectivity: Patient structural connectomes were built from single‑subject tractography (MRtrix3 pipeline). Preprocessing included denoising, eddy and susceptibility correction, multi‑tissue constrained spherical deconvolution, probabilistic tractography with 20 million streamlines per subject, SIFT2 weighting, and mapping to the AAL parcellation to produce weighted 90×90 SC matrices. Fractional anisotropy (FA) and graph‑theory measures (graph strength, global/local efficiency, characteristic path length, betweenness centrality) were computed. Computational model and fitting: Each brain region was modelled as a Stuart‑Landau (normal form of a supercritical Hopf) oscillator coupled according to SC and scaled by a global coupling parameter G. Local bifurcation parameters a j determined node dynamics (a < 0 noisy fixed point; a = 0 at bifurcation; a > 0 oscillatory). For group models where individual SC was unavailable, the average healthy SC from the DoC control group was used. Model fitting proceeded in two steps: (1) optimise global coupling G by minimising the Kolmogorov–Smirnov (KS) distance between empirical and simulated FC across G values (500 simulations per G), and (2) fit local a parameters constrained to a nine‑dimensional space defined by nine resting‑state networks (derived from ICA of the DoC healthy controls). A genetic algorithm was used to find the optimal nine network a values (range −0.2 to 0.2), and each region's a j was the linear combination of the networks it belongs to. Perturbation and PILI: The perturbation protocol switched a single region's bifurcation parameter to a synchronous regime (a = 0.6) for 100 s, repeating this for each of the 90 regions and over 100 trials per region. Integration at each time point was calculated from instantaneous phases (Hilbert transform) by binarising phase‑locking across thresholds and measuring the size of the largest connected component; integration is the integral across thresholds. The perturbative integration latency index (PILI) was computed as the area between the perturbed integration curve and the maximal basal integration until return to baseline, capturing both amplitude and latency of the return to baseline. Regional PILI values were averaged across trials; a global mean PILI was computed as the average across regions. Virtual pharmacology: To simulate LSD and psilocybin administration in patient models, the authors derived parameter shifts (changes in G and the nine network a parameters) from group‑level drug versus placebo models (healthy subjects). These extracted shifts were applied to each patient model to produce a simulated drug state. The method was validated with two cases of empirical fMRI from patients scanned before and during a light propofol dose: parameter shifts derived from group models of light propofol were applied to patient models and compared with empirical propofol‑state models. Statistics: At the group model level, Cohen's d was used to compare PILI distributions. For individual models, conventional frequentist tests were used: one‑way ANOVA and post‑hoc Tukey HSD for three‑group comparisons, Mann–Whitney U tests, Wilcoxon signed‑rank tests for paired comparisons, paired t‑tests for region‑wise changes (Bonferroni correction across 90 regions), and correlation analyses with Bonferroni correction where indicated. The authors note using each simulation as a datapoint for group‑level models would artificially inflate sample size, hence reliance on effect sizes there.
PILI distinguishes states of consciousness: Using group‑level Hopf models fitted to FC for multiple states, PILI discriminated reduced‑consciousness states from richer phenomenological states. DoC groups and anaesthesia showed reduced PILI compared with their controls, while psychedelic conditions showed increased PILI relative to placebo. Reported effect sizes (Cohen's d) included UWS versus CNT = −0.52 and MCS versus CNT = −0.48; propofol versus wake = −0.54 and dexmedetomidine versus wake = −0.32. Psychedelic comparisons yielded smaller positive effects (LSD Cohen's d = 0.11; psilocybin Cohen's d = 0.28). Network‑wise PILI changes were tabulated in the Supporting Information. Individualised models and group differences: Individualised patient models (CNT N = 35, UWS N = 20, MCS N = 26) showed a significant main effect of diagnostic group on mean regional PILI (one‑way ANOVA F(2,78) = 8.84). Post‑hoc Tukey tests indicated higher PILI in CNT versus UWS (mean difference = 6.80, 95% CI [2.94, 10.66], p < 0.001) and higher PILI in MCS versus UWS (mean difference = 4.11, 95% CI [0.09, 8.12], p = 0.044). No significant difference was found between CNT and MCS (mean difference = 2.69, 95% CI [−0.92, 6.30], p = 0.182). All p‑values reported had multiple‑comparison correction applied. Virtual clinical trial—simulated psychedelics: Applying parameter shifts derived from healthy drug datasets to individual DoC models, simulation of LSD and psilocybin increased PILI across patients with DoC, with larger effects in MCS than UWS. For LSD, increases were greater in MCS (z = 4.13, p < 0.001) than UWS (z = 2.8, p < 0.001); for psilocybin, UWS z = 3.71, p < 0.001 and MCS z = 3.92, p < 0.001. Region‑wise analyses showed that, in MCS patients, the largest network‑level PILI increases after simulated LSD and psilocybin occurred in the default mode network (DMN), attention and frontoparietal networks. In UWS patients, the sensorimotor network exhibited the highest PILI increase for both drugs. Region‑wise PILI did not show significant correlations with regional structural connectivity (SC) or functional connectivity (FC) after Bonferroni correction (reported uncorrected correlations: e.g. LSD MCS r = 0.21, p = 0.04; others non‑significant). Validation with propofol: As a validation of the virtual pharmacology method, parameter shifts for a light propofol dose were applied to two patient models with empirical pre/post propofol scans. The mean differences in residual per‑region PILI between simulation‑derived and empirical‑derived models were reported as ~0.0002 in one patient and ~1 in another; one patient with low baseline PILI showed a small simulated effect, while a patient with higher baseline PILI showed a larger decrease after propofol. The authors interpret this as the method capturing differential individual effects, with further discussion in Supporting Information. Baseline predictors of simulated treatment response: The simulated treatment effect was defined as change in PILI from baseline to simulated drug. Correlational analyses revealed that, in UWS patients, mean SC correlated with simulated treatment effect for both LSD (r = 0.66, p = 0.001) and psilocybin (r = 0.66, p < 0.002), surviving Bonferroni correction. In contrast, in MCS patients, mean FC correlated with simulated treatment effect (LSD r = 0.65, p < 0.001; psilocybin r = 0.57, p = 0.002), and network‑specific FC (notably DMN FC) also correlated with treatment effect for both drugs; structural measures did not significantly correlate with outcome in MCS. Representative patient examples illustrated that high baseline mean FC predicted strong simulated response in MCS, whereas in UWS high mean SC distinguished those with larger simulated responses. Other graph metrics showed group differences (e.g. characteristic path length lower in UWS than MCS; global efficiency and graph strength differed between controls and UWS), but mean SC was the dominant predictor in UWS.
L. and colleagues interpret their findings as demonstrating that PILI, a model‑based perturbational biomarker, differentiates states of consciousness: diminished states (DoC, anaesthesia) had lower PILI consistent with dynamical rigidity and sub‑critical regimes, while psychedelic states showed higher PILI consistent with increased sensitivity and proximity to criticality. Using personalised whole‑brain models (digital twins), simulated administration of LSD and psilocybin increased the brain's perturbational response in patients with DoC, shifting dynamics closer to criticality; this effect was larger in MCS than UWS. The authors emphasise that increased PILI indicates greater spatially extensive and sustained integration after perturbation — properties associated with complexity — but caution that a change in PILI need not directly equate to improved behavioural consciousness or beneficial phenomenology. The authors note dissociations in baseline predictors: in UWS patients structural integrity (mean SC) correlated with simulated response, implying that when structural scaffolding is severely compromised it limits capacity for psychedelic‑induced complexity; in MCS patients, functional connectivity — particularly within the DMN, attention and frontoparietal networks — predicted response, suggesting preserved functional organisation is necessary to support dynamical changes. They present this as a hierarchical relationship where white‑matter architecture provides the scaffold on which functional dynamics and complexity emerge. The authors discuss practical and conceptual implications: virtual pharmacology and digital twins could enable mechanistic investigations, drug discovery and personalised treatment testing without exposing vulnerable patients to risk. They acknowledge important limitations: PILI shows substantial within‑subject variability and currently lacks single‑subject diagnostic precision comparable to measures like the perturbational complexity index (PCI); cross‑dataset differences in acquisition complicate comparisons; the Hopf model is phenomenological and not mechanistic at a biophysical level; and the virtual pharmacology relies on a superposition assumption that parameter shifts observed in healthy subjects generalise linearly to patients with DoC. The authors also recognise reliance on parameter estimates from single datasets as a source of potential bias. They stress that empirical studies in DoC patients would be required to validate simulated predictions and refine parameter sets. Overall, the authors present their virtual clinical trial as a computational proof‑of‑principle that psychedelics can increase modelled brain complexity in DoC and as a step towards computational personalised medicine and drug discovery, while explicitly noting many open questions about phenomenological and behavioural consequences and the need for empirical validation.
We built whole-brain computational models, optimized at the individual patient level, to explore using psychedelic drugs as a treatment for DoC. We assessed brain dynamics using an in silico perturbation protocol, which measures the return to baseline dynamics after a modeled perturbation. Firstly, using group-level models, we showed how the response to perturbation can discriminate between different states of consciousness (i.e., DoC, anesthesia from propofol, anesthesia from dexmedetomidine, and the psychedelic state under LSD and psilocybin). Then, we simulated the administration of LSD and psilocybin separately upon individualized models of patients with DoC and assessed the modeled dynamics via response to perturbation at baseline and at the simulated psychedelic state. Lastly, we characterized the empirical baseline predictive biomarkers of simulated treatment effect.
We first sought to establish an in silico perturbational biomarker that can distinguish between states of consciousness, based on the level of consciousness, and phenomenological richness of experience. This could later serve as a proxy to assess the simulated treatment response in our individual patient models. To achieve this, we created Hopf whole-brain computational models optimized on the FC matrix obtained from the BOLD fMRI data corresponding to each state of consciousness, averaged at the group level (Figure). We used previously published datafrom several states of consciousness to create 13 different models: three for the DoC dataset (see Table, Supporting Information, for the demographics of the patients with DoC) (healthy controls (CNT) N = 35, UWS N = 20, and MCS N = 26), four models based on healthy subjects under anesthesia (propofol N = 13, dexmedetomidine N = 11, as well as a separate model of wakefulness for each dataset), four models based on healthy subjects in the psychedelic state (LSD N = 12, psilocybin N = 15), and the two related placebo conditions. We also performed a supplementary analysis using a model built on another psilocybin dataset, presented in Supporting Information. Briefly, each model possesses two sets of parameters: the global coupling parameter (G) and the local bifurcation parameters (a). The global coupling parameter linearly scales the coupling between brain regions in the SC matrix. The local bifurcation parameters denote the dynamics of each node: negative values produce noisy activity, while positive values lead to stable oscillatory dynamics. After retrieving the optimal global coupling value, the local parameter optimization was restricted to a nine-parameter space representing nine resting state networks (attention, auditory, DMN, frontoparietal, limbic, sensorimotor, precuneus, thalamus, visual) defined by an independent components analysis on the fMRI data from the healthy controls of the DoC dataset (see Model Fitting to Empirical Data in the experimental section for more details on optimization). Optimal global parameters and local parameters for each model used in this study can be found in Tablesand(Supporting Information). We then investigated the dynamical properties of the different states of consciousness following an in silico perturbation. Each model was exposed to a perturbation protocol, which consisted of shifting the dynamics of one of the 90 brain regions to a more stable state for 100 s, repeated for each brain region. The integration, a measure calculated from the phase coherence synchronization between brain regions (see Intergration in experimental section for further details) over time was then determined immediately after the perturbation, and in a baseline state without perturbation. We then calculated PILI for each model after perturbation by summing the area beneath the perturbational integration curve until it returned to the baseline state (Figure). We determined the PILI for each node by selectively perturbing one node while maintaining the dynamics of the other nodes unperturbed, in their basal state, across 100 trials (see Model Perturbation Protocol and Perturbative Integrative Latency Index in the experimental section for more details). These trials were then averaged to produce a single PILI value for each node. Finally, a singular, mean PILI value representing the average across all brain regions was computed. To statistically compare the PILI values between states of consciousness, we took each simulation as a separate datapoint, thus having 90 brain regions and 100 simulations for each region. We calculated the Cohen's d effect size to compare the distributions between each condition and its associated comparison. In states of reduced consciousness, DoC patients and healthy subjects under anesthesia had decreased PILI values compared to the comparison conditions (CNT and wakefulness, for DoC and anesthesia groups, respectively) (Figure). Specifically, UWS (Cohen's d = -0.52) and MCS (Cohen's d = -0.48) patients had lower PILI values compared with CNT. Similarly, decreases in PILI values were also observed for propofol (Cohen's d = -0.54) and dexmedetomidine anesthesia (Cohen's d = -0.32). In contrast, we found that PILI values in the psychedelic conditions were higher compared to the placebo condition, for both LSD (Cohen's d = 0.11) and psilocybin (Cohen's d = 0.28). Network-wise changes in PILI are listed in Table(Supporting Information).
In dynamical-systems theory, the critical point denotes the state that sits exactly between order and chaos. At this phasetransition point, several interesting properties emerge. These include scale-free distributions of activity (spanning all spatial and temporal scales), long-range correlations, and high sensitivity to perturbations, which result in extensive and sustained responses.The healthy brain is thought to operate in a zone near the point of criticality or marginally sub-critical.Moving away from criticality in either direction is defined by changes in two properties of the system. Stability increases when the system shifts toward the sub-critical, highly ordered regime. Stability is inversely related to chaos and describes how rapidly a system returns to baseline following perturbations. Highly stable (subcritical) systems produce brief, localized perturbational responses. In contrast, systems nearer to criticality show more sustained and widespread integration; thus, a greater response to perturbation. Shifting beyond the critical point into the chaotic (super-critical) domain renders the system hypersensitive, but at the cost of losing the differentiated activity. Complexity is inherently related to criticality. Both terms are multidimensional concepts that cannot be captured with a single metric. Complexity, as defined by Tononi et al.,refers to a balance between global integration and local differentiation of neural activity. Measures like Lempel-Ziv complexity quantify some facets of this concept, but are not all-encompassing. For example, a signal of "chaotic" noise would have a high Lempel-Ziv complexity, despite having a suboptimal balance depending on the value of the bifurcation parameter, can be at a stable fixed-point (a < 0), a stable limit cycle (a > 0), and a bifurcation between both regimes (a = 0). The SC provided by white matter tractography from DWI modulates the coupling between brain regions, and the global coupling parameter scales the coupling between brain regions in the SC. The local bifurcation parameters are optimized by minimizing the Euclidean distance of the simulated FC towards the FC of the empirical fMRI data, restricted to a parameter space representing nine resting state networks obtained from an independent component analysis (one a parameter per network, with the a parameters of every region in the atlas being the linear combination of the a parameters of the networks the region belongs to). Brain dynamics are assessed by simulating a perturbation and observing the resultant return to a dynamical baseline. b) PILI protocol. We applied the Hilbert transform to the simulated BOLD signals to obtain the instantaneous phases. We constructed and binarized a phase locking matrix at each time point and calculated the number of regions in the largest connected component over thresholds. Integration was defined as the integral over all thresholds. The perturbation protocol consisted of modifying the bifurcation parameter of one brain region to the stable regime (a = 0.6) for 100 s. Integration was then calculated over 300 s in the basal unperturbed state and immediately after the model perturbation. PILI was computed as the integral between the curves of integration values over time for the perturbed dynamics (orange) compared to the maximum of the basal state dynamics (blue). of differentiation and integration. Because a system at criticality simultaneously preserves local differentiation and global integration, it is also where complexity, in the sense of "integrated differentiation," naturally peaks. The Perturbative Integration Latency Index (PILI) captures key dimensions of criticality and complexity. PILI measures how extensively and for how long perturbation propagates through the network: smaller, transient responses reflect lower complexity (and sub-critical states), while widespread and sustained perturbational responses indicate higher complexity and proximity to criticality.
We then constructed individualized whole-brain models for each subject in our DoC dataset (20 UWS, 26 MCS, 35 CNTs). The parameters of each model were optimized towards the individual BOLD fMRI data. To define the coupling between nodes, we used the individual patient SC, generated by white matter tractography from the DWI data. From this point onwards, all the analyses were performed using these digital twins of each subject. We performed the perturbation protocol on each patient individually and calculated PILI values for each of the 90 brain regions. In all our analyses on individual models, we decided to consider a single PILI value in each brain region as the mean across all 100 simulations. The results from a one-way ANOVA between the PILI values across MCS patients, UWS patients, and CNT showed a significant main effect of group (F(2, 78) = 8.84). Post-hoc comparisons using Tukey's honestly significant difference (HSD) test indicated that PILI values were higher in the CNT group compared to the UWS group (mean difference = 6.80, 95% CI [2.94, 10.66], p < 0.001). Additionally, PILI values were higher in the MCS group compared to the UWS (mean difference = 4.11, 95% CI [0.09, 8.12], p = 0.044). There was no significant difference between the CNT and MCS groups (mean difference = 2.69, 95% CI [-0.92, 6.30], p = 0.182) (Figure). All reported p-values were corrected for multiple comparisons.
The next steps were to simulate the administration of LSD and psilocybin on the individual patient models. This was achieved through developing a virtual pharmacology method (see Virtual Pharmacology Approach in the Experimental Section for more details). Briefly, this involved creating whole-brain models optimized at the group level using the BOLD fMRI data from healthy individuals being administered a drug and at the respective placebo or comparison conditions. For each drug condition individually, we extracted the changes in the global coupling and local bifurcation parameters that manifested upon comparing the placebo/comparison model to the drug model. We then applied the shift in parameters representing the virtual administration of the drug separately to each patient's model (Figure). The parameter changes representing the simulation of LSD and psilocybin can be seen in Table. The resulting stimulatory druginduced changes in brain dynamics were then characterized by calculating PILI at the baseline state (parameters in their optimized values for that specific patient) and after the simulation of the drug (after shifting the parameters that represent the drug effects). While direct data from DoC patients under psychedelics are unavailable, we sought to provide validation for our virtual pharmacology framework using data from propofol anesthesia. We had access to two unique cases in which BOLD fMRI data were acquired before and during the administration of a light dose of propofol (average effect site concentration 1.80 μg mL -1 ). Our propofol dataset in healthy individuals had two levels of propofol dose: a light dose of propofol (1.75 μg mL -1 average blood plasma concentration), resulting in a Ramsay score of 3, and a moderate dose (3.20 μg mL -1 average blood plasma concentration), resulting in an unresponsive state (Ramsay score of 5/6). Therefore, we set out to validate our model of a virtual administration of LSD and psilocybin by simulating a virtual administration of a light dose of propofol and comparing the resulting simulated brain dynamics to the models built on the empirical data of the patient under a light dose of propofol. Using our virtual pharmacology approach, we optimized models based on empirical BOLD fMRI for normal wakefulness and a light dose of propofol and extracted the changes in parameters between these conditions before applying them to the individual patient models. We then assessed the brain dynamics through PILI at the baseline, after empirical propofol, and after the simulation of propofol (Figure). The mean differences in residual per-region PILI between the simulationderived model and empirical-derived model were 0.0002 in one patient and around 1 in another. One patient with a low baseline PILI had a relatively small effect from the simulation of propofol, whereas another patient with a relatively high baseline PILI had a somewhat larger decrease in PILI as a result of propofol. The simulation of propofol using our virtual pharmacology method seems to capture these differential effects. Please consult Supporting Information and Figure(Supporting Information) for a longer discussion about the validity of our virtual pharmacology approach. Subsequently, we simulated the administration of LSD and psilocybin to each individual patient model. Simulating LSD significantly increased the PILI values of patients with DoC, with a greater effect on patients in the MCS (z = 4.13, p < 0.001) compared to patients with UWS (z = 2.8, p < 0.001). This was also the case after simulating psilocybin, as the increases in PILI values observed in patients with UWS (z = 3.71, p < 0.001) were lower than those observed in patients in the MCS (z = 3.92, p < 0.001). Next, we examined how the global response to perturbation differed across initial stimulation site, calculating the PILI values in each region, averaged across subjects at the baseline, and after simulating LSD and psilocybin. To better contextualize the regional changes in PILI, we calculated the differences in PILI in the nine resting state networks, which define our local parameter space (Table). In the MCS group, after the simulation of both LSD and psilocybin, the networks with the largest differences in PILI values were the DMN, the attention, and the frontoparietal network. In the UWS group, the sensorimotor network had the highest increase in PILI values for both the LSD and psilocybin conditions. For both LSD and psilocybin, region-wise PILI values did not correlate with regional SC after Bonferroni correction for the four comparisons (LSD: MCS r = 0.21, p = 0.04, UWS r = 0.05 p = 0.56; psilocybin: MCS r = 0.13, p = 0.23; UWS r = 0.02, p = 0.84), nor FC (LSD: MCS r = 0.17, p = 0.11, UWS r = 0.12, p = 0.84; psilocybin: MCS p = 0.26, r = 0.02, UWS r = 0.05, p = 0.61).
We then sought to identify which structural and functional capacities of the brains of patients with DoC at baseline were associated with a high simulated treatment effect. In this way, we could ascertain the most promising candidates for potential psychedelic treatment. We defined the simulated treatment effect to be the difference in PILI between the baseline state and after the simulation of LSD and psilocybin. Within each of the MCS patients, UWS patients, and the associated control group, we quantified the strength of FC across the nine resting state networks that define our local bifurcation parameter space. Additionally, we computed several graph theory metrics from the patients' SC matrices, including global efficiency, local efficiency, centrality, and graph strength. We also determined the mean FA values from DWI data. See TablesandHere, we present four representative example patients and their structural and functional characteristics at baseline to elucidate insights into which characteristics could predict treatment efficacy (Figure). These include one patient from each diagnostic group who had the largest increases and decreases in PILI, averaged between psilocybin and LSD. MCS S41 showed the largest increase in PILI values for both LSD and psilocybin. At baseline, MCS S41 had one of the highest mean FC strengths of the group (z = 4.46, p < 0.001). However, their mean SC strength was not significantly different from the group mean (z = 0.4, p = 0.74). UWS S07 had one of the highest simulated treatment effects of the group. This patient had a mean FC which was not significantly different from the group mean (z = 0.4291, p = 0.70). However, the mean SC was significantly higher than the group mean (z = 3.16, p = 0.002). Examining those subjects with the worst simulated treatment outcome, UWS S14 had a relatively high FC (z = 2.65, p = 0.008), yet one of the worst SC of the group (z = -4.20, p < 0.001). Whilst MCS S40 had a SC not significantly different from the mean of the group (z = 0.79, p = 0.43), yet a low mean FC (z = -2.98, p = 0.002). To further investigate these relationships and identify baseline biomarkers of predicted treatment efficacy, we performed correlation analyses between each structural and functional capacity and the simulated treatment effect. In the UWS group, only the mean SC (LSD: r = 0.66, p = 0.001; psilocybin: r = 0.66, p < 0.002) correlated with the simulated treatment effect and remained significant after Bonferroni correction (Figure). In the MCS group, there were no significant correlations between structural capacities and the simulated treatment effect for either drug (Table, Supporting Information). However, the mean FC strongly correlated with the simulated treatment effect for both LSD (r = 0.65, p < 0.001) and psilocybin (r = 0.57, p = 0.002) in the MCS group. Additionally, in both the LSD and psilocybin conditions, FC in the DMN significantly correlated with the simulated treatment effect, and for LSD, FC in the attention, frontoparietal network, and sensorimotor network also correlated with the simulated treatment effect (Table, Supporting Information).
We used individualized whole-brain computational models in patients with DoC based on empirical DWI and fMRI data. Using each digital twin, we conducted a virtual clinical trial to explore a novel treatment paradigm: using psychedelic drugs, specifically LSD and psilocybin, to increase the complexity of brain activity in patients with DoC. To this end, we first set out to validate PILI, a modeling-based perturbational biomarker based on the integration of a modeled perturbation in distinguishing between states of consciousness, namely, DoC, anesthesia, and the psychedelic state. Second, we simulated the administration of LSD and psilocybin on each patient and assessed the resulting brain dynamics as a proxy for simulated treatment response. Third, we characterized the baseline predictive biomarkers associated with increases in PILI, ultimately identifying predictive biomarkers of potential treatment efficacy. Assessing the brain's response to a perturbation allows us to understand how it integrates external information and provides insights into to the criticality of brain dynamics and their implied complexity. Our results showed that PILI could distinguish between states of consciousness according to the phenomenological richness of experience. States of diminished consciousness, such as DoC (UWS and MCS) and anesthesia (propofol and dexmedetomidine), showed lower PILI values compared to the control subjects and wakefulness conditions, respectively. States with a higher richness of experience, ergo the psychedelic state (LSD and psilocybin), possessed higher PILI values compared to their placebo conditions. This supports research linking the complexity of brain activity and criticality to consciousness,ideas which are at the core of the integrated information theory of consciousness.States in which consciousness is diminished, such as in DoC or anesthesia, possess more stable brain dynamics, with lower complexity.This dynamic rigidity, manifesting through a faster return to the baseline following a perturbation, exemplifies the shift towards a sub-critical regime. Support from TMS-EEG studies on anesthesia and DoC shows lower evoked complexity in these states compared to normal wakefulness.During the psychedelic state, the increase in PILI epitomizes the brain's shift even closer to criticality. This echoes several previous lines of investigation showing that the psychedelic state is associated with an increase in brain complexity, an increased repertoire of available states, and more critical dynamics.Both propofol and dexmedetomidine were administered at unresponsive doses. However, dexmedetomidine often induces dreaming, which occurs less often with propofol.This likely contributes to dexmedetomidine's intermediate PILI value and underscores the importance of considering phenomenological differences between unresponsive states.When comparing PILI values from diagnostic groups at the single-subject level, there was a difference in mean PILI values between the UWS and MCS groups, but not between the MCS and CNT groups. This mirrors findings from studies on the PCI, where, below a numerical cut-off, a patient is considered unconscious, and above this level, there is less discrimination between conscious states like MCS and CNT.It also aligns with clinical assessments of consciousness at the bedside, where UWS patients do not show any signs of consciousness, yet MCS patients exhibit signs of residual consciousness. The established threshold allows PCI to function as a diagnostic tool at the single-subject level. While PILI currently lacks this discriminative capacity at the single-subject level, future developments in computational an independent set of two patients receiving a light dose of propofol. Baseline and empirical propofol show the PILI values before and after receiving propofol. Simulated propofol shows the PILI values after simulating propofol. c) Individual patient models before and after simulation of LSD and psilocybin. Blue circles represent each patient at baseline, orange diamonds represent LSD, and green diamonds represent psilocybin. d) Left: Plotted t-statistic values of region-wise t-tests displaying brain regions with significant increases, in PILI from simulating LSD and psilocybin in patients with UWS (top) and MCS patients (bottom), Bonferroni corrected for multiple comparisons for the 90 brain regions. Right: Average PILI values for each of the 90 brain regions, averaged across patients at the baseline state (X-axis) and after the simulation of LSD (orange) and psilocybin (green) (Y-axis). ) Examples of patients and their functional and structural capacities at baseline. One example is chosen from each diagnostic group that underwent the highest and lowest average simulated treatment effect between psilocybin and LSD conditions. Left: Absolute regional changes in PILI values for each subject. Warmer colors (red and yellow) indicate higher PILI changes, while cooler colors (green and blue) indicate lower PILI changes. Right: Radar plots comparing the functional and structural capacities of each patient (filled grey shape) to the mean of the associated diagnostic group (blue line). The left radar plot illustrates the average FC across all regions and within the networks that define the local parameter space: frontoparietal, default mode, auditory, attention, visual, thalamus, sensorimotor, precuneus, and limbic networks. The right radar plot shows average SC metrics across all regions, and the average global efficiency, graph strength, local efficiency, betweenness centrality, fractional anisotropy, and characteristic path length. modeling should aim to create a diagnostic tool of comparable precision. We also developed a method to administer virtual pharmacological treatments based on extracting the changes in parameters observed due to the drug effects and applying them to models that represent other states of consciousness, here, patients with DoC. In principle, this approach could be generalized to extract parameter changes that represent any pharmacological or neuromodulatory intervention that acutely results in changes in brain connectivity, such as anesthetics, stimulation techniques, or even simulating other conditions, such as listening to music. These changes could then be applied to any model at the group, or single-subject level to simulate such conditions. Crucially, the success of this method depends on the parameter extraction originating from a dataset that estimates each state within the same subjects. This ensures the most accurate estimation of the treatment effects by imposing the specificity of functional and structural changes. To illustrate the importance of this, we present (Figure, Supporting Information) and discuss (Supporting Information) analyses using a separate fMRI dataset from healthy subjects under psilocybin, acquired with a between-subjects design. Here, the PILI was slightly lower in the psilocybin condition compared to the placebo condition, and the simulated treatment response in patients with UWS was not associated with the SC. The PILI increases in response to the administration of LSD and psilocybin to patients with DoC. This represents a shift in brain dynamics towards criticality, with greater effects observed in the MCS compared to UWS groups. The simulated psychedelic-induced changes in brain dynamics in MCS patients show a similar absolute shift towards criticality as those seen between UWS and CNT, and between anesthesia and wakefulness. Whilst at face value, this may imply that the absolute size of the change in complexity could be enough to support a conscious state, the phenomenological implications and behavioral consequences of these changes in PILI remain an open question. It is conceivable that the administration of a psychedelic to patients with DoC may invoke a supraphysiological state, which does not improve consciousness. Predicting such phenomenological effects is outside of the remit of this study, and early clinical work, such as a recently published case study, can begin to approach these questions.In patients in the MCS, after both LSD and psilocybin, some of the largest network-wise increases in PILI were observed in the DMN, attention, and frontoparietal. Notably, the DMN and frontoparietal networks also showed the largest absolute increases in PILI in group-level models of LSD and psilocybin in healthy subjects. Each of these higher-order networks is associated with supporting facets of consciousness and cognition and shows diminished connectivity in patients with DoC. Particularly, patients with DoC show diminished DMN connectivity, with patients in the MCS having stronger DMN connectivity compared to patients with UWS.Therefore, alterations in the dynamics within these networks suggest that they may impart some changes in experience. Interestingly, in patients in the MCS, the DMN, frontoparietal, and attention networks showed the highest correlation between baseline connectivity and the simulated treatment effect. This implies that the strong connectivity in these networks at baseline is necessary to support the functional and dynamical changes entailed by the administration of a psychedelic drug, thus any potential treatment response. In the UWS group, the only network with significant changes in PILI for both psilocybin and LSD was the sensorimotor network. This is likely due to the large cortical structural damage present in other networks, restricting the propagation of activity into such areas. The lower-order role of the sensorimotor network in initiating and coordinating movements rather than being associated with consciousness per se, decreases the likelihood of treatment outcomes positively augmenting consciousness and cognition in patients in the UWS. Whilst generally consistent, psilocybin and LSD had some differential effects, with psilocybin showing more variability and a greater effect in patients with UWS. Although both drugs produce similar phenomenological effects,they have distinct pharmacological profiles: LSD targets D1-3 receptors, while psilocin, the active metabolite of psilocybin, inhibits the serotonin transporter.Psilocybin and LSD have been shown to have unique patterns of functional reorganization dependent on distinct neurotransmitter mechanisms.Also, recent modeling work simulating the stimulation of different receptors in patients with DoC has demonstrated the differential effects these receptors have on brain dynamics. Compared to 5-HT2A receptor stimulation, activation of D1 and D2 receptors had a minimal impact on the modeled trajectory of brain dynamics in patients with DoC towards that of controls.These differences in pharmacology and connectivity likely contribute to the variations in their simulated effects on brain dynamics and highlight their potential differing utility in treatments. In patients with UWS, the SC was correlated with the simulated treatment response, whereas in patients in the MCS, it was the FC that was predictive. Neither the network FC connectivity, nor average FC correlated with simulated treatment response in patients with UWS, while in patients in the MCS, none of the structural characteristics correlated with treatment response. Despite there being no significant difference in mean SC between MCS and UWS patients, characteristic path length was significantly lower in patients with UWS compared to patients in the MCS. Additionally, global efficiency and average graph strength were significantly different between controls and patients with UWS, yet not between controls and patients in the MCS. Graph strength refers to the sum of connection weights linked to each node. Global efficiency reflects the average inverse shortest path length between all node pairs, while characteristic path length is the average shortest path length across the network. Each of these metrics reflects different aspects of network efficiency. Reductions in connectivity within hub-regions would result in decreases across these measures. Thus, decreases in efficiency measures without a significant drop in mean SC may indicate regional disconnection, particularly affecting hub regions, while maintaining connectivity or even overcompensating in other regions in patients with UWS. Since efficient organization is already severely impaired inpatients with UWS, the primary limiting factor for enabling a dynamic response to psychedelics becomes whether any structural scaffolding remains at all, rather than how efficiently it is organized. In contrast, in MCS patients, sufficient structural efficiency is preserved, such that the total strength of SC becomes less important. Once a minimal level of efficiency of communication is maintained, FC rather than SC becomes the primary determinant of capacity for psychedelic-induced changes in complexity. Our results point to a hierarchical relationship between structure, function, and complexity. Efficient white-matter architecture provides the first-layer scaffold that allows spatially differentiated regions to exchange information, upon which richer functional dynamics can emerge. If that scaffold falls below a minimum level of efficiency, complexity cannot be supported. Consistent with this view, mean SC predicted treatment outcome in UWS patients more strongly than any specific graph-theoretic metric. In practical terms, topological details of network efficiency contribute relatively little once damage is extensive, since the sheer number of intact white-matter fibers cap a patient's capacity for complex brain activity. This implies a bleak prognosis for patients with UWS whose hub regions have suffered major structural dam-age. However, there is a growing interest in another use case of psychedelics in those with severe brain injury that seeks to harness the neuroplastic effects to augment recovery immediately in the acute phase of brain injury.Future work utilizing models that account for plasticity could explore this. Taken together, our correlational analyses suggest that a patient in the MCS with high average FC and specifically within the DMN and attention networks, has the greatest likelihood to have the largest changes in brain dynamics following a psychedelic drug. In silico experiments with personalized generative wholebrain computational models hold substantial potential for mechanistic investigations, drug discovery, and ultimately personalized medicine. Creating a unique digital twin of each patient offers a powerful methodological approach for simulating pharmacological treatments or neuromodulation techniques, without exposing patients to actual risks, in effect creating a virtual clinical trial.The considerable heterogeneity in the etiology, lesion site, and symptoms of patients with DoC makes this approach particularly appealing, simulating interventions and subsequently tailoring them for each patient with respect to drug choice, dosage, stimulation type, etc. Notable progress has already been made in epilepsy research, where personalized models have provided a valuable means to test stimulations and simulate post-surgical outcomes.As computational resources grow and methodologies refine, digital twins are poised to become highly valuable in clinical neuroscience, enabling personalized treatment strategies, and providing mechanistic insights. This study has several limitations. The inherent range in PILI values within each subject at the single-subject level underscores the limited applicability of PILI for between-subject comparisons, emphasizing its more appropriate use for characterizing states within the same subjects. Similarly, differences in acquisition protocols and scanning parameters make comparing PILI values between datasets challenging, as evidenced by the variability in PILI values across different control group models (see Figureand Supporting Information). The Hopf model is a phenomenological model that does not aim to replicate underlying biology. However, recent studies have shown that many complex bottom-up models at the regional level merely replicate a Hopf bifurcation.Therefore, while the Hopf model is limited in its mechanistic explanatory power compared to biophysical models, it still similarly captures regional and global brain dynamics. This work relies on the superposition assumption, being that a linear subtraction of the effects of psychedelics on controls well predicts the action of the drugs in patients with DoC. Rather than a uniform shift in the local/global parameters models, it is likely that patients with DoC have a different response to psychedelics compared to healthy individuals. Despite this significant assumption, it provides a useful approximation in the absence of empirical data on patients under a psychedelic. Future empirical data may be used to refine these parameter sets for patients with DoC. Additionally, our method is based on parameter estimations from a single dataset, which could have bias. Here, we ran a virtual clinical trial by simulating the administration of LSD and psilocybin in individualized whole-brain models of patients with DoC. We showed that psychedelics increase the brain's response to perturbation and globally bring the brain dynamics of patients with DoC closer to criticality. We also reveal that patients in the MCS are more likely to increase in complex- ity compared to patients with UWS and that this is dependent upon the strength of the baseline FC. This work builds on the idea of using psychedelics for DoC and contributes toward computational personalized medicine and drug discovery.
Dataset Descriptions: This study included previously published fMRI data from six datasets, three from the University and University Hospital of Liège (DoC, propofol, dexmedetomidine), two from Imperial College London (LSD, psilocybin), and one from Maastricht University (psilocybin). For the DoC dataset, DWI acquired in the same patients and healthy controls was also used. An auxiliary dataset of DWI and fMRI data from two patients was also used: One UWS and one MCS patient, before and after they received a light dose of propofol. See Dataset Descriptions in the Supporting Information for further details of the acquisition and demographic information. Figurebriefly describes each step of analysis in the study and where each dataset was used. The DoC, propofol, and dexmedetomidine datasets received approval from the Ethics Committee of the Faculty of Medicine at the University of Liège (DoC: 2009/241, propofol: 2007/191, dexmedetomidine: 2012/135). Written informed consent was obtained from the healthy participants and, in the case of patients with DoC, from their legal representatives. The LSD and psilocybin studies from Imperial College London were approved by the UK National Health Service Research Ethics Committees in West-London (LSD: 13/LO/0229) and Bristol (psilocybin: 08/H0101/144), respectively. All participants provided informed consent. The psilocybin study from Maastricht Univer-sity obtained ethical approval from Maastricht University's Medical Ethics Committee (NL60352.068.17 / METC173006). Informed consent was obtained from all participants. Functional Magnetic Resonance Imaging and Functional Connectivity Estimation: Since the datasets were acquired at different locations, we endeavored to mitigate scanner and acquisition differences by repreprocessing all data with the same pipeline and following uniform denoising procedures. The preprocessing of fMRI images was conducted utilizing FSL (FMRIB Software Library v6.0, Analysis Group, FMRIB, Oxford, UK). Structural T1-weighted images were processed using the fsl_anat function, which involved brain extraction, bias field correction, and normalization. Functional MRI data underwent preprocessing with MELODIC (Multivariate Exploratory Linear Optimized Decomposition into Independent Components). The Maastricht psilocybin fMRI data uniquely underwent realignment and susceptibility distortion correction with FSL_topup to correct for magnetic field distortion caused by the high-strength 7T field. Main preprocessing steps for all datasets included discarding the initial five volumes, motion correction with MCFLIRT, brain extraction using BET (Brain Extraction Tool), spatial smoothing with a 5 mm FWHM Gaussian kernel, rigid body registration, and applying a high-pass filter with a cutoff of 100.0 s. Additionally, single-session ICA with automatic dimensionality estimation was performed in order to identify noise-driven components. Each component was visually inspected using the FSLeyes package in Melodic mode to classify the single-subject independent components as either signal or noise/lesion-driven artifacts, based on the spatial map, time series, and temporal power spectrum. Noise components were then regressed out using fsl_regfilt to obtain a cleaned version of the functional data. Furthermore, BOLD data were detrended, demeaned, and band-pass filtered between 0.01 and 0.08 Hz. To obtain the BOLD time series for the 90 cortical and subcortical brain regions defined by the AAL atlas (excluding the cerebellum), FSL tools were used in each individual's native EPI space. This parcellation was found to be particularly suitable for studying whole-brain spatiotemporal dynamics.Given that the computational models were personalized at the individual level, utilizing 90 regions provided a balanced approach between spatial resolution and computational feasibility. The cleaned functional data were registered to the T1-weighted structural image using FLIRT. The T1-weighted image was then normalized to the standard MNI space using FLIRT (12 DOF) and FNIRT. The resulting transformations were concatenated, inverted, and applied to warp the resting-state atlas from MNI space to the cleaned functional data. Nearest-neighbor interpolation ensured the preservation of labels. The BOLD time series for each of the 90 brain regions was extracted for each subject in their native space using fslmaths to create a binary mask of each brain region and fslmeants to obtain the time series of each mask. Each of the 90 regional average BOLD signals was then filtered in the 0.04-0.07 Hz range, which has previously been used due to its functional relevance and resistance to noise.Pairwise Pearson correlation coefficients were then computed between all 90 brain regions for each subject. For the group-level matrices, the Fisher transform was applied to the r-values to derive z-values for the 90×90 functional connectivity (FC) matrices before being averaged and subsequently back-transformed. Anatomical Connectivity: The patient structural connectomes were created using advanced, single-subject analyses of the diffusion MRI (dMRI) dataset of patients with DoC, including 35 age-matched controls. MRI scans were performed using a Siemens 3T Trio scanner (Siemens Inc., Munich, Germany) with a 64-channel head coil. Using tools from MR-trix3, raw dMRI data were corrected,denoised, and corrected for Gibbs ringing artifacts,as well as distortions induced by motion, eddy currents, and EPI/susceptibility artefacts.The response functions from the pre-processed diffusion-weighted images were estimated using the Dhollander algorithm and used to calculate fiber orientation distributions (FODs) via constrained spherical deconvolution (CSD).Here, MRtrix3Tissue (v5.2.8;) was utilized, a variant of MRtrix3,which enabled three-tissue constrained spherical deconvolution, generating separate FODs for white matter (WM), grey matter (GM), and cerebrospinal fluid (CSF) SS3T-CSD,from single-shell (+b = 0) diffusion MRI data. Whole-brain tractography was then performed,generating 20 million streamlines per subjectusing dynamic seeding. The spherically informed filtering of tractograms (SIFT2) algorithm was applied to align the fiber density of the reconstructed streamlines with the underlying white matter structures.Compared to tractograms constructed solely by the number of streamlines, SIFT2 adjusted the weight of individual streamlines to ensure alignment with the underlying image data.SIFT2 is highly reproducibleand enhances the biological interpretability of the estimated white matter tracts.Probabilistic tractography was then performed using the iFOD2 (improved 2nd order integration over fiber orientation distributions) algorithm, which enhanced anatomical plausibility. During tracking, the direction for each step was obtained by sampling from the FOD at the current position, with the probability of a particular direction being proportional to the FOD amplitude. Additional tractography settings included a default step size of 0.8 mm, a maximum angle between successive steps of 45°, a maximum length of 250 mm, a minimum length of 5 mm, a cutoff fractional anisotropy (FA) value of 0.08, and the use of b-vectors and b-values from the diffusion-weighted gradient scheme in FSL format. Concurrently, the b0 image in native space was registered to the T1 anatomical image using FSL FLIRT.The T1 structural image was then normalized to standard space using FLIRT and FNIRT. The resulting transformations were reversed and applied to the AAL atlas in standard MNI space using a nearest neighbor algorithm to convert the atlas to each subject's native space. The number of streamlines connecting each pair of regions was then estimated to produce the final SC matrix for each patient. To summarize, for each participant, a 90×90 symmetric weighted network of SC was constructed, which represented the density of white matter tracts between anatomical regions of the brain. The FA values were generated by firstly estimating the diffusion kurtosis tensor map from the pre-processed DWI image. A map of tensor-derived FA for each patient was then averaged across voxels to produce one FA value per subject. Graph theory metrics were calculated using The Brain Connectivity Toolboximplemented in MATLAB (version R2017a, MathWorks, Natick, MA). Graph strength, global efficiency, characteristic path length, average local efficiency, and betweenness centrality were calculated to capture different aspects of network structure. These measures provided insights into node importance, communication efficiency, and overall network topology. Computational Whole-Brain Model: A whole-brain model was used to simulate the brain activity of each of the 90 cortical and subcortical brain regions in the AAL atlas. In particular, the dynamics of each region were described by means of coupled Stuart-Landau oscillators in the normal form of supercritical Hopf bifurcations.This model simulated spontaneous signals akin to the average BOLD signal in a given brain region. The coupling between oscillators representing each region was driven by a given SC. The model at node j was described by the following set of coupled differential equations: where z is a complex variable z j = x j + iy j , w j is the natural node frequency (estimated as the average peak frequency of each region of the BOLD recordings of each group or subject, for the group and single-subject models, respectively), 𝜂 j is additive Gaussian noise, 𝛽 is the scaling factor for the Gaussian noise (set to 0.04), C is the SC matrix, where C jk is the coupling between nodes j and k, G is the global coupling that acts as a scaling factor for the SC matrix for all nodes, and a j is the bifurcation parameter controlling the behaviour of the oscillator. For a ≈ 0 the system is at a supercritical Hopf bifurcation, and the Gaussian noise induces complex dynamics as the system oscillates between both sides of the bifurcation. For a < 0, the system decays to a stable fixed point at z j = 0 and noise-induce oscillations appear due to the Gaussian noise. For a > 0, the system dynamics set into a stable limit cycle with a frequency of f j = 𝜔 j 2𝜋 . The SC matrix was scaled to a maximum value of 0.2 to ensure oscillatory dynamics bifurcations.The fMRI BOLD signal simulation was obtained as the real part of z (x j ). The stochastic differential equations were solved by means of the Euler-Maruyama method. For all the group models involving healthy subjects, including those for psychedelics, anesthesia, and the associated wakefulness and placebo comparisons, the SC matrix representing the average of the healthy controls in the DoC dataset was used. This approach was necessary as the authors did not have access to the individual SC data for each subject in each dataset. In previous modeling research, a template SC obtained from a commonly used dataset was usually employed, which represented the white matter connectivity of healthy subjects. This was assumed to well represent the SC of healthy controls in general.Here, the same approach was adopted by using the SC estimate from the healthy control group of the DoC dataset. This ensured that the strength of the connectivity values, and thus the optimal parameters for the models employing a healthy structure, were in the same range as those of patients with DoC. Crucially, models of patients with DoC included the empirical structural data for each patient. The diverse etiology of patients with DoC often resulted in significant structural differences between patients. By using personalized SC, the inherent bias of assuming a healthy white matter structure in DoC patients was avoided, which enabled building more robust models incorporating more empirical information. Individualized singlesubject whole-brain models were created using the individual structural and functional data. Model Fitting to Empirical Data: This work involved the creation of seven models at the group level and 46 models at the single-subject level. For the group-level analyses, the FC averaged over each subject in each condition was used to obtain the parameters for each group model. For the single-subject analyses, the individual FC of each patient was used to fit their parameters. Two parameters were fitted to the models in two steps. First, the global coupling G was fitted to the FC matrix by exploring Adv. Sci. 2025, e11780 e11780 (13 of 17) a parameter range in steps of 0.01 starting from 0 until the fitting curve reached an optimum value. Three hundred seconds of BOLD activity were simulated each time. The same filtering applied to the empirical data was applied to the simulated time series, and their FC matrices were obtained. Afterwards, G was optimized based on the Kolmogorov-Smirnov distance (KS distance) between the empirical and simulated FC. For each G value, 500 simulations were performed, and the curve of the average KS distance was used to determine the optimal G value at the local minimum. Using the KS distance, the distribution of the FC values of the simulations was brought to a similar level to that of the empirical data. After fixing the G values for each group/subject as a starting point, in a second step, the a values were also fitted in order to fine-tune their dynamics by adjusting the bifurcation parameters of the regions. For each region j, the bifurcation parameter a j was obtained as the linear combination of nine resting state networks including that region. These resting state networks were obtained by means of a group independent component analysis (ICA) of the fMRI dataset from the healthy controls in the DoC dataset. This was implemented through the MELODIC ICA functionality of FSL, with automatic dimensionality estimation. Each component was manually inspected and identified as either noise or a resting state network. This resulted in the identification of nine resting state networks. These were then named according to their similarity with canonical resting state networks (default mode [DMN], auditory, attention, visual, thalamus, sensorimotor, precuneus, frontoparietal, limbic). Where there was not a corresponding canonical resting state network, the network was named based on the brain areas involved, e.g., thalamus. These IC maps were then thresholded and binarized to produce masks. Each network mask was then compared to the AAL atlas to see which regions overlapped with each component mask at the voxel level. This produced network components in the AAL space. The IC maps and corresponding AAL expression can be seen in Figure(Supporting Information). Thus, the parameter space was reduced to a dimension of 9 (one for each network). The Euclidean distance was used as the goodness of fit parameter to optimize the fit for each region, and a genetic algorithm to find the optimal values of the a parameters of each network as in refs. [43, 94]. Nine initial values (one for each network) were simulated, the default parameters of the MATLAB function ga were used, the algorithm was run ten times, and the average a values of the networks over the course of the ten runs were used as the final optimal value. Ten was chosen as a balance between computational power requirements and finding the true optimal value. For each network, the parameter range of its associated a values to explore by the genetic algorithm was set to -0.2 to 0.2.Integration: Integration refers to the capacity of the brain to assimilate communication between discrete processing units. The integration measure used in this study evaluates the synchronization of brain activity measured by the BOLD signal and was previously defined in refs. [31, 45]. The Hilbert transform was applied to filtered BOLD signals to extract the instantaneous phase 𝜑 j (t) for each region j. The Hilbert transform derives the analytic representation of a real-valued signal, in this case, the BOLD time series. The synchronisation between pairs of brain regions was characterized as the difference between their instantaneous phases. At each time point, the phase difference P jk (t) between two regions j and k was calculated as Here, P jk (t) = 1 when the two regions are in phase, 𝜑 j (t) = 𝜑 k (t) (full synchronized). At any time t, the phase interaction matrix P(t) represents the instantaneous phase synchrony of each pair of brain regions. The phase interaction matrix was then binarized over 100 evenly spaced thresholds between 0 and 1. Then, the size of the largest connected component across all brain regions was obtained at each threshold, thus creating a function of connected components as a function of the threshold. The integration at time t was then denoted as the integral of this curve (Figure). The integration measure reflected the synchronization of brain region dynamics at different thresholds. This reflected communication at all possible scales afforded by the atlas, including long-range communication, a hallmark of criticality. A larger connected component signified extensive synchronization across brain regions, with enhanced capacity for flexible, long-range communication. Conversely, smaller components reflected localized, less integrated dynamics. Model Perturbation Protocol and Perturbative Integrative Latency Index: The perturbation protocol has been previously described in refs. [31, 45]. Briefly, the local bifurcation parameters of the Hopf model were used to simulate a perturbation through switching each a from its original base value to a more synchronous regime (a = 0.6) for 100 s. This value was based on previous studies employing the PILIand has been shown to well differentiate conscious states. The integration was computed over 400 s of simulated BOLD time series in the basal state and immediately after the perturbation offset. This protocol was performed at every brain region in the 90 AAL parcellation and repeated 100 times. After perturbation, the bifurcation parameter was reset to its original value obtained from the model fitting procedure. Plotting the integration at each time point in the perturbed state created a curve showing how the integration changes over time in response to a perturbation; the basal integration line represented the integration in the absence of any perturbation, with the alpha parameters kept at their optimized values. PILI in a particular region was calculated as the area under the perturbational integration curve from perturbation offset to the interception point with the basal state line of integration. This was repeated for each brain region for 100 trials before being averaged to obtain a regional PILI value. The singular mean PILI represented the average PILI over all brain regions. Since the PILI was based on the integral of the integration curve, it characterized both the strength and latency of the return to its basal state after a model perturbation. PILI quantified how extensively perturbations propagate through the network. A small response to perturbation reflected lower complexity, whereas a large response to perturbation indicated widespread integration of multiple regions, consistent with higher complexity and a system closer to criticality.Slowness of recovery to perturbation is a key property of critical systems, where it is sometimes referred to in the literature as "critical slowing".A biological mechanism underlying increases in the local bifurcation parameter was suggested to be akin to increasing the local excitability.With this extension, such kinds of perturbations resembled the effects of those seen with transient increases in cortical excitability induced experimentally by methods like anodal transcranial direct current stimulationVirtual Pharmacology Approach: To simulate the administration of psilocybin and LSD to the individualized patient models, a virtual pharmacology method was developed. This was based upon group-level models optimized to the BOLD fMRI data from healthy individuals being administered a drug (LSD, psilocybin) and at the respective placebo conditions. To simulate the administration of the target drug, the global coupling parameters and local bifurcation parameters representing each condition were then extracted and compared within each condition to attain a shift in the parameters that represented that particular condition. This change in parameters thus represented a set of instructions on how to shift the model dynamics from one condition to another. This shift was then applied in parameters representing the virtual administration of the drug (LSD and psilocybin) separately to each patient's model. As a validation step, the shift in parameters was also extracted from group models optimized to the BOLD fMRI data from healthy individuals before and after a light dose of propofol, and a virtual administration of propofol in two patients was simulated. The effect of modeled brain dynamics after the virtual administration of propofol would then be compared with modeled brain dynamics from these two unique cases, where BOLD fMRI data before and after the administration of a light dose of propofol were present. Statistical Analyses: When comparing the mean PILI values between whole-brain models built at the group level, Cohen's d effect size was used to compare the distributions of PILI values produced by each model. Since each group level model can in theory, produce an infinite number of simulations, thus artificially inflating the sample size, any p-values produced by frequentist statistical tests were less applicable in these cases. For all analyses using individual models, a standard frequentist approach was employed since each model essentially can act as a unique datapoint with its own distribution, akin to a subject. To compare the PILI values between diagnostic groups, from the personalized models of DoC patients, Mann-Whitney U tests were performed. To compare the mean PILI values obtained from simulating the psychedelic drug to the baseline PILI values, Wilcoxon signed-rank tests were performed. Paired t-tests were used to evaluate the region-wise changes in PILI as a result of simulating the psychedelic with Bonferroni correction for the 90 regions. Wilcoxon signed rank tests were used to compare the network-wise PILI values between the baseline state and simulated drug state with Bonferroni correction for each of the nine resting state networks.
Create a free account to open full-text PDFs.