https://orcid.org/0000-0003-3084-0177
Figures
Abstract
A topic of growing interest in computational neuroscience is the discovery of fundamental principles underlying global dynamics and the self-organization of the brain. In particular, the notion that the brain operates near criticality has gained considerable support, and recent work has shown that the dynamics of different brain states may be modeled by pairwise maximum entropy Ising models at various distances from a phase transition, i.e., from criticality. Here we aim to characterize two brain states (psychedelics-induced and placebo) as captured by functional magnetic resonance imaging (fMRI), with features derived from the Ising spin model formalism (system temperature, critical point, susceptibility) and from algorithmic complexity. We hypothesized, along the lines of the entropic brain hypothesis, that psychedelics drive brain dynamics into a more disordered state at a higher Ising temperature and increased complexity. We analyze resting state blood-oxygen-level-dependent (BOLD) fMRI data collected in an earlier study from fifteen subjects in a control condition (placebo) and during ingestion of lysergic acid diethylamide (LSD). Working with the automated anatomical labeling (AAL) brain parcellation, we first create “archetype” Ising models representative of the entire dataset (global) and of the data in each condition. Remarkably, we find that such archetypes exhibit a strong correlation with an average structural connectome template obtained from dMRI (r = 0.6). We compare the archetypes from the two conditions and find that the Ising connectivity in the LSD condition is lower than in the placebo one, especially in homotopic links (interhemispheric connectivity), reflecting a significant decrease of homotopic functional connectivity in the LSD condition. The global archetype is then personalized for each individual and condition by adjusting the system temperature. The resulting temperatures are all near but above the critical point of the model in the paramagnetic (disordered) phase. The individualized Ising temperatures are higher in the LSD condition than in the placebo condition (p = 9 × 10−5). Next, we estimate the Lempel-Ziv-Welch (LZW) complexity of the binarized BOLD data and the synthetic data generated with the individualized model using the Metropolis algorithm for each participant and condition. The LZW complexity computed from experimental data reveals a weak statistical relationship with condition (p = 0.04 one-tailed Wilcoxon test) and none with Ising temperature (r(13) = 0.13, p = 0.65), presumably because of the limited length of the BOLD time series. Similarly, we explore complexity using the block decomposition method (BDM), a more advanced method for estimating algorithmic complexity. The BDM complexity of the experimental data displays a significant correlation with Ising temperature (r(13) = 0.56, p = 0.03) and a weak but significant correlation with condition (p = 0.04, one-tailed Wilcoxon test). This study suggests that the effects of LSD increase the complexity of brain dynamics by loosening interhemispheric connectivity—especially homotopic links. In agreement with earlier work using the Ising formalism with BOLD data, we find the brain state in the placebo condition is already above the critical point, with LSD resulting in a shift further away from criticality into a more disordered state.
Author summary
In this study, we aim to characterize two brain states (psychedelics-induced and placebo), as captured in functional magnetic resonance imaging (fMRI) data, with features derived from the Ising model formalism (system temperature, critical point, susceptibility) and from algorithmic complexity. Under the hypothesis that psychedelics drive the brain into a more disordered state, we study criticality features of brain dynamics under LSD in a within-subject study using the Ising model formalism and algorithmic complexity using Lempel-Ziv and the Block Decomposition methods. Personalized Ising models are created by first using BOLD fMRI data from all the subjects and conditions to create a single Ising “archetype” model—which we can interpret as the average model of the data at unit temperature—and then by adjusting the model temperature for each subject and condition. We find that the effects of LSD translate into increased BOLD signal complexity and Ising temperature, in agreement with earlier findings and predictions from existing theories of the effects of psychedelics, such as the relaxed beliefs under psychedelics (REBUS), the anarchic brain hypothesis, and the algorithmic information theory of consciousness (KT). However, in contrast with some of the previously cited theories, we find that the system in the placebo condition is already in the paramagnetic phase—above the critical point—with ingestion of LSD resulting in a shift away from Ising criticality into a more disordered state. Finally, we highlight the fact that the structural connectome can be recovered to a good degree by fitting an Ising model and that the reduction of homotopic links (direct or indirect) appears to play an important role in the slide to disorder under psychedelics.
Citation: Ruffini G, Damiani G, Lozano-Soldevilla D, Deco N, Rosas FE, Kiani NA, et al. (2023) LSD-induced increase of Ising temperature and algorithmic complexity of brain dynamics. PLoS Comput Biol 19(2): e1010811. https://doi.org/10.1371/journal.pcbi.1010811
Editor: Marcus Kaiser, University of Nottingham, UNITED KINGDOM
Received: September 1, 2022; Accepted: December 11, 2022; Published: February 3, 2023
Copyright: © 2023 Ruffini et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The multimodal neuroimaging data from the experiment is available at https://openneuro.org/datasets/ds003059/versions/1.0.0 (doi:10.18112/openneuro.ds003059.v1.0.0). Python code for Ising modeling of BOLD data can be found at GitHub, https://github.com/giulioruffini/LSD-paper-code-Ruffini-et-al.-2022, for LZW analysis at https://github.com/giulioruffini/StarLZW, and for fMRI analysis at https://github.com/decolab/cb-neuromod. BDM calculation programs in different languages are available at the Algorithm Dynamics Lab website, https://www.algorithmicdynamics.net/software.html.
Funding: GR, GiD, DL-S, ND have received funding from the European Research Council (ERC Synergy Galvani) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 855109). GR, GD, GiD, DL-S and ND have received funding from European Union’s Horizon 2020 research and innovation programme under grant agreement No 101017716 (Neurotwin). FER is supported by the Ad Astra Chandaria foundation. NAK has received funding from the European Union under the Horizon 2020 programme (MultipleMS grant agreement 733161 and DECISION Grant agreement 847949). AP-A was supported by the Human Brain Project SGA3 (945539) and the Ramón y Cajal fellowship (RYC2020-029117-I). RC-H is supported by the Ralph Metzner Distinguished Professorship endowment. GD has received funding from the Spanish national research project (ref. AEI-PID2019-105772GB-I00/AEI/10.13039/501100011033) funded by the Spanish Ministry of Science, Innovation and Universities (MCIU), State Research Agency (AEI). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: I have read the journal’s policy and the authors of this manuscript have the following competing interests: authors with a Neuroelectrics affiliation (GR, GiD, DL-S and ND) work for a company creating brain stimulation solutions. GR is a co-founder of Neuroelectrics. The other authors have declared that no competing interests exist.
Introduction
The complex features displayed by the emergent neural dynamics of the human brain are in many ways reminiscent of behaviors of paradigmatic systems in statistical physics [1–9]. For example, just as condensed materials can transit between ordered (e.g., solid) or disordered (e.g., gas) phases, neural networks can undergo phase transitions between highly regular and highly unsynchronized activity. Of particular interest are the dynamics of systems poised close to phase transitions, where so-called critical points display a balance of order and disorder that has been conjectured to be crucial for living organisms [10]. At critical points, non-trivial collective patterns spanning all scales are observed, giving rise to a rich repertoire of short- and long-range correlations. The fluctuations in critical systems are highly structured, obeying fundamental physical principles rooted in the system’s symmetries—which are shared by many systems despite their underlying constituency, establishing ‘universality classes’ of large-scale coordinated activity [11].
Theoretical work indicates that computation and information dynamics in complex systems display special features at critical points [12–16]). Building on these insights, criticality has been suggested as a useful principle to account for the brain’s inherent complexity that is required to guide behavior in rich environments by processing various sources of information [2, 17–23]. Furthermore, it has been found that systems in nature often tune themselves into a critical state—a principle known as ‘self-organized criticality’ [24–26]. It has been proposed that neural systems may self-organize, evolutionary and developmentally, through homeostatic plasticity rules that result in the system arriving near critical states [27–29].
Building on this pioneering work, this paper aims to leverage tools from the statistical physics literature to study brain dynamics that support different states of consciousness. Here we focus on the altered state of consciousness induced by LSD, a psychedelic substance that is recently harnessing attention for its therapeutic potential for addressing mental disorders such as anxiety, depression, or substance addiction [30–34]. While our understanding of how LSD works at the neuronal and pharmacological levels is relatively well-developed [29, 35, 36], its effects at the whole-brain level are less understood. Yet, this is paramount to explain better and predict the outcomes of psychedelic interventions—such as psychedelic-assisted therapy [30, 37, 38]. It has been argued that the observed expansion of the repertoire of functional patterns elicited by hallucinogenic substances can be understood as being driven by an enhancement of brain entropy [29]. In turn, this may take place while the brain is moving closer to criticality, with this being related to the relaxation of high-level cognitive priors [39, 40], which in turn may lead to a favorable context for conducting psychotherapy [39, 41]. Studies on functional neuroimaging regarding LSD effects have shown initial evidence of the mechanistic alterations on brain dynamics at the network level, with the majority of the findings suggesting a relative weakening of usual functional configurations giving place to an increase of brain entropy, global functional integration, and more flexible brain dynamics [29, 30, 42–47].
Given the theoretical and empirical evidence suggesting that LSD induces a shift towards a broader repertoire of neural patterns enabled by more flexible dynamics, the present study aims to characterize the underlying neural mechanisms responsible for these effects via mathematical and computational modeling tools from statistical physics and the study of phase transitions. Specifically, we leverage maximum entropy models—the most parsimonious models one can build following statistical physics principles [48, 49]—to describe the coordinated activity of various brain areas as captured by BOLD signals. By building individualized computational models for LSD and placebo conditions, we aim to analyze differences in terms of the model’s temperature, a parameter that describes the system’s randomness. In particular, a low temperature amplifies distinctions between energy levels and hence a higher dissimilarity between the likelihood of occupancy of different states; conversely, a high system temperature ‘flattens’ the energy landscape and allows access to more states. Moreover, the system’s temperature situates the state in relation to the critical boundaries defining phase transitions. Finally, given the intuitive relationship between signal diversity, algorithmic complexity, and criticality, we explore complexity metrics characterizing the two conditions in the data with the hypothesis that complexity correlates with system temperature and hence experimental condition [40, 50, 51].
Materials and methods
Ethics statement
The study was approved by the National Research Ethics Service committee London-West London and was conducted in accordance with the revised declaration of Helsinki (2000), the International Committee on Harmonization Good Clinical Practice guidelines, and National Health Service Research Governance Framework. Imperial College London sponsored the research, which was conducted under a Home Office license for research with Schedule 1 drugs. Written informed consent was given, and data was de-personalized (anonymized) prior to analysis.
Experimental protocol and BOLD fMRI data
Our dataset consists of BOLD time series obtained under LSD effects and in a placebo condition from fifteen participants (N = 15). Details of the data acquisition process can be found in detail in [30, 52], and we only provide a summary here. A group of 15 participants (4 women; mean age, 30.5±8.0) were recruited by personal spoken communication and provided with written consent to participate in the experiment after a briefing and screening study for physical and mental health. This included an electrocardiogram (ECG), urine and blood tests, plus a psychiatric test and drug history. Participants under 21 years, with psychiatric diseases, family history of psychotic disorders, experience with psychedelics, pregnancy, problems with alcohol (>40 units per week), or any other medically significant condition that could affect the study were excluded. Participants attended twice the laboratory (once per condition), with a separation of 2 weeks. The order of these was alternated across participants without giving the information about the experimental condition used in each session. With the help of a medical doctor, a cannula was inserted into a vein in the antecubital fossa after the tests. The dose consisted of 75 mg of LSD via a 10 ml solution proportioned over two minutes with an infusion of saline afterward. Similarly, a placebo (10 ml saline) was injected intravenously over two minutes. After this, participants were suggested to lie relaxed with closed eyes inside the MRI scanner. Participants provided subjective reports of the drug effects after a 5 to 15 minute period from the administration of the drug. The peak effects were reported to occur 60–90 minutes after the dose. Subsequently, effects were generally constant for four hours. The scanning was done approximately 70 minutes after the dose and lasted a full hour. BOLD scanning was composed of three resting state scans of seven minutes. The middle one incorporated listening to two excerpts of music from two songs by ambient artist Robert Rich. For our analysis here, we exclude the scan with music. Thus, the dataset includes four scans, each of 7 minutes (two under LSD, two in the placebo condition) for each of the fifteen subjects (for a total of sixty). The data from the two scans in each condition are combined in the following analysis.
Parcellated data series
The fMRI BOLD data was preprocessed using FSL tools (FMRIB’s Software Library, www.fmrib.ox.ac.uk/fsl) with standard parameters and not discarding any Independent Component Analysis (ICA) components. This included corrections for head motion during the scans, which was within the normal, acceptable range for fMRI sessions. Also, FSL was used to generate parcellated (region-averaged) BOLD signals for each participant and condition, i.e., the signal was averaged over all voxels within each region defined in the automated anatomical labeling (AAL) atlas. Only the 90 AAL atlas cortical and subcortical non-cerebellar brain regions of interest (ROIs or parcels) were considered [53] (see Table A in S1 Text). For each subject and condition (LSD or placebo), we obtained a BOLD time series matrix of dimensions 90 x 217 (AAL regions x time points, where the sampling cadence was 2 s) from each session. The preprocessing pipeline of the BOLD data and the generation of the BOLD time series is described in more detail in [52]. Fig J in S1 Text provides plots of the FC for each condition, and Fig L group provides the FC per subject.
Data binarization
To build an Ising model from the data or carry out complexity analysis, we transformed each parcel’s BOLD fMRI data series into a binary format. This was done by using a threshold: each value of the region-averaged BOLD time series was assigned a value of + 1 if its value was greater than the threshold and −1 otherwise. The threshold was set to the median of the time series in the parcel of interest. Hence the median of the thresholded time series of each voxel was zero, and the entropy of each time series was maximal. The binarization was carried out independently for each condition, session, subject, and parcel. The binarized data from the two scans in each session were concatenated.
Data analysis of the fMRI series
Lempel-Ziv complexity.
Lempel-Ziv (LZ) refers to a class of adaptive dictionary compression algorithms that parse a string into words and use increasingly long reappearing words to construct a dictionary. Compression is achieved by pointing to the identifiers of words in the dictionary [54–56]. LZ algorithms are universally optimal—their asymptotic compression rate approaches the entropy rate of the source for any stationary ergodic source [56], i.e., . where ℓ is the length of the compressed string, n is the length of the string and
is the entropy rate.
LZ has been used across many studies to characterize brain state—from neurodegeneration [57–60] to anesthesia [61], disorders of consciousness [62, 63] to fetal consciousness [64]. The entropic brain hypothesis [43] proposes that within a critical zone, the entropy of spontaneous brain activity is associated with various dimensions of subjective experience and that psychedelics acutely increase both [39]. More recently, it has been proposed that algorithmic versions of entropy, which can be estimated by LZ and other methods such as BDM (discussed below), provide a more direct link with the subjective structure of experience [40, 50].
After applying the algorithm to a string of length n, we obtain a set of words in a dictionary of length c(n) that can be thought of as a first complexity metric [65]. The description length of the sequence encoded by LZ can be approximated by the number of words seen times the number of bits needed to identify a word, and, in general, . We note that LZ provides an upper bound on algorithmic complexity, but, given its reduced programming repertoire (LZ is the Kolmogorov complexity computed with a limited set of programs that only allow copy and insertion in strings [66]), it fails, in general, to effectively compress random-looking data generated by simple, but highly recursive programs, e.g., an image of the Mandelbrot set (deep programs [40, 67]). As a simple example of these limitations, consider a string concatenated with its bit-flipped, “time-reversed”, or dilated version. Such simple algorithmic manipulations will not be detected and exploited by LZ. Despite these limitations, LZ can be useful for studying the complexity of data in an entropic sense. To use LZW, we first transform numeric data into a list of symbols. Here we binarize the data to reduce the length of strings needed for the algorithm to stabilize. A reasonable strategy in the study of algorithmic information content in data is to preserve as much information as possible in the resulting transformed string. In this sense, using methods that maximize the entropy of the resulting series is recommended. Here, as discussed above, we use the median of the data as the threshold.
Two associated metrics derived from LZ compression are commonly used: c(n) and ℓ. Of the two, the latter is more closely related to Kolmogorov complexity or description length. A natural way to normalize this metric is to divide the description length by the original string length n, with units of bits per character. Here, we analyze synthetic lattice data generated by the Ising model by computing ρ0 for each lattice realization. We also analyze ρ0 by concatenating them for a global measure (with real data). The former allows for a plot of the mean and standard deviation of the measure with synthetic data generated using the Metropolis algorithm as a function of temperature.
To calculate a global measure of LZ complexity using the experimental data, we first concatenated the data from all subjects, conditions, and sessions, then, flattened it along the spatial dimension first and then time. Secondly, we compressed the flat binarized data with the LZ algorithm and saved the output “archetype dictionary”, which was used as the initial dictionary to compress the data from all the subjects in the two conditions, LSD and placebo (comparison by condition), and also per subject (comparison by condition and subject). The archetype dictionary was built by concatenating first the LSD data for all subjects and sessions and then the placebo data because we were forced to make a choice when flattening the array. However, we checked whether concatenating the placebo data before the LSD changed the results, and the difference was negligible. For details and code used to compute LZ complexity—we use the Lempel-Ziv-Welch (LZW) variant—see [55, 68].
Block decomposition method.
The Block Decomposition Method (BDM) [69] combines the action of two universally used complexity measures at different scales by dividing data into smaller pieces for which the halting problem involved can be partially circumvented in exchange for a huge calculation based upon the concept of algorithmic probability and extending the power of the so-called Coding Theorem Method (CTM). The calculation, however, has been precomputed and hence can be reused in future applications by exchanging time for memory in the population of a precomputed look-up table for small pieces. BDM is used to account for both statistical regularities and algorithmic ones and is sensitive enough to capture small changes in complexity while invariant to different descriptions of the same object. It, therefore, complements the use of lossless compression algorithms to calculate an upper bound of algorithmic complexity. We carry out the same analysis as for ρ0 described above, i.e., both on synthetic and experimental data.
Ising spin model
Abstract frameworks from statistical physics can shed light on understanding emerging phenomena in large networks, such as phase transitions in systems where nodes—neurons or cortical columns here—interchange information under the assumptions of the maximum entropy principle [48, 49, 70]. The description of systems with many degrees of freedom can be summarized by coarse-graining variables (describing macrostates), which introduces statistics into modeling. At the so-called critical points, observable quantities such as diverging correlation length or susceptibility to external perturbations reveal singularities in the limit as the number of degrees of freedom goes to infinity. At these transitions, from order to disorder, there is a loss of sense of scale, with fractal properties in energy and information flow.
In this context, elements such as neurons, columns, or brain regions are modeled by spins (i.e., with two states, up or down, on or off) with pair interactions. The emerging statistical properties of large networks of these elements are studied under different conditions (temperature or an external magnetic or electric field). The prototypical simplest system in this context is the classical 2D Ising model, which features nearest-neighbor interactions and a phase transition. This model has been shown to be universal, i.e., that all the physics of every classical spin model (with more general types of interactions) can be reproduced in the low-energy sector of certain “universal” models such as the 2D Ising model [71]. This fact reflects the intrinsic computational power of near-neighbor interactions. In fact, Ising models have been shown to be universally complete [72, 73], with a map between any given logic circuit to the ground states of some 2D Ising model Hamiltonian. On the other hand, the fact the brain exhibits characteristics of criticality that may be modeled by systems such as the Ising model is now well established, with ideas that go back to pioneers such as Turing [74], Bak [1, 75], and Hopfield [76]. There is further evidence that the dynamics of the healthy brain occupy a sub-critical zone ([77], and see also [29] and references therein).
The Ising spin model used here arises naturally from model fitting of the binarized data (discussed in the next section) and is slightly more general than the original one. It is essentially the Sherringon-Kirkpatrick spinglass model [78] with an external field and allows for arbitrary pair, arbitrary sign interactions. It is defined by the energy or Hamiltonian of the lattice of N spins,
(1)
where σi denotes the orientation of each spin in the lattice (±1), Jij is the coupling matrix or Ising connectivity (with pairs counted once), and hi is an external magnetic field applied independently at each site. In the context of analysis of parcellated BOLD data, σi in this formalism represents the binarized state of each brain parcel and is equal to + 1 (−1) when the parcel is active (inactive), hi is a parameter that modulates the mean activity in a single parcel, and
is a parameter that accounts for the interactions between parcels i and j (correlations), and which can be positive or negative.
Maximum entropy principle derivation of Ising model.
The (spinglass) Ising model, which was originally proposed to describe the physics of ferromagnetic materials, can be seen to arise naturally from the Maximum Entropy Principle (MEP) [48]. The MEP is used to find an appropriate probability distribution given some constraints, e.g., derived from data. It states that the probability of finding a given spin configuration σ can be inferred from a simplicity criterion: given some data-derived constraints—which here will be that the mean spin and spin correlation values of the model should match that of the data—the probability function we should use is the one that fits the constraints but has otherwise maximal Shannon entropy. It turns out that the probability distribution function of observing a given spin configuration σ (bold symbols indicate an array) in thermal equilibrium (at unit temperature without loss of generality) is given by the Boltzmann distribution [5]
(2)
with
given in Eq 1, and where ∑{σ′} indicates a sum over all the 2N possible spin configurations.
Estimation of archetype parameters using maximum pseudo-likelihood.
To create the global model archetype, we concatenate the binarized data from the fifteen participants and the two conditions (a total of four sessions) to generate a single dataset. We similarly produce condition (LSD and placebo) archetypes by concatenating the data in each condition.
We estimate the model parameters J and h using an approximation to maximum likelihood as described in Ezaki et al. [5]. Briefly, we find the parameters that maximize the probability of observing the data given the model,
(3)
with
(4)
The likelihood maximum can be found by gradient ascent, with
(5)
(6)
where the old/new superscripts refer to the values before and after updating, and with η a small positive constant.
When the number of nodes in the system is large, the calculation of the likelihood from the model because computationally intractable. For this reason, it is approximated by the pseudo-likelihood, a mean-field approximation that approaches the true likelihood as the number of time points goes to infinity [5, 79],
(7)
with
the modeled probability distribution for a given spin given the states of all the others, a quantity much easier to compute,
(8)
with
(9)
Using this approximation, the gradient ascent rule becomes
(10)
(11)
where
and
are the 1-point and 2-point correlation functions with respect to distribution
[5],
(12)
Personalization of model—Individualized temperature.
With h and J fixed in the archetype using the entire dataset, we adapt the archetype model for each subject and condition by changing the model temperature T = 1/β, that is, by writing
(13)
In this case, the gradient ascent algorithm becomes
(14)
with a fixed point at
. This, as in the prior equations, can be seen by taking the derivative of the approximate log-likelihood
with respect to the inverse temperature β, with the partition function, Z, given by
The first term becomes the empirical average of the Hamiltonian, and the second the (data dependent) model mean (both up to a constant since we are not dividing by the number of measurements). It can be computed from Eq 12 and from
(15)
Metropolis algorithm and Ising model observables.
The Metropolis algorithm is used to compute the observables of the Ising model built from the experimental data. A loop is performed where a random spin is chosen and flipped with a probability , where
is the change in energy the flipped spin would cause in comparison to the previous configuration and T is the system temperature. This process provides the means to sample the system probability distribution at a given temperature. The method is applied for a sufficient number of iterations to reach a steady state, and the macroscopic variables are obtained as a mean across iterations over a large number of configurations in steady state. To evaluate
, we rewrite the energy of the system (Hamiltonian) as
(16)
(17)
where in the second line, we separate out the contribution to the energy from spin σk (which receives two terms from J, which is symmetric). The energy associated with the k-th spin is
(18)
After a flip of this spin, its energy contribution becomes
(19)
and
(20)
The main global observable is the average lattice magnetization, M, over the spin lattice, S, of N spin configurations, σi,
(21)
Then, the magnetic susceptibility, χ, and heat capacity, Cv, can be computed as follows. Let us apply a uniform external field h,
(22)
and let
(23)
be the partition function. The average magnetization, 〈M〉, is the ensemble average of lattice magnetization,
(24)
Then it can be verified that (see, e.g., [80])
(25)
Following similar considerations (with T = 1/β), the global susceptibility, χ, and the system heat capacity, Cv, are given by
(26)
where the standard deviation computation is taken over the ensemble. We can extend the definition of global susceptibility (obtained from average quantities of the system) to local susceptibility, χn, (specific to each spin or brain parcel). Let
so
(27)
The global and local node susceptibility quantify, respectively, the sensitivity of the lattice state (magnetization) to an external uniform or local perturbation (represented by a magnetic field). They are interesting quantities in the context of brain dynamics from various perspectives, including the effects of external perturbations such as sensory inputs or brain stimulation. Similarly, the link susceptibility
(28)
provides information about the sensitivity to additive changes to a particular link in the system. To understand the impact of changes of scale of a particular link (Jij → Jij + αJij), we define the corresponding scaled quantity Lij = lijJij.
To smooth the results and deal with slow lattice fluctuations at cold temperatures in finite-size systems (the system is bistable for the finite case), we work with modified variables to compute the magnetization and susceptibility statistics using Metropolis-generated data. Selection of one of the metastable branches is achieved by flipping, as needed at each sample point, the entire lattice so that the magnetization is positive, σn → σn sign(M). E.g., M → |M| = M sign(M) for the evaluation of magnetization and global susceptibility. Even n-point correlation functions such as 〈Mσn〉 are invariant. The results away from the critical temperature differ by a constant from the ideal one (the effect is one of smoothing near the critical point) [81].
Statistical analysis
For the analysis of the Ising models, we performed a paired Wilcoxon test of the temperature estimates. To further test whether the observed increase in temperature under LSD relative to placebo was specific to the architecture parameters h and J, we conducted a non-parametric permutation test on the mean temperature difference [82]. In each of one thousand permutations, we randomly shuffled the parcellation labels and hence the h and J components of the archetype parameters (see Eq 3). These permuted archetypes (hp, Jp) were then used to fit the system temperature of each participant and condition (see Eq 13). The null hypothesis is that the measured BOLD activity patterns are drawn from the same probability distribution, regardless of the drug manipulation in which it was recorded. Once the personalized βp estimates were obtained, we performed the temperature contrast between LSD and placebo. This was carried out by computing a one-sample t-test for every permutation. The level of statistical significance was estimated by the fraction of t-values derived for every random shuffle that was higher than the observed t-values obtained from our original temperature contrast obtained from the non-permuted h and J. Similarly, we performed a paired one-tailed Wilcoxon test with the LZW and BDM complexity estimates.
To compare archetypes, we used paired t-tests (Wilcoxon signed-rank tests provided very similar results). We used Pearson correlation to study correlations of connectivity matrices, metrics, and receptor maps.
dMRI tractography
We will find it useful to compare Ising connectivity with structural connectivity obtained from dMRI tractography in a separate cohort as described in Cabral et al., 2014 [83]. Briefly, the data was obtained from 21 healthy, normal participants (11 males and 10 females, age: 22–45 years) different from the healthy subjects in the LSD cohort. All scans were performed on the same Philips Achieve 1.5 Tesla Magnet. Diffusion MRI was acquired by using a single-shot echo-planar imaging-based sequence with coverage of the whole brain with 33 optimal nonlinear diffusion gradient directions (b = 1200 s/mm2) and 1 non-diffusion weighted volume (b = 0), repetition time (TR) = 9390 ms; echo time (TE) = 65 ms. T1-weighted structural images with a three-dimensional ‘FLASH’ sequence (TR = 12 ms, TE = 5.6 ms, flip angle = 19°, with elliptical sampling of k-space, giving a voxel size of 111 mm in 5.05 min) were also acquired. The AAL template was used to parcellate the brain into 90 regions (45 for each hemisphere), which define the network nodes. For each participant, parcellation was conducted in the diffusion-MRI native space. The b0 image in diffusion-MRI space was linearly co-registered to the T1-weighted structural image using the Flirt tool (FMRIB, Oxford). The transformed T1-weighted image was then mapped to the T1 template of ICBM152 in MNI space, inversed and further applied to warp the AAL mask from MNI space to the diffusion-MRI native space. Interpolation using the nearest-neighbor method ensured the preservation of labeling values. The links between nodes were weighted in proportion to the number of white matter tracts detected between brain areas. The processing of the diffusion-MRI data was performed using the Fdt toolbox in FSL (www.fmrib.ox.ac.uk/fsl, FMRIB). Preprocessing involved the co-registration of the diffusion-weighted images to a reference volume using an affine transformation to correct head motion and eddy current-induced image distortion. Subsequently, the local probability distribution of fiber direction at each voxel was estimated. The probtrackx algorithm was used, allowing for automatic estimation of two fiber directions within each voxel, which can significantly improve the tracking sensitivity of non-dominant fiber populations in the human brain. For each voxel in the brain, probabilistic tractography was used to sample 5000 streamline fibers passing through that voxel. The connectivity probability from voxel i to another voxel j was defined by the proportion of fibers passing through voxel i that reach voxel j. This was then extended from the voxel level to the region (parcel) level. The connectivity Cnp from region n to region p is calculated as the number of fibers passing through any voxel in parcel n that connect to any voxel in parcel p, divided by 5000N, where N is the number of voxels in parcel n. For each brain region, the connectivity to each of the other 89 regions was calculated. Since the connectivity from n to p is not necessarily equal to that from p to n but highly correlated for all subjects, the undirected connectivity Cnp between regions n and p was defined by averaging the two.
Results
Ising connectome vs. dMRI connectome
The first step in our analysis is to identify the parameters h and J that could best fit the fMRI data from all subjects in both conditions—which leads to what we call the ‘archetype model.’ Additionally, in order to check for network changes due to the effect of LSD, we also estimated the parameters h and J that could best reproduce the data from the LSD and placebo conditions separately. Gradient descent resulted in very good fits of the functional connectivity data in each condition (r ∼ 1, see Fig J in S1 Text).
Fig 1 (top) provides the resulting connectivity matrix (i.e., the values of Jij) estimated from all the data and also the changes in connectivity between the LSD and placebo conditions. We can observe that structural intrahemispheric connectivity in the Ising model is mostly positive. In contrast, interhemispheric connectivity is mostly negative and weaker except for homotopic links (strong and positive).
A) Global archetype (combining all data) Ising connectivity. The down-trending diagonal values reflect strong inter-hemispheric homotopic connectivity. B) Difference in connectivity magnitude between the placebo and LSD archetypes, with a visible loss in homotopic connectivity (see also figures in section 5, S1 Text). C) Unsigned global archetype Ising connectivity. D) AAL connectome.
Fig 1 (bottom) provides an unsigned version of the global archetype connectivity matrix J and, for reference, the structural connectome (obtained with dMRI tractography from a separate cohort as described in Cabral et al., 2014 [83]—see Methods). We can see that the two are very similar (r = 0.60). The correlation between the two intrahemispheric submatrices is higher (r = 0.70), while for the interhemispheric ones is strong (r = 0.54) but mostly accounted for by the homotopic diagonal (removing it reduces r to 0.14).
Differences between LSD and PCB archetype connectivities
To analyze the archetypes, we introduce further metrics. We define the unsigned Global Connectivity (GC) as the array resulting from summing the rows of the magnitude (absolute value) of the connectivity adjacency matrix,
(29)
We also introduce the interhemispheric and intrahemispheric submatrices, Jinter and Jintra, where only the respective links are included (others set to zero), and their respective global connectivities GCinter and GCintra. Finally, we also use the notion of homotopic connectivity (HC). In our parcellation, we label homologous areas symmetrically (see labeling of parcels in section 1, S1 Text), and
(30)
The change in link strength between |JLSD| and |Jplacebo| was significant (statistic=-4.4, p = 1.09 × 10−5, treating the array elements in the two archetypes as samples). No significant changes were found in the signed version of global connectivity between the two archetypes, but they were strong in the unsigned version of this quantity (GC) (t-test statistic=-7.6, p = 3.1 × 10−11, Wilcoxon p = 2.3 × 10−9). The changes in h exhibited a trend to weaker external uniform drive (t-test statistic=−1.54, p = 0.13) on the system. The main effect, which can be observed in Fig 1, was a robust decrease in homotopic connectivity across parcels (connectivity between homologous areas in the two hemispheres, HC) in the LSD condition (t-test statistic=−7.8, p = 7.6 × 10−10, Wilcoxon p = 1.1 × 10−7). These changes at the level of the two archetypes reflect differences in the total homotopic FC for each subject across the two conditions (paired Wilcoxon statistic=-4.38, p = 6 × 10−4). In more detail, interhemispheric (t-test statistic=−6.9, p = 6.1 × 10−10, Wilcoxon p = 8.2 × 10−9) but not so much intrahemispheric GC (t-test statistic=−2.9, p = 0.005, Wilcoxon p = 0.003) was significantly reduced in the LSD condition. The decrease in homotopic connectivity scaled with the homotopic connectivity estimated in the placebo condition (r = −0.4, p = 1.1 × 10−5).
Archetype connectomes and 5-HT2A distribution
As the LSD effects should be related to the density of the 5-HT2A receptors [39, 84], we analyzed the correlation between the receptor density distribution, Ising connectivity, and the changes in the archetype parameters obtained separately from the LSD and placebo data. The 5-HT2A receptor map in AAL parcellation was obtained as in Deco et al. (2018) [84]. We found a weak negative correlation between homotopic connectivity and receptor density (r = -0.27, p = 0.01) and a positive correlation between GCintra and receptor density (r=0.37, p = 3 × 10−4), both especially evident in the LSD condition. For the overall GC and receptor density, the correlation was weaker (r=0.24, p = 0.02 with data from both conditions). No correlation was observed between receptor density and GCinter. Finally, we did not observe correlations between receptor density and changes in GC or homotopic connectivity.
Archetype phase diagrams
To find the critical temperature of the archetype model, we explored its phase space by calculating key quantities belonging to the Ising formalism and related measures of signal and algorithmic complexity. Figs 2 and 3 provide, respectively, the local susceptibility and local magnetization at each node in the lattice, and the global magnetization, susceptibility, energy, heat capacity, LZW, and BDM complexity estimates and their standard deviation, as a function of temperature. As can be observed, the critical point (i.e., the peak of the standard susceptibility) is well below T = 1, the nominal temperature of the archetype model (which corresponds to the global fit) and of the individualized models. This indicates that archetype and individualized models are all in the paramagnetic, supercritical phase. Susceptibility, heat capacity, and complexity peak at approximately the same temperature (T ∼ 0.7–0.8). Finally, Fig 2 displays the peak susceptibility of nodes in the parcellation and the temperatures at which they occur.
A) Local magnetization and B) local susceptibility (i.e., specific to each spin or brain parcel).
A) Global susceptibility and magnetization, B) energy and heat capacity. C) LZW complexity (ρ0) and its standard deviation, and D) BDM complexity and its standard deviation. The Metropolis algorithm has been run for 108 spin flip tests for each temperature to obtain these plots.
Individual temperature shift with LSD
The next step in our analysis was to analyze the personalized temperatures of each subject during the LSD and placebo conditions. These temperatures can be compared to the system’s critical temperature to locate where the temperatures of the two conditions and different subjects sit in the archetype phase space. Fig 4 displays, on the top, personalized temperatures of the fifteen subjects in the LSD and placebo conditions and the temperature difference between conditions for each subject.
The LZW and BDM estimates here are computed from the experimental fMRI BOLD data times series. A, C, E) Box plots of the metrics (temperature, LZW and BDM complexities) for each subject under the effect of LSD and placebo (one star: p < 0.05, four stars: p < 0.0001). B, D, F) Bar plots of the differences in metric (temperature, LZW and BDM complexities) between LSD and placebo for each subject. The dashed line in each plot represents the mean of the metric.
Temperatures in the LSD state were found to be substantially higher than in placebo condition (6.7% ± 5.1%; mean ± standard deviation), with this increase being statistically significant (two-tailed t-test p = 6.1 × 10−5; Cohen’s d = 1.35, one-tailed Wilcoxon test p = 9 × 10−5). The non-parametric permutation test shuffling ROI labels in the model to test sensitivity to archetype parameters confirmed these were significant results (p = 1 × 10−3).
LZW and BDM complexity from empirical and model data
We also extracted features related to the algorithmic complexity both from the empirical fMRI data and the synthetic data generated with the model with the Metropolis algorithm. We would expect these features to show a similar trend to the one obtained with the Ising temperatures. Figs 4 and 5 display the individual LZW and BDM complexity of 15 subjects under the effect of LSD and placebo (left) and the difference in LZW and BDM complexity between the LSD and placebo conditions for each subject (right), calculated respectively using experimental and simulated data generated by the personalized Ising models.
The LZW complexity in the LSD state from the data was found to be higher than in the placebo condition for 2/3 of the subjects (Fig 4, middle). The mean and median relative LZW complexity shift with LSD were respectively 0.008 and 0.005, with a standard deviation of 0.013. The paired one-tailed Wilcoxon test with the LZW complexity estimates resulted in p = 0.04 for the Wilcoxon test (p = 0.04 for the t-test). The probability of finding this or a larger increase of ρ0 with LSD ingestion was estimated at p < 0.04 (one-tailed Wilcoxon test).
As shown in Fig 4 (bottom), the BDM complexity in the LSD state inferred directly from the data was higher for more than 70% of subjects. The mean and median BDM complexity shift with LSD were respectively 338 and 298, with a standard deviation of 804. Also, the BDM complexity from empirical data correlated with the condition (p = 0.04, one-tailed Wilcoxon test).
As can be observed in Fig 5, the increase in LZW and BDM complexity in the LSD state, using synthetic data generated by the personalized models, was statistically significant (one-tailed Wilcoxon test p = 9 × 10−5 and p = 2 × 10−4, respectively).
Here we derive complexity estimates from synthetic data generated by the individualized archetype models for each subject. A, C: Box plots of the metrics for each subject under the effect of LSD and placebo. B, D: Bar plot of the differences in metrics between LSD and placebo for each subject. The dashed line represents the mean of the metric (three stars: p < 0.001, four stars: p < 0.0001).
We also analyzed the correlation between the features extracted from the data and from the model. A summary of the statistics of the different metrics and the Pearson correlation coefficients between them is given in Tables 1 and 2, whereas the correlation scatter plots can be found in Figure O in section 10, S1 Text. The difference in LZW and BDM complexity from the model correlates strongly with the delta Ising temperature (r(13) = 0.97, p = 3.46 × 10−9 for LZW and r(13) = 0.95, p = 1.1 × 10−7 for BDM). The changes in LZW from experimental data do not display a significant relationship with the changes in Ising temperature (r(13) = 0.06, p = 0.83), nor with the delta LZW complexity from synthetic data (r(13) = 0.10, p = 0.73), presumably due to the short length of the time series and the limited number of subjects. The difference in BDM complexity from the experimental data correlates with the delta Ising temperature (r(13) = 0.56, p = 0.03), whereas it does not show a strong statistical relationship with the delta BDM from the model (r(13) = 0.42, p = 0.12), with the delta LZW from the data (r(13) = 0.35, p = 0.19), and with the delta LZW complexity from the model (r(13) = 0.48, p = 0.07).
Change in Ising temperature and complexity metrics in the LSD condition relative to the Placebo condition.
Correlation between questionnaire scores and data features
In the LSD trial, participants were asked to answer two types of questionnaires to give feedback on their psychedelic experience. They completed VAS-style ratings via button-press and a digital display screen after each scan and the 11-factor altered states of consciousness (ASC) questionnaire at the end of each dosing day [30]. The most important ASC scores reported by the authors were elementary imagery (“I saw regular patterns in complete darkness or with closed eyes.”, “I saw colors before me in total darkness or with closed eyes.”, “I saw lights or flashes of light in total darkness or with closed eyes.”), complex imagery (“I saw scenes rolling by in total darkness or with my eyes closed.”, “I could see pictures from my past or fantasy extremely clearly.”, “My imagination was extremely vivid.”), and audio-video synesthesia (“Noises seemed to influence what I saw.”, “The shapes of things seemed to change by sounds and noises.”, and “The colors of things seemed to be changed by sounds and noises.”).
We calculated the questionnaire scores difference between the LSD and placebo sessions for each category. Then, we computed the average scores of each questionnaire to have a single measure of the overall intensity of the subjective experience. We checked the Pearson correlation between the subjects’ questionnaire scores and the features extracted from the data, i.e., Ising temperature, LZW, and BDM from experimental data.
Tables B and C in section 11, S1 Text, provide the questionnaire scores for all subjects, their average, and the Pearson correlation coefficients between the scores and the features extracted from the data.
No strong statistical relationships were seen between the subjective reports and the features. However, trends were observed between complex imagery and Ising temperatures (r = 0.42) and elementary imagery and LZW complexity (r = 0.48), both obtained from the ASC questionnaire. From the VAS-style ratings, trends were observed between emotional arousal and all three features (positive trends for Ising temperatures and BDM complexity, whereas negative ones with LZW complexity).
Discussion
We have developed here a framework for the generation and analysis of personalized Ising spin (spinglass) models starting from group archetypes. We discuss our results next.
Ising and structural connectome
A first finding from fitting the (maximum entropy) Ising model to functional BOLD data is that it is possible to recover, to a good extent, the underlying structural connectome as inferred from dMRI data. In fact, what is recovered is a signed approximation to the structural connectome, which may be useful for whole-brain modeling using more sophisticated models such as neural mass models where it is not easy to decide where long-range connectivity should be inhibitory. In addition to negative couplings, the biggest discrepancy between the two connectomes is the homotopic component of connectivity (v. Fig G in S1 Text), both of which reflect features of the FC (v. Fig K and the figures in section 7, S1 Text). This is not surprising since one of the most interesting differences between functional and structural connectomes is precisely the homotopic component, which is especially salient in resting-state fMRI functional connectivity [85] despite not being structurally as strong as detected by dMRI tractography [86, 87]. This discrepancy suggests that homotopic functional connectivity may arise from indirect links through other areas, e.g., by a synchronizing role of subcortical areas such as the thalamic nuclei, achieved by interacting symmetrically with both hemispheres [88], or from direct neuromodulation of structural elements (the tractogram provides an estimation of the number of fibers connecting parcels, but it does not provide direct information about the number of synapses associated with the fibers or the connectivity state as induced by neuromodulatory agents, for example). In either case, the Ising model must reproduce homotopic functional connectivity clearly present in the resting state data (parcel state correlations), and its structural elements (J) reflect this. An additional observation is that the Ising connectome was mostly dominated by positive intrahemispheric and homotopic links. Interhemispheric connectivity was otherwise weak and mostly negative in our dataset.
Ising archetypes for LSD and placebo conditions. 5-HT2A distribution
With regard to the archetypes built separately in each condition, the main difference between the global LSD and placebo archetypes was a decrease in connectivity under LSD, mostly associated with homotopic links. We observed the largest losses in homotopic connectivity in the calcarine, hippocampus, thalamus, and fusiform parcels (see Figs A and B in S1 Text). The increase of disorder under the LSD condition, as quantified by Ising temperature or complexity metrics, was thus associated with a loosening of interhemispheric connectivity driven by a loss of homotopic connections (Fig 1, right). The decrease in homotopic connectivity was, in turn, proportional to the placebo condition archetype homotopic connectivity. Our results suggest a transient reduction, induced by LSD, of Ising structural homotopic connectivity, where strong homotopic connections see the largest decrease in strength. Such a reduction has a strong impact on the system because the Ising connectome can be seen to consist of two networks (one per hemisphere) with weaker interhemispheric connections that are dominated by homotopic links providing “long-range” connectivity. The effect of reducing these connections can be simulated in the Ising formalism (v. Fig M in S1 Text). E.g., reducing homotopic connectivity in the placebo archetype by twenty percent leads to a decrease of its critical temperature of ten percent (from T = 0.82 to 0.74), making the system at nominal temperature (T = 1) more disordered (reducing an equal number of random connections has a smaller impact with a mean T = 0.80 with a standard deviation of 0.02). Another way to study the sensitivity of the system to changes in its connections is through the concept of linked susceptibility as defined above (Lij). In particular, this displays increased values along the homotopic line (Fig N in S1 Text). With regard to the potential mechanism for these changes, we found that the density of 5-HT2A receptors was positively correlated with the strength of intrahemispheric global connectivity and negatively correlated with the homotopic connectivity of each parcel, suggesting that the relative functional weight of homotopic links may be reduced by the increased intrahemispheric 5-HT2A agonistic effects of LSD.
Global archetype critical point and local susceptibility
We found it especially useful to generate a global archetype (v. Fig F in S1 Text) assimilating the data from all participants and conditions (see also Fig D in S1 Text for the empirical FC of all the data). This is convenient first because of the relative scarcity of data produced by BOLD fMRI sequences and, second, because it provides a common framework for comparison of individualized models (in this case, characterized by their temperatures relative to the archetype). The archetype generates behavior as a function of temperature similar to the standard, nearest-neighbor Ising model, with a clearly defined critical point. However, unlike the original Ising model, this model does not display translational symmetry because the model parameters (Jij and hi) vary with each parcel. As a consequence, some of the nodes of the archetype model display higher sensitivity to perturbations than others (see figures in section 3, S1 Text). The brain parcels corresponding to these nodes are natural areas for further study of the effects of stimulation (invasive or non-invasive) in more sophisticated computational brain models [89, 90].
Personalized temperatures and complexity
With regard to the personalized temperatures, we found a strong correlation between individualized temperature and condition. The individualized temperatures derived from BOLD fMRI data almost uniformly increased with the LSD condition relative to the placebo. Statistically speaking, the results are very robust. Moreover, both the archetype temperature (T = 1) and the individualized ones for both conditions were found to be above the critical point of the model. The state in all cases was found to be in the paramagnetic phase, as was found in Ezaki et al. (2020) [6].
Since Ising temperature was expected to correlate with entropy and disorder in the models, and because of the relevance of the notion of complexity in the analysis of brain data, we studied two proxies (upper bounds) of algorithmic complexity, one closely related to Shannon entropy (LZW), the other more directly connected with the concept of Kolmogorov complexity (BDM). We did this first, starting from the data itself and then also working with synthetic data generated from the personalized Ising models. Complexity metrics derived from the data related only weakly (but significantly at the group level) with temperature and condition, presumably due to the limited amount of data per participant and condition, and we cannot discount the possibility of paradoxical effects of LSD in some subjects leading to decreases of complexity. However, it is well known that measures that estimate complexity are sensitive to data length [66, 91], and this is a limitation in this study because BOLD data is intrinsically slowly varying and the time series are not very long. To compensate for this, we explored model-derived LZW and BDM complexity estimates using synthetic data. We found these to be monotonic functions of temperature (Figs 3 and 5) that correlated strongly with the condition.
Fig 3 (bottom right) displays the monotonic behavior of BDM complexity with respect to temperature, with a very sharp increase and then almost a plateau. One limitation in the BDM calculation is that, as the size of matrices increases, BDM gets gradually closer and closer to an entropy measurement as we use precomputed CTM for smaller matrices. Moreover, it is important to note that, in our Ising and complexity analysis, we are bounded to binary data, which makes the analysis contingent on how meaningful the binarized version of the data is.
Overall, our results are consistent with the notion that psychedelics drive brain dynamics into a more disordered state and, in our modeling framework, away from criticality. The RElaxed Beliefs Under pSychedelics (REBUS) framework [39]) addresses the phenomenology and physiology of psychedelic states by integrating mechanistic insights with the free-energy principle (FEP, [92, 93]). The FEP is a formulation based on causal statistical theory that accounts for action, perception, and learning, and the entropic brain hypothesis (EBH [29]), which associates entropy of spontaneous brain activity with the richness (vividness and diversity) of subjective experience (and posits that psychedelics increase both).
Relation to the REBUS model and criticality
The mechanisms of action of psychedelics are believed to begin with stimulation of a specific serotonin receptor subtype (5-HT2A) on cortical layer V pyramidal neurons, i.e., this is the control signal that results in increased neuronal excitability and dysregulated cortical population-level activity [94]. Through this action, psychedelics are believed to disrupt the regular functioning of various high-level system properties—e.g., major cortical oscillatory rhythms and the integrity of large-scale networks—that are hypothesized to encode the precision-weighting of internal models—i.e., priors, beliefs, or assumptions. According to the so-called ‘REBUS’ model [39], psychedelics disrupt models sitting atop the modeling hierarchy, with direct consequences on experience. Interference with microcircuitry associated with high-level models may be expected to release lower-level models that would otherwise be suppressed.
In [40], we hypothesized that increased apparent (Shannon) complexity (entropy) of spontaneous activity is a logical corollary of this action, as apparently complex, bottom-up uncompressed data more strongly dominates brain activity. Such increases have indeed been observed using measures such as LZ [95] and are also seen in our data. The present study extends on previous work, however, by incorporating algorithmically-oriented methods such as BDM. Our results agree with earlier work pointing to the usefulness of complexity metrics and the hypothesis that signal complexity should increase with increased system temperature, and that psychedelics should increase both [50]. With regard to the criticality perspective, in Ruffini & Lopez-Sola (2022), we proposed that a modeling system tracking structured data in the world (a wakeful brain) will display dimensionality reduction and criticality features that can be used empirically to quantify the structure of the program run by an agent [40]. Criticality features are seen to result from the collapse of a dynamical system (the dynamical brain) into a low dimensional space associated with identified regularities in complex data—the hallmark of an algorithmic compressive system. If tracking world data using models keeps the modeling dynamics on a reduced, lower-dimensional manifold, psychedelics may be said to “lift” up the enslaved dynamics into higher dimensions. Our finding of the resting brain state being above criticality and driven further away from the critical point under the influence of LSD is consistent with this view.
As discussed in the context of the algorithmic information theory of consciousness (KT) [40, 50], psychedelics will shift the dynamics of an agent system tracking world data to a less constrained state further away from criticality, and hence produce more complex signals. The above is in line with the idea that the brain operates near a critical point (v. [21, 96, 97] and references therein) and that psychedelics move dynamics into a more disordered state [39]. However, at least in our model and as already described in previous work using similar methods [5, 6], we found that the resting brain in the placebo condition was already above the critical point—that is, the resting, the wakeful brain is in a supercritical state as observed through the Ising framework lens. This is consistent with the idea that, while mutual information peaks at the critical temperature, information flow in such systems peaks in the disordered (paramagnetic) phase [98], which is in apparent contradiction with other studies suggesting a subcritical nature of healthy brain dynamics [39]. We note here in passing that the term “criticality” is used differently by different authors, e.g., dynamical vs. statistical criticality, whose relation is the subject of current research [97]. At least with the data in this study, the statistical criticality metric developed here—the supercritical system temperature of an Ising model—, provided the best correlate of experimental condition.
Finally, with regard to the relation of these metrics with subjective experience as measured by questionnaires, some weak correlations were found at the trend level with all metrics (temperature and complexity metrics), presumably because of the limited number of data and the inherent heterogeneity of self-reports and that these particular ratings were done post-hoc without precise reference to the scanning periods. Specifically designed neurophenomenological methods probing the structure of experience during altered states of consciousness may provide better correlates in future studies [40].
Conclusions
In this paper, we have studied criticality features of brain dynamics under LSD in a within-subject study using the Ising model formalism and algorithmic complexity using Lempel-Ziv and the Block Decomposition methods. Personalized Ising models were created using BOLD fMRI data by adjusting the model temperature for each subject and condition in relation to the group average model. We found that the effects of LSD translate into increased BOLD signal complexity and Ising temperature, in agreement with earlier findings and predictions from existing theories of the effects of psychedelics, such as the relaxed beliefs under psychedelics (REBUS), the anarchic brain hypothesis [39], and the algorithmic information theory of consciousness (KT) [40, 50]. More precisely, in agreement with earlier work [6], we found that the system in the placebo condition was already in the paramagnetic phase—above the critical point—, with ingestion of LSD resulting in a shift away from criticality into a more disordered state. These results must be interpreted with care, as the Ising framework used here provides a simple one-parameter model for fitting the individual data. Other studies using different methods suggest that LSD shifts system dynamics closer to criticality [99, 100], while, for example, using similar methods to ours, higher fluid intelligence has been associated with closer to critical resting-state neural dynamics [6]. This highlights the question of whether modeling choices such as what data and data features to model (e.g., BOLD FC or electrophysiological data), modeling scale (micro, meso or macroscale), and how to model it (e.g., Ising, Hopf, or neural mass models) can impact conclusions on the critical features of the system. Further work is needed to reconcile these results, which may stem from both these modeling choices and associated criticality metrics [101]. Finally, the Ising formalism provides a parsimonious method to infer structural connectivity from functional data. The Ising connectome, we found, is related to the structural (dMRI) connectome but differs in two important ways. First, it displays stronger homotopic connectivity, which, as we saw, is weakened during LSD. Second, it is a signed connectome, and the negative connections play an important dynamical role. Thus, Ising connectomics may provide complementary insights into the general problem of whole brain computational modeling.
Supporting information
S1 Text. Supplementary information file.
This file contains figures and tables in the following sections: 1) List of AAL areas, 2) Changes in connectivity, 3) Susceptibility per parcel, 4) Empirical FC, 5) Ising and dMRI connectomes, 6) Network plots, 7) Group FC for each condition (empirical vs. Ising), 8) Homotopic FC and FC negativity. FC for each subject and condition, 9) Homotopic and link susceptibility, 10) Correlation between metrics, 11) Subject questionnaire scores and their correlation with metrics.
https://doi.org/10.1371/journal.pcbi.1010811.s001
(PDF)
References
- 1. Chialvo DR. Critical brain networks. Physica A. 2004;340:756–765.
- 2. Chialvo DR. Emergent complex neural dynamics. Nat Phys. 2010;6(10):744–750.
- 3. Das TK, Abeyasinghe PM, Crone JS, Sosnowski A, Laureys S, Owen AM, et al. Highlighting the Structure-Function Relationship of the Brain with the Ising Model and Graph Theory. BioMed Research International. 2014;(237898). pmid:25276772
- 4. Cocchi L, Gollo LL, Zalesky A, Breakspear M. Criticality in the brain: A synthesis of neurobiology, models and cognition. Prog Neurobiol. 2017;158:132–152. pmid:28734836
- 5. Ezaki T, Watanabe T, Ohzeki M, Masuda N. Energy landscape analysis of neuroimaging data. Philos Trans R Soc A. 2017;375(20160287). pmid:28507232
- 6. Ezaki T, Fonseca Dos Reis E, Watanabe T, Sakaki M, Masuda N. Closer to critical resting-state neural dynamics in individuals with higher fluid intelligence. Commun Biol. 2020;3(1):52. pmid:32015402
- 7. Fraiman D, Balenzuela P, Foss J, Chialvo DR. Ising-like dynamics in large-scale functional brain networks. Phys Rev E Stat Nonlin Soft Matter Phys. 2009;79(6 Pt 1):061922. pmid:19658539
- 8.
Niebur E, Plenz D, editors. Criticality in Neural Systems. Annual Reviews of Nonlinear Dynamics and Complexity (VCH). Weinheim, Germany: Wiley-VCH Verlag; 2014.
- 9. Yu S, Yang H, Shriki O, Plenz D. Universal organization of resting brain activity at the thermodynamic critical point. Front Syst Neurosci. 2013;7:42. pmid:23986660
- 10.
Lewin R. Complexity. 2nd ed. Chicago, IL: University of Chicago Press; 1992.
- 11.
Goldenfeld N. Lectures on Phase Transitions and the Renormalization Group. London, England: Addison Wesley; 1991.
- 12. Shew WL, Plenz D. The functional benefits of criticality in the cortex. Neuroscientist. 2013;19(1):88–100. pmid:22627091
- 13. Mediano PAM, Shanahan M. Balanced information storage and transfer in modular spiking neural networks. arxiv. 2017;.
- 14. Erten E, Lizier J, Piraveenan M, Prokopenko M. Criticality and information dynamics in epidemiological models. Entropy (Basel). 2017;19(5):194.
- 15. Li M, Han Y, Aburn MJ, Breakspear M, Poldrack RA, Shine JM, et al. Transitions in information processing dynamics at the whole-brain network level are driven by alterations in neural gain. PLOS Computational Biology. 2019;15(10):1–27. pmid:31613882
- 16. Mediano PAM, Rosas FE, Farah JC, Shanahan M, Bor D, Barrett AB. Integrated information as a common signature of dynamical and information-processing complexity. Chaos. 2022;32(1):13115. pmid:35105139
- 17. Deco G, Jirsa VK. Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J Neurosci. 2012;32(10):3366–3375. pmid:22399758
- 18. Deco G, Jirsa V, McIntosh AR, Sporns O, Kötter R. Key role of coupling, delay, and noise in resting brain fluctuations. Proceedings of the National Academy of Sciences. 2009;106(25):10302–10307. pmid:19497858
- 19. Cieri F, Zhuang X, Caldwell JZK, Cordes D. Brain Entropy During Aging Through a Free Energy Principle Approach. Frontiers in Human Neuroscience. 2021;15:139. pmid:33828471
- 20. Tagliazucchi E, Balenzuela P, Fraiman D, Chialvo DR. Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis. Front Physiol. 2012;3:15. pmid:22347863
- 21. Haimovici A, Tagliazucchi E, Balenzuela P, Chialvo DR. Brain organization into resting state networks emerges at criticality on a model of the human connectome. Phys Rev Lett. 2013;110(17):178101. pmid:23679783
- 22. Marinazzo D, Pellicoro M, Wu G, Angelini L, Cortés JM, Stramaglia S. Information transfer and criticality in the Ising model on the human connectome. PLoS One. 2014;9(4):e93616. pmid:24705627
- 23. Hancock F, Rosas F, Mediano P, Luppi A, Cabral J, Dipasquale O, et al. May the 4c’s be with you: An overview of complexity-inspired frameworks for analyzing resting-state neuroimaging data; 2022.
- 24.
Jensen HJ. Self-Organized Criticality: Emergent Complex Behavior in Physical and Biological Systems. Cambridge Lecture Notes in Physics. Cambridge University Press; 1998.
- 25.
Pruessner G. Self-Organised Criticality. Cambridge, England: Cambridge University Press; 2012.
- 26.
Christensen K, Moloney NR. Complexity and criticality. London, England: Imperial College Press; 2005.
- 27. Plenz D, Ribeiro TL, Miller SR, Kells PA, Vakili A, Capek EL. Self-Organized Criticality in the Brain. Frontiers in Physics. 2021;9.
- 28. Hesse J, Gross T. Self-organized criticality as a fundamental property of neural systems. Frontiers in Systems Neuroscience. 2014;8. pmid:25294989
- 29. Carhart-Harris R. The entropic brain—revisited. Neuropharmacology. 2018;142(167–178). pmid:29548884
- 30. Carhart-Harris RL, Muthukumaraswamy S, Roseman L, Kaelen M, Droog W, Murphy K, et al. Neural correlates of the LSD experience revealed by multimodal neuroimaging. Proceedings of the National Academy of Sciences. 2016;113(17):4853–4858. pmid:27071089
- 31. Gasser P, Kirchner K, Passie T. LSD-assisted psychotherapy for anxiety associated with a life-threatening disease: a qualitative study of acute and sustained subjective effects. J Psychopharmacol. 2015;29(1):57–68. pmid:25389218
- 32. Griffiths RR, Johnson MW, Carducci MA, Umbricht A, Richards WA, Richards BD, et al. Psilocybin produces substantial and sustained decreases in depression and anxiety in patients with life-threatening cancer: A randomized double-blind trial. J Psychopharmacol. 2016;30(12):1181–1197. pmid:27909165
- 33. Johnson MW, Hendricks PS, Barrett FS, Griffiths RR. Classic psychedelics: An integrative review of epidemiology, therapeutics, mystical experience, and brain network function. Pharmacol Ther. 2019;197:83–102. pmid:30521880
- 34. Krebs TS, Johansen PØ. Lysergic acid diethylamide (LSD) for alcoholism: meta-analysis of randomized controlled trials. J Psychopharmacol. 2012;26(7):994–1002. pmid:22406913
- 35. McCall MA, Tieman DG, Hirsch HV. Chronic intraventricular administration of lysergic acid diethylamide (LSD) affects the sensitivity of cortical cells to monocular deprivation. Brain Res. 1982;250(2):301–308. pmid:7171990
- 36. Passie T, Halpern JH, Stichtenoth DO, Emrich HM, Hintzen A. The pharmacology of lysergic acid diethylamide: a review. CNS Neurosci Ther. 2008;14(4):295–314. pmid:19040555
- 37. Müller F, Lenz C, Dolder P, Lang U, Schmidt A, Liechti M, et al. Increased thalamic resting-state connectivity as a core driver of LSD-induced hallucinations. Acta Psychiatr Scand. 2017;136(6):648–657. pmid:28940312
- 38. Preller KH, Razi A, Zeidman P, Stämpfli P, Friston KJ, Vollenweider FX. Effective connectivity changes in LSD-induced altered states of consciousness in humans. Proc Natl Acad Sci U S A. 2019;116(7):2743–2748. pmid:30692255
- 39. Carhart-Harris RL, Friston KJ. REBUS and the Anarchic Brain: Toward a Unified Model of the Brain Action of Psychedelics. Pharmacological Reviews. 2019;71(3):316–344. pmid:31221820
- 40. Ruffini G, Lopez-Sola E. AIT foundations of structured experience. J AI Consci. 2022; p. 1–39.
- 41. Geyer MA, Vollenweider FX. Serotonin research: contributions to understanding psychoses. Trends Pharmacol Sci. 2008;29(9):445–453. pmid:19086254
- 42. Carhart-Harris RL, Erritzoe D, Williams T, Stone JM, Reed LJ, Colasanti A, et al. Neural correlates of the psychedelic state as determined by fMRI studies with psilocybin. Proc Natl Acad Sci U S A. 2012;109(6):2138–2143. pmid:22308440
- 43. Carhart-Harris RL, Leech R, Hellyer PJ, Shanahan M, Feilding A, Tagliazucchi E, et al. The entropic brain: a theory of conscious states informed by neuroimaging research with psychedelic drugs. Front Hum Neurosci. 2014;8:20. pmid:24550805
- 44. Singleton SP, Luppi AI, Carhart-Harris RL, Cruzat J, Roseman L, Nutt DJ, et al. LSD and psilocybin flatten the brain’s energy landscape: insights from receptor-informed network control theory; 2021.
- 45. Preller KH, Burt JB, Ji JL, Schleifer CH, Adkinson BD, Stämpfli P, et al. Changes in global and thalamic brain connectivity in LSD-induced altered states of consciousness are attributable to the 5-HT2A receptor. Elife. 2018;7. pmid:30355445
- 46. Tagliazucchi E, Roseman L, Kaelen M, Orban C, Muthukumaraswamy SD, Murphy K, et al. Increased global functional connectivity correlates with LSD-induced ego dissolution. Curr Biol. 2016;26(8):1043–1050. pmid:27085214
- 47. McCulloch DEW, Knudsen GM, Barrett FS, Doss MK, Carhart-Harris RL, Rosas FE, et al. Psychedelic resting-state neuroimaging: A review and perspective on balancing replication and novel analyses. Neurosci Biobehav Rev. 2022;(104689):104689. pmid:35588933
- 48. Jaynes ET. Information theory and statistical mechanics. Physical Review. 1957;106(4):620–630.
- 49.
Jaynes ET. Probability theory—the logic of science. Cambridge; 2003.
- 50. Ruffini G. An algorithmic information theory of consciousness. Neuroscience of Consciousness. 2017;3(1). pmid:30042851
- 51. Ruffini G, Deco G. The 2D Ising model, criticality and AIT; 2021.
- 52. Deco G, Cruzat J, Cabral J, Knudsen G, Carhart-Harris R, Whybrow P, et al. Whole-Brain Multimodal Neuroimaging Model Using Serotonin Receptor Maps Explains Non-linear Functional Effects of LSD. Curr Biol. 2018;28(19):3065–3074. pmid:30270185
- 53. Tzourio-Mazoyer N, Landeau B, Papathanassiou D, Crivello F, Etard O, Delcroix N, et al. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage. 2002;15(1):273–289. pmid:11771995
- 54. Lempel A, Ziv J. On the complexity of finite sequences. IEEE Transactions on Information Theory. 1976;IT-22(1):75–81.
- 55. Welch T. A Technique for High-Performance Data Compression. Computer. 1984;17(6):8–18.
- 56.
Cover TM, Thomas JA. Elements of information theory. 2nd ed. John Wiley & sons; 2006.
- 57.
Gomez C, Hornero R, Abasolo D, Lopez M, Fernandez A. Decreased Lempel-Ziv complexity in Alzheimer’s disease patients’ magnetoencephalograms. In: 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference. IEEE; 2005.
- 58. Abásolo D, Hornero R, Gómez C, García M, López M. Analysis of EEG background activity in Alzheimer’s disease patients with Lempel-Ziv complexity and central tendency measure. Med Eng Phys. 2006;28(4):315–322. pmid:16122963
- 59.
Zhao P, Van-Eetvelt P, Goh C, Hudson N, Wimalaratna S, Ifeachor E. Characterization of EEGs in Alzheimer’s disease using information theoretic methods. In: 2007 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE; 2007.
- 60. Ruffini G, Ibañez D, Kroupi E, Gagnon JF, Montplaisir J, Postuma RB, et al. Algorithmic complexity of EEG for prognosis of neurodegeneration in idiopathic rapid eye movement behavior disorder (RBD). Ann Biomed Eng. 2019;47(1):282–296. pmid:30167913
- 61. Hudetz AG, Liu X, Pillay S, Boly M, Tononi G. Propofol anesthesia reduces Lempel-Ziv complexity of spontaneous brain activity in rats. Neurosci Lett. 2016;628:132–135. pmid:27291459
- 62. Casali AG, Gosseries O, Rosanova M, Boly M, Sarasso S, Casali KR, et al. A Theoretically Based Index of Consciousness Independent of Sensory Processing and Behavior. Sci Transl Med. 2013;5(198). pmid:23946194
- 63. Tononi G, Boly M, Massimini M, Koch C. Integrated information theory: from consciousness to its physical substrate. Nature Reviews Neuroscience. 2016;17:450–461. pmid:27225071
- 64. Moser J, Bensaid S, Kroupi E, Schleger F, Wendling F, Ruffini G, et al. Evaluating complexity of fetal MEG signals: A comparison of different metrics and their applicability. Front Syst Neurosci. 2019;13:23. pmid:31191264
- 65. Zhang XS, Roy RJ, Jensen EW. EEG Complexity as a Measure of Depth of Anesthesia for Patients. IEEE Transactions on Biomedical Engineering. 2001;48(12):1424–1433. pmid:11759923
- 66. Kaspar F, Schuster HG. Easily calculable measure for the complexity of spatiotemporal patterns. Phys Rev A. 1987;36(2):842–848. pmid:9898930
- 67. Ruffini G, Ibañez D, Castellano M, Dunne S, Soria-Frisch A. EEG-driven RNN classification for prognosis of neurodegeneration in at-risk patients. ICANN 2016. 2016;.
- 68.
Ruffini G. Lempel-Ziv Complexity Reference. arXiv:170709848 [csIT]. 2017;.
- 69. Zenil H, Hernández-Orozco S, Kiani NA, Soler-Toscano F, Rueda-Toicen A, Tegnér J. A Decomposition Method for Global Evaluation of Shannon Entropy and Local Estimations of Algorithmic Complexity. Entropy. 2018;20(8). pmid:33265694
- 70. Schneidman E, II MJB, Segev R, Bialek W. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature. 2006;440(20):1007–10112. pmid:16625187
- 71. las Cuevas GD, Cubitt TS. Simple universal models capture all classical spin physics. Science. 2016;351 (6278).
- 72. Gu M, Perales A. Encoding universal computation in the ground states of Ising lattices. Phys Rev E. 2012;86(011116):1–6. pmid:23005377
- 73. Kohler T, Cubitt T. Translationally Invariant Universal Classical Hamiltonians. Journal of Statistical Physics. 2019;176:228–261.
- 74. Turing AM. I.—computing machinery and intelligence. Mind. 1950;LIX(236):433–460.
- 75. Bak P, Tang C, Wiesenfeld K. Self-organized criticality: An explanation of the 1/f noise. Phys Rev Lett. 1987;59(4):381–384. pmid:10035754
- 76. Hopfield JJ. Neural networks and physical systems with emergent collective computational abilities. PNAS. 1982;79:2554–2558. pmid:6953413
- 77. Priesemann V, Wibral M, Valderrama M, Pröpper R, Le Van Quyen M, Geisel T, et al. Spike avalanches in vivo suggest a driven, slightly subcritical brain state. Front Syst Neurosci. 2014;8:108. pmid:25009473
- 78. Sherrington D, Kirkpatrick S. Solvable model of a spin-glass. Phys Rev Lett. 1975;35(26):1792–1796.
- 79. Besag J. Statistical analysis of non-lattice data. J R Stat Soc D. 1975;25:179–195.
- 80.
Pathria RK, Beale PD. Statistical Mechanics. Elsevier Ltd; 2011.
- 81. Ibarra-García-Padilla E, Malanche-Flores CG, Poveda-Cuevas FJ. The hobbyhorse of magnetic systems: the Ising model. Eur J Phys. 2016;37(6):065103.
- 82. Holmes AP, Blair R, Watson J, Ford I. Nonparametric analysis of statistic images from functional mapping experiments. Journal of Cerebral Blood Flow & Metabolism. 1996;16(1):7–22. pmid:8530558
- 83. Cabral J, Luckhoo H, Woolrich M, Joensson M, Mohseni H, Baker A, et al. Exploring mechanisms of spontaneous functional connectivity in MEG: How delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations. NeuroImage. 2014;90:423–435. pmid:24321555
- 84. Deco G, Cabral J, Saenger VM, Boly M, Tagliazucchi E, Laufs H, et al. Perturbation of whole-brain dynamics in silico reveals mechanistic differences between brain states. NeuroImage. 2018;169:46—56. pmid:29225066
- 85. Toro R, Fox PT, Paus T. Functional coactivation map of the human brain. Cereb Cortex. 2008;18(11):2553–2559. pmid:18296434
- 86. Rosen BQ, Halgren E. A whole-cortex probabilistic diffusion tractography connectome. eNeuro. 2021;8(1):ENEURO.0416–20.2020. pmid:33483325
- 87. Rosen BQ, Halgren E. An estimation of the absolute number of axons indicates that human cortical areas are sparsely connected. PLoS Biol. 2022;20(3):e3001575. pmid:35286306
- 88. Mancuso L, Costa T, Nani A, Manuello J, Liloia D, Gelmini G, et al. The homotopic connectivity of the functional brain: a meta-analytic approach. Sci Rep. 2019;9(1):3346. pmid:30833662
- 89. Deco G, Kringelbach ML. The economy of anatomy: Discovering the turbulent homogeneous isotropic functional core organisation of the human brain. bioRxiv. 2019;.
- 90. Ruffini G, Wendling F, Sanchez-Todo R, Santarnecchi E. Targeting brain networks with multichannel transcranial current stimulation (tCS). Current Opinion in Biomedical Engineering. 2018;.
- 91. Zhang XS, Roy RJ, Jensen EW. EEG complexity as a measure of depth of anesthesia for patients. IEEE Trans Biomed Eng. 2001;48(12):1424–1433. pmid:11759923
- 92. Friston K. The free-energy principle: a unified brain theory? Nat Rev Neurosci. 2010;11(2):127–138. pmid:20068583
- 93.
Parr T. Active inference: the free energy principle in mind, brain, and behavior. Cambridge, Massachusetts: The MIT Press; 2022.
- 94. Nutt D, Erritzoe D, Carhart-Harris R. Psychedelic psychiatry’s brave new world. Cell. 2020;181(1):24–28. pmid:32243793
- 95. Schartner MM, Carhart-Harris RL, Barrett AB, Seth AK, Muthukumaraswamy SD. Increased spontaneous MEG signal diversity for psychoactive doses of ketamine, LSD and psilocybin. Scientific Reports. 2017;7(46421). pmid:28422113
- 96. R Chialvo D. Critical brain networks. Physica A. 2004;340(4):756–765.
- 97. Zimmern V. Why brain criticality is clinically relevant: A scoping review. Front Neural Circuits. 2020;14:54. pmid:32982698
- 98. Barnett L, Lizier JT, Harré M, Seth AK, Bossomaier T. Information flow in a kinetic Ising model peaks in the disordered phase. Phys Rev Lett. 2013;111(17):177203. pmid:24206517
- 99. Atasoy S, Roseman L, Kaelen M, Kringelbach ML, Deco G, Carhart-Harris RL. Connectome-harmonic decomposition of human brain activity reveals dynamical repertoire re-organization under LSD. Sci Rep. 2017;7(1):17661. pmid:29247209
- 100. Varley TF, Carhart-Harris R, Roseman L, Menon DK, Stamatakis EA. Serotonergic psychedelics LSD & psilocybin increase the fractal dimension of cortical brain activity in spatial and temporal domains. Neuroimage. 2020;220(117049):117049. pmid:32619708
- 101. Toker D, Pappas I, Lendner JD, Frohlich J, Mateos DM, Muthukumaraswamy S, et al. Consciousness is supported by near-critical slow cortical electrodynamics. Proc Natl Acad Sci U S A. 2022;119(7):e2024455119. pmid:35145021