Paths to Oblivion: Common Neural Mechanisms of Anaesthesia and Disorders of Consciousness

The human brain generates a rich repertoire of spatiotemporal dynamics during normal wakefulness, supporting a wide variety of conscious experiences and cognitive functions. However, neural dynamics are reconfigured, in comparable ways, when consciousness is lost either due to anaesthesia or disorders of consciousness (DOC). Here, leveraging a neurobiologically realistic whole-brain computational model informed by functional MRI, diffusion MRI, and PET, we sought to identify the neurobiological mechanisms that explain the common reconfiguration of neural dynamics observed both for transient pharmacological intervention and chronic neuroanatomical injury. Our results show that, by incorporating local inhibitory action through a PET-based GABA receptor density map, our model can reproduce the brain dynamics of subjects undergoing propofol anaesthesia, and that this effect depends specifically on the spatial distribution of GABA receptors across cortical regions. Additionally, using a structural connectome obtained from DOC patients, we demonstrate how the dynamics that characterise loss of consciousness can emerge from changes in neuroanatomical connectivity. Crucially, we find that each of these two interventions generalises across datasets: a model with increased GABA-mediated inhibition can reproduce the dynamics of DOC patients’ brains, and a model with a DOC connectome is also compatible with brain dynamics observed during propofol anaesthesia. These results demonstrate how increased inhibition and connectome randomisation represent different neurobiological paths towards the characteristic dynamics of the unconscious brain. Overall, the present findings begin to disentangle the neurobiological mechanisms by which highly dissimilar perturbations of the brain’s neurodynamics can lead to unconsciousness.


53
A central challenge of contemporary neuroscience is the quest to understand how the 54 neurobiology and function of the human brain give rise to conscious experience 1,2 . One way 55 to address this question is to identify changes in brain function that accompany changes in 56 conscious state. However, the brain is a paradigmatic example of a complex system 3 , and In order to shed light on how the human brain supports consciousness, two fundamental 64 challenges must be addressed. The first is to identify neurobiological signatures of 65 consciousness and its loss that are generalisable across different paths to unconsciousness, 66 rather than being specific to any of them 13 . The second is to provide a mechanistic  (Figure 1). 143 We used this modelling approach to simulate the empirical fMRI macroscale brain dynamics 144 observed in the same N=16 subjects at baseline and during loss of consciousness induced by 145 the intravenous anaesthetic, propofol. We also studied the fMRI dynamics of a cohort (N=21) 146 of patients suffering from chronic disorders of consciousness (DOC) as a result of severe brain 147 injury (traumatic or anoxic), comparing them with a group of N=20 healthy controls. By

164
The neurobiologically realistic dynamic mean-field (DMF) model only one free parameter: a 165 global coupling parameter, denoted by G, which scales the excitatory-to-excitatory coupling 166 between brain regions, as established by the empirical structural connectome. Thus, calibrating 167 the model corresponds to finding the value of G that allows the model to best simulate observed 168 fMRI dynamics of the human brain at rest (which is assessed in terms of the KS-distance 169 between real and simulated data, see 86 ). We followed this procedure for each of our two

196
Inhibitory modulation from GABA-A receptor distribution reveals a 197 shared mechanism for loss of consciousness 198 The effects of propofol anaesthesia on the brain were modelled by capitalising on the recently 199 built whole-brain map of GABA-A regional receptor density, generated on the basis of  approach, here we show for the first time that the influence of regional GABA-A receptor 209 density on functional dynamics can be modelled using a DMF model informed by regional 210 GABA-A receptor density.

212
The strategy followed in 86 was to first calibrate the model on baseline data to obtain a global 213 coupling value, and then fit a secondary excitatory parameter separately on baseline and post-214 dose data. Our approach follows Deco's but differs in one key respect: given the inhibitory 215 nature of GABA, we modulated the inhibitory (rather than excitatory) local gain. 1 To do so, we 216 introduced an inhibitory gain scaling parameter in the model, denoted by sI. This parameter 217 allowed us to scale the inhibitory gain at each region according to the empirical local density 218 of GABA-A receptors, as quantified based on PET-derived maps of receptor density 94 .

220
This procedure allowed us to ask whether adjusting the value of inhibitory gain sI according to 221 local GABA-A receptor density would allow the model to simulate the characteristic dynamics 222 of acute propofol-induced unconsciousness. A positive answer to this question would implicate 223 regional GABA-ergic inhibition as a neurobiological mechanism behind the action of propofol 224 (a known GABA-ergic agonist) in inducing the characteristic macroscale dynamics observed 225 during loss of consciousness due to propofol anaesthesia 13 information about the regional distribution of GABA-A receptors) might produce generally 240 better-fitting simulations just in virtue of its increased complexity. To control for this 241 possibility, we compared the propofol model with an analogous model, which also incorporated 242 regional GABA-A receptor density, but whose inhibitory gain scaling parameter was chosen 243 as the one that best fitted the empirical dynamics observed in awake subjects. We refer to this 244 as the baseline model.

246
Having completed the fitting procedure for our models, we then proceeded to analyse the 247 models' performance. To this end, we generated 100 simulations from each of the propofol and 248 baseline models. For both models, we then computed the KS distance between each simulation, 249 and the empirical dynamics observed during wakefulness, and during anaesthesia. This 250 provided us with a way to quantify the ability of each model (in terms of goodness of fit, i.e. 251 low KS distance) to simulate the empirical brain dynamics observed during wakefulness, and 252 the empirical brain dynamics observed during propofol-induced loss of consciousness.

254
Results show that incorporating the regional distribution of GABA-A receptors as local 255 modulators of the global inhibitory gain, substantially improved the DMF model's fit to 256 empirical brain dynamics observed during propofol anaesthesia ( Figure 3B and Table 1). This 257 improvement could not be explained as a generic effect due to additional model complexity: if 258 this had been the case, improvements should have been observed for the model's ability to fit 259 both awake and anaesthetised dynamics. Instead, the improvement was specific to 260 anaesthetised dynamics, and the propofol model actually had poorer ability to fit the awake 261 dynamics. In other words, the addition of this inhibition parameter in accordance with its 262 empirical distribution across brain regions, makes the model capable of switching between 263 simulating awake or anaesthetised brain dynamics -and when it is simulating anaesthetised 264 brain dynamics, the model's fit to awake brain dynamics deteriorates. Since propofol is a well-265 known GABA-ergic agonist, these results confirm that taking into account GABA agonism 266 (local modulation of inhibitory gain by regional GABA-A receptor density) is sufficient to 267 recapitulate the known effects of the GABA-ergic agent propofol on empirical brain dynamics, 268 leading to dynamics that are known to characterise the state of unconsciousness.    Crucially, we also confirmed that the improved fit to anaesthetised dynamics is not merely the  Figure 4B). In both cases, the model's ability to fit anaesthetised dynamics is 296 not improved with respect to the baseline model ( Figure 4A,B and Table 1). Therefore, our 297 results show that the specific regional distribution of GABA-A receptors across the cortex plays    Whole-brain computational models provide a unique tool to understand the effects of 316 connectome alterations on macroscale brain dynamics 70,71,83 . Connectome replacement allows 317 us to determine which of two conditions is more compatible with a given perturbation of the 318 connectome, in terms of the connectome's capacity to support the corresponding brain 319 dynamics. Specifically, a smaller impact of connectome perturbation on the model's fit to 320 condition X than condition Y, indicates that the perturbed connectome is better suited to 321 supporting the dynamics of condition X than Y. We term this procedure "Connectome 322 Replacement Analysis", which is illustrated in Figure 5A.

323
Leveraging this capability, we subjected the DMF model to the virtual equivalent of severe 324 brain injury: namely, we replaced the underlying connectivity matrix governing the long-range 325 interactions between brain regions with a connectome obtained from diffusion-weighted 326 imaging of N=21 patients with chronic DOC due to severe brain injury. This procedure imparts 327 on the model with effects akin to what severe brain injury does on anatomical connectivity.

328
This "virtual DOC" provides a way to isolate the effects over brain dynamics of connectivity 329 disruptions that result in loss of human consciousness.

330
Note that such substantial perturbations are expected to deteriorate the model's ability to 331 replicate the dynamics of awake healthy brains (i.e. increasing the KS distance, corresponding 332 to decreased goodness-of-fit): our models were optimised with biophysical parameters 333 pertaining to healthy brains, and using a healthy connectome. However, our hypothesis was 334 that the dynamics generated by the model with DOC connectome should be more similar (lower 335 KS distance, indicating a better fit) to the empirical dynamics of DOC patients' brains, than to 336 the dynamics of conscious, healthy brains.

337
Our results supported these predictions. As expected, connectome replacement led to a 338 reduction in the model's ability to fit control brain dynamics ( Figure 5B). Also as expected, the 339 model with the healthy connectome was better able to simulate conscious than unconscious 340 dynamics ( Figure 5B, blue plots). Crucially, however, this pattern reversed following 341 replacement of the healthy connectome with the DOC connectome ( Figure 5B, red plots): 342 dynamics generated using the DOC connectome were more similar to the empirical brain 343 dynamics of DOC patients, than to conscious controls (Table 2). This observation supports our 344 hypothesis, demonstrating that unconscious brain dynamics are more compatible with the DOC 345 connectome than conscious dynamics; below, we also demonstrate that this result is not 346 specific to the chronic unconsciousness that characterises disorders of consciousness, but rather  Remarkably, these results could be replicated by replacing the original healthy structural 363 connectome with a randomised version having the same average connectivity 95 . After 364 perturbation, the model's fit to DOC brain dynamics became better than the model's fit to 365 conscious dynamics ( Figure 6A and Table 2), suggesting that unconscious dynamics are more 366 compatible than conscious dynamics with a randomised connectome. In contrast, the opposite 367 effect was observed when the original connectome was rewired into a regular (lattice) network, 368 which had a larger impact on the model's fit to DOC than control brain dynamics ( Figure 6B 369 and Table 2). Crucially, the fact that lattice networks further degrade the fit of DOC dynamics

384 385
Generalisation across datasets 386 Having identified the role of GABA-mediated inhibition for propofol anaesthesia, we next 387 sought to determine to what extent inhibition can also explain the dynamics of unconsciousness 388 arising from severe brain injury. Our rationale was that, even though these patients have not 389 been exposed to GABA-ergic agents but rather owe their condition to severe brain injury, propofol anaesthesia, but more broadly as a general mechanism responsible for the 394 characteristic dynamics of unconscious states -whether due to anaesthesia or brain injury.

396
Therefore, we followed the same "virtual anaesthesia" procedure with empirical data from 397 healthy controls and DOC patients. Intriguingly, we observed analogous results: local 398 modulation of inhibitory gain based on GABA-A receptor density allowed the model to 399 substantially improve its fit to DOC patients' brain dynamics ( Figure 7A and Table 3).

400
However, in contrast with propofol anaesthesia, the improvements were also observed when 401 the regional receptor map was scrambled, or replaced by a uniform map ( Figure 7B,C and 402 Table 3). Thus, whereas propofol anaesthesia depends on the specific distribution of GABA-A   should also lead to a model that is better able to simulate the dynamics of an anaesthetised 431 brain, than an awake brain -thereby recapitulating what we previously observed with DOC 432 brain dynamics. Remarkably, results show that -as previously observed with DOC patients -433 simulated dynamics generated from a model using the DOC connectome are more compatible 434 (lower KS distance) with the dynamics of propofol anaesthesia than with the dynamics of 435 awake subjects' brains ( Figure 8A and Table 4).

437
Furthermore, the importance of topology of the perturbed connectome was also observed for 438 propofol dynamics, with randomisation of the connectome having a smaller effect on the 439 model's ability to fit propofol dynamics than awake dynamics ( Figure 8B and Table 4).

440
Conversely, the decrease in the model's ability to fit empirical data was exacerbated for the 441 propofol condition, when the original connectome was replaced with a lattice network ( Figure   442 8C and Table 4).

444
These findings generalise our DOC results to propofol anaesthesia, indicating that the DOC 445 connectome is not only more compatible with DOC patients' brain dynamics than with healthy 446 subjects' dynamics. Rather, the generalisation to propofol anaesthesia suggests that the DOC 447 connectome may be more compatible with unconscious dynamics in general: whether arising 448 from brain injury or pharmacological intervention.  chronic lesions to brain anatomy and connectivity -can give rise to loss of consciousness and 474 its characteristic brain dynamics 14 .

475
To this end, we employed a whole-brain Dynamic Mean Field model that simulates macroscale suggests that GABA receptor density at these regions may be an especially promising candidate 500 for predicting anaesthetic effects.

502
Remarkably, our PET-informed results showed that considering GABA-mediated scaling of 503 regional inhibitory gain also improved the model's ability to simulate the characteristic 504 dynamics of DOC patients' brains, even though these patients owe their chronic condition to 505 severe brain injury rather than pharmacological intervention. This observation suggests a 506 change of excitatory-inhibitory balance in favour of inhibition, not just in the generation of 507 dynamics pertaining to propofol anaesthesia, but more broadly as a general neurobiological 508 mechanism for the dynamics that characterise unconsciousness -whether due to anaesthesia or 509 brain injury. Indeed, there is evidence that physiologically awake but unconscious DOC  world" organisation that characterises many real-world networks, including the human brain Thus, thanks to connectome replacement we can infer that the increased neuronal inhibition 561 that characterises both disorders of consciousness and anaesthesia, is functionally equivalent 562 to randomisation of the connectome. However, propofol's anaesthetic effects are mediated by 563 GABA-A receptors according to their specific regional distribution, whereas disorders of 564 consciousness can be explained in terms of a more generic increase in global inhibition -565 possibly arising from randomisation of the connectome due to anatomical lesions, whose extent 566 and location are do not follow uniform patterns. Anaesthesia may be expected to operate 567 similarly across individuals, in terms of which regions are more or less affected by propofol.

568
In contrast, each DOC patient is unique in the cause, extent and location of their brain injury.

569
As a result, whereas anaesthesia may depend on specific localised patterns, it stands to reason influencing the nodes' activity (through inhibitory modulation) or the connectivity between 581 them (through connectome randomisation). As befits such a complex dynamical system as the 582 human brain, it is likely that other paths to unconsciousness will also exist, explaining 583 phenomena such as regular sleep-wake alternation, epileptic seizures, and the effects of non-

584
GABAergic anaesthetics such as ketaminesome of which have already started to be explored 585 using whole-brain computational modelling 68,76,78,79,81 . Extending the present framework to 586 account for additional ways of losing consciousness will be a crucial endeavour, informing 587 both our understanding of brain network function and of human consciousness. Likewise, it is 588 vital to combine multimodal neuroimaging and whole-brain modelling to identify paths from 589 unconsciousness back to consciousness, using our understanding of post-anaesthetic recovery 590 to restore consciousness in DOC patients, whether by means of custom-designed drugs or deep 591 brain stimulation 57,58,86 .

593
Limitations 594 A well-known adage asserts that "All models are wrong", and the present work is no exception.   As previously reported 13 : "In the scanner, subjects were instructed to relax with closed eyes, 740 without falling asleep; 8 minutes of fMRI scan without any task ("resting-state") were acquired 741 for each participant. Additionally, a separate 5-minute long scan was also acquired while a 742 plot-driven story was presented through headphones to participants, who were instructed to 743 listen while keeping their eyes closed" 13 . The present analysis focuses on the resting-state data 744 only; the story scan data have been published separately 84 and will not be discussed further 745 here.

756
The DOC patient functional data employed in this study have been published before 13,41,115 .

757
For clarity and consistency of reporting, where applicable we use the same wording as our 758 previous study 13 .  (Table 5).

848
For the DOC patients, due to the presence of deformations caused by brain injury, rather than 849 relying on automated pipelines, each patient's brain was individually preprocessed using 850 SPM12, with visual inspections after each step. To further reduce potential movement artifacts, 851 data underwent despiking with a hyperbolic tangent squashing function. Since the controls had 852 a shorter scan duration than DOC patients, to ensure comparability between the two cohorts 853 the DOC functional scans were truncated to be of the same length as the control subjects' 854 functional scans (after removal of the initial scans).

855
To reduce noise due to cardiac and motion artifacts, we applied the anatomical CompCor    Whole-brain computational modelling 922 Whole-brain spontaneous brain activity (as quantified using blood oxygen level dependent 923 (BOLD) signal data from functional MRI) was simulated using a neurobiologically realistic  Additional factors that can influence the long-range excitatory-to-excitatory coupling between 939 brain regions, such as neurotransmission but also synaptic plasticity mechanisms, are 940 accounted for by a global coupling parameter, G. Since conductivity of the white matter fibers 941 is assumed to be constant across the brain, G constitutes the only free parameter in the model.

942
The following differential equations therefore govern the model's behaviour: Gaussian noise injected to region n". The model's fixed parameters are reported in Table 6 959 37,86,88 . Additionally, is the neuromodulatory scaling factor modulating the transfer 960 function for each cortical region in the model as a function of , the regional density of  A Balloon-Windkessel (BW) hemodynamic model 92 was then used to turn simulated regional 973 neuronal activity into simulated regional BOLD signal. The Balloon-Windkessel model 974 considers the BOLD signal as a nonlinear function of the normalized total deoxyhemoglobin 975 voxel content, normalized venous volume, resting net oxygen extraction fraction by the 976 capillary bed, and resting blood volume fraction. The BOLD-signal estimation for each brain 977 area is computed from the level of neuronal activity in that particular area. Finally, simulated 978 regional BOLD signal was bandpass filtered in the same range as the empirical data (0.008-979 0.09 Hz). The code used to run all the simulations in this study was written in optimised C++ using the 983 high-performance library Eigen. The C++ core of the code, together with Python and

986
To simulate BOLD data, FastDMF splits the problem in two steps: integrating the coupled 987 differential equations underlying the DMF model, to obtain excitatory firing rates in each brain 988 region; and using these firing rates to integrate the (uncoupled) differential equations of the 989 BW hemodynamic model and obtain BOLD timeseries.