Unconsciousness reconfigures modular brain network dynamics

The dynamic core hypothesis posits that consciousness is correlated with simultaneously integrated and differentiated assemblies of transiently synchronized brain regions. We represented time-dependent functional interactions using dynamic brain networks, and assessed the integrityof the dynamic core by means of the flexibility and largest multilayer module of these networks. As a first step, we constrained parameter selection using a newly developed benchmark for module detection in heterogeneous temporal networks. Next, we applied a multilayer modularity maximization algorithm to dynamic brain networks computed from functional magnetic resonance imaging (fMRI) data acquired during deep sleep and under propofol anesthesia. We found that unconsciousness reconfigured network flexibility and reduced the size of the largest spatiotemporal module, which we identified with the dynamic core. Our results present a first characterization of modular brain network dynamics during states of unconsciousness measured with fMRI, adding support to the dynamic core hypothesis of human consciousness.


INTRODUCTION
The dynamic nature of consciousness is already manifest in William James' 1890 metaphor of the "stream of consciousness" [1]. According to James, consciousness "does not appear to itself chopped up in bits" but instead flows and reconfigures itself over time, much like the stream of a river. This metaphor has informed several scientific hypothesis concerning the relationship between conscious experience and neural processes in the human brain, perhaps the most influential being Edelman and Tononi's dynamic core hypothesis [2][3][4]. According to this hypothesis, consciousness must be understood as a process that unfolds over time (the "dynamic core") comprising an ever-changing network of regions that exchange information over relatively short time spans. This position is in sharp contrast with the proposal that activity within certain specific and localized neuroanatomical structures comprise the neural basis for consciousness [5].
The properties of the dynamic core are shaped by the tension between two observations concerning the nature of subjective experience [6]. On the one hand, conscious experience is extraordinarily rich, in the sense that admits countless different contents and variations. On the other hand, each individual conscious experience is integrated and cannot be decomposed into more elementary independent subexperiences (e.g. the senses of smell, sight and hearing are never experienced separately, but always fused together into the ongoing contents of consciousness). Thus, the dynamic core should present a very large number of possible configurations, corresponding to the multitude of available conscious experiences. However, these configurations must be constrained to represent highly integrated brain states. While a sequence of randomly connected brain regions is a highly differentiated process, each individual state in the sequence lacks integration (Fig. 1A, first row). Conversely, a stable sequence of fully connected regions lacks differentiation (Fig. 1A, second row). In between these two examples, the dynamic core consists of a sequence exploring an ample repertoire of highly integrated brain states (Fig. 1A, third row).
The presence of simultaneously integrated and differentiated dynamics can be captured by different metrics of neural complexity [7]. For instance, unconscious brain states have been characterized by their entropy [8], algorithmic complexity [9] and long-range information integration [10]. Furthermore, the persistence and complexity of activity evoked by external stimulation depends on the level of consciousness for different brain states such as deep sleep, general anesthesia, and disorders of consciousness [11]. These changes are consistent with the identification of subcritical dynamics as a signature of unconsciousness, which results in stable collective neural behaviour that is incompatible with the rapid reconfiguration of metastable states postulated by the dynamic core hypothesis [12][13][14].
Perhaps the most straightforward method to assess the dynamic core hypothesis is also the most underexplored to date. The theory of complex networks provides a framework to directly assess the presence of integration and segregation in neuroimaging data acquired during different states of consciousness [15]. A sequence of brain states can be represented as a multilayer network, with each layer encoding transient functional interactions between brain regions during a given time period, while the dynamic core can be represented as a time-dependent module evolving in this network (Fig. 1B) [16]. Over the last years, modularity maximization algorithms have been applied to multilayer networks to reveal the rapid and transient structure of whole-brain dynamic networks. These dynamics are sensitive to factors such as cognitive performance [17][18][19][20] and neuropsychiatric impairments [21]. In spite of these advances, however, the relationship between consciousness and the modular structure of multilayer brain networks remains to be investigated.
To clarify this relationship, we constructed multilayer connectivity networks from functional magnetic resonance imaging (fMRI) recordings acquired during the different stages of first row depicts a sequence of highly differentiated but segregated states (green assemblies). The second row depicts a sequence of highly integrated states lacking differentiation (blue assemblies).
Finally, the third row depicts a highly integrated dynamic core shifting through several different configurations (red assemblies). B. These temporally evolving assemblies can be represented as modules unfolding in a multilayer network (layers S 1 , S 2 , S 3 ), with inter-layer couplings for node k given by ω kS 1 S 2 , ω kS 1 S 3 , ω kS 2 S 3 , and resolution parameters γ S 1 , γ S 2 ,γ S3 . In this example, all interlayer couplings and resolution parameters are equal, and interlayer connections only appear between consecutive slices (e.g. ω KS1S3 = 0). Figure adapted from Cavanna et al, 2018 [22] human non-rapid eye movement (NREM) sleep, and under the effect of propofol, a general anesthetic which increases inhibitory neurotransmission. Our main purpose was to obtain the time-dependent modular structure of these networks using the multilayer Louvain algorithm, a method with multiple free parameters related to the connectivity strength between temporal layers, and the characteristic size of the detected modules (Fig. 1B) [23]. Previous uses of this algorithm either employed an ad-hoc choice of parameter values, or performed an exhaustive exploration of parameter space [16][17][18][19][20][21]24]. We avoided these suboptimal solutions by introducing a new benchmark for the detection of modules in time-dependent networks with scale-free degree and module size distributions, adapted from a benchmark developed for static networks [25]. After parameter selection, we applied the multilayer Louvain algorithm to obtain the time-dependent modular structure of fMRI functional connectivity networks. We focused our analysis on two metrics related to the dynamics of modules over time: the flexibility (number of times a node switched its module assignment over time) and the size of the largest multilayer module, which we identified with the dynamic core.

MATERIALS AND METHODS fMRI data sets
We analyzed two different fMRI data sets comprising healthy subjects, both during conscious wakefulness and during states of reduced consciousness.   [27]. Pulse oxymetry and respiration were recorded via sensors from the Trio [sampling rate 50 Hz]) and MR scanner compatible devices (BrainAmp MR+, BrainAmpExG; Brain Products, Gilching, Germany), facilitating sleep scoring during fMRI acquisition. Further details on data acquisition and preprocessing can be found in a previous publication [28].

Module detection in multilayer networks
We consider a multilayer network with adjacency matrix given by A ijs , were i, j index the network node and s indexes the layer, which is here interpreted as a temporal dimension.
Given a certain partition, its multilayer modularity (Q) is computed as, where k js = i A ijs , µ = 1 2 jr (k jr + s ω jrs ), m s = j k js and δ(g is , g jr ) equals 1 if node i of layer s belongs to the same module as node j of layer r. γ s is the resolution parameter for layer s, and ω jrs represents the interlayer connectivity of node j between layers r and s. Here we consider the same γ s for all layers, and ω jrs = 0 only if r and s are consecutive layers, furthermore, all non-zero entries of ω jrs are equal.
Based on this target modularity function, we applied a generalized multilayer version of the Louvain algorithm to detect and track modules over time (http://netwiki.amath. unc.edu/GenLouvain/GenLouvain). The output of this algorithm was a 2D matrix G it presenting the module assignment of nodes (i, rows) across time (t, columns).

Benchmark for time-dependent module detection
To approach the problem of constructing a benchmark for module detection in complex time-dependent networks, we first reproduced a benchmark for static complex modular networks introduced by Lancichinetti et al. [25].
This benchmark assumes that both the degree and the module size distributions follow power laws, with exponents γ and β, respectively. The number of nodes is N and the average degree is < k >. A mixing parameter is defined as µ = #intra-modular links # total of links , so each node shares a fraction 1 − µ of its links with nodes of its own module and a fraction µ with the other nodes of the network. Both the degrees of nodes, k, and the size of the modules, s, are constrained between predefined minimum and maximum values. The values of k min and k max are chosen such that the average degree is < k > and s min and s max are chosen so that any node can belong to at least one module, s min > k min (1 − µ) and s max > k max (1 − µ).
The proposed temporal evolution for our benchmark consisted of two different dynamic processes adapted from Granell et al., 2015 [29]: dynamics of merge-split and grow-shrink of modules. The combination of these two processes allowed us to represent the most frequent behaviours seen in the dynamics of real modular systems.

Null model based on division of communities: merge-split
We present an algorithm to induce this type of dynamics in heterogeneous networks of N nodes constructed according to the algorithm presented by Lancichinetti et al. [25]. We start from a network with N nodes and apply the following steps: 1. A module named C, sufficiently large for the division into smaller sub-modules, is chosen at random.
2. The nodes belonging to C are assigned to two sub-modules, C1 and C2. A fraction x of the nodes in C belongs to C1, and (1 − x) belongs to C2. 3. For each node i in C2, we apply the following steps: (a) µ 21i of node i is calculated as: µ 21i = #Links with C1 nodes #Total links (b) A mixing parameter between communities is chosen per node, µ 12 .
(c) While µ 21i > µ 12 nodes a, b and c are searched such that they meet the conditions set forth in the division step presented in Fig 3A (division) : node a and node b belong to C1 and are connected to node i, node c belongs to C2 and is not connected to node i.
(d) The link between node i and node b is deleted and a new link is created connecting it to node a. This is repeated until there are no more nodes a, b and c meeting these conditions or until µ 21i < µ 12 . The rewiring scheme is presented in Fig. 3A.
(e) The final adjacency sub-matrix is saved as the dynamic network at time t = i.

Null model based on contraction of communities: grow-shrink
We also present an algorithm that selects two random modules from a certain complex network with N nodes, and then grows the size of one module at the expense of the other.
The steps are as follows: 1. Two modules C1 and C2 are chosen at random.
2. The adjacency matrix is reordered so that the two modules appear consecutively. The adjacency submatrix containing modules C1 and C2 is selected.
3. C1 nodes are re-tagged and rewired, changing their membership to C2: Let n 1 be the number of C1 nodes, ranging from 1 to the X − th node of C1, where X is the integer part of n 1 x. The following steps are applied to each node i of this set: (a) µ 12i of node i is calculated as, µ 12i = #Links with C2 nodes #Total links (b) A mixing parameter between modules is chosen per node, µ 12 (c) While µ 12i > µ 12 , nodes a, b, c and d are searched such that they meet the conditions set forth in Fig 3A (contraction).
(d) The link from node i to node b is deleted, and a new link is created between node i and node a. This is repeated until there are no more nodes a and b fulfilling these conditions, or until µ 12i < µ 12 . The rewiring scheme is presented in Fig.   3A.
(e) Node i is removed from C1 and added to C2.
(f) The average degree of intercommunity links of C2 nodes (< k iC2 >) and the degree of node i (k i ) are calculated.
(g) While k i < < k iC2 > pairs of nodes whose C2 intramodular degrees are between < k iC2 > and k max are selected. The links between those nodes are deleted, and new links between node i and other nodes in C2 with k < < k iC2 > are added.
(h) The final adjacency submatrix is saved as the dynamic network at time t = i.

Construction of dynamic networks from fMRI data
After preprocessing, we extracted the average BOLD signals from 401 regions of interest, for all conditions and participants. These regions correspond to a sub-parcellation of the automated anatomic labeling (AAL) atlas [30], resulting in more detailed anatomical subdivisions [31]. Next, we applied a sliding window procedure to compute Pearson's linear correlation coefficient between all pairs of ROIs during 60 s windows, moving each window for one volume to compute the next correlation value. This resulted in a weighted correlation matrix C ijs containing in its i, j, s the correlation coefficient between ROIs i and j during a 60 s window starting at volume s. Since the multilayer Louvain algorithm requires positive link weights, we opted to threshold the correlation matrix to produce a binary adjacency matrix A ijs which equals 1 only if correlation between nodes i, j at volume s is larger than a certain threshold, and zero otherwise. Since we used our dynamic benchmark to estimate to optimal γ and ω values, we chose a threshold to result in the same average degree used for the benchmark (< k >= 20), so that both networks were comparable.

Node flexibility and the largest multilayer module
We defined two metrics based on the module membership matrix G it given by the multilayer Louvain algorithm. First, we defined the node flexibility F i as the normalized number of module membership changes for that node: where T is the total number of volumes in the recording.
We defined the largest multilayer module size LM M as the normalized size of the largest module in G it , i.e.
where N is the number of nodes and T the total number of volumes in the recording, so N T is the maximum possible size for the giant module.

Time-dependent benchmark and parameter selection
We investigated the performance of the multilayer Louvain algorithm based on introducing equally weighted (ω) connections between consecutive temporal layers (see Fig. 1B).
Also, we considered equal resolution parameters across layers (γ). Thus, the module detection algorithm depended only on these two parameters.
We applied our benchmark to determine the optimal parameter values. For this purpose, we introduced a grid of γ and ω values, and for each pair of values we measured the Rand index between the ground-truth modules and those detected using the multilayer Louvain algorithm, averaging the results over 500 independent realizations (see Fig. 3B). The ground truth used to determine the performance of the multilayer Louvain algorithm was determined by applying the standard Louvain algorithm to each time point and then selecting the partition maximizing the overlap with the heuristic expectation. It is important to note that this could only be done for networks with a priori information about the allocation of nodes in temporal modules, such as in the merge-split and grow-shrink dynamics.
The optimal parameters obtained following this procedure were γ = 0.55 and ω = 1 for both benchmarks (values are indicated as black boxes in Fig. 3B).
Finally, while the merge-split and grow-shrink dynamics preserved the amount of total links in the network, the degree distributions could be changed during temporal evolution.
To quantify this possibility we measured α(t), i.e. the time-dependent exponent of the power law degree distribution, and then computed its standard deviation across time for both dynamics (see Fig. 3D). We observe that α(t) changed for the contraction dynamics, while it remained stable for the division dynamics.

Modular structure of dynamic brain connectivity networks
We applied the multilayer Louvain algorithm using the optimal parameters inferred from the benchmarks to dynamic functional connectivity networks obtained from fMRI data.
First, we investigated healthy subjects during conscious wakefulness, and during the wakesleep cycle, consisting of four stages of increasing depth and progressive loss of conscious content (N1, N2 and N3, respectively).  thalamic and anterior cingulate nodes, consistent with the observation that the thalamus becomes deactivated and disconnected from sensory cortices during early sleep [28]. Conversely, flexibility increased during sleep in parietal, frontal and prefrontal regions associated with higher cognitive functions, with the strongest increases seen during N3 sleep.
Finally, we compared the regional probability of belonging to the largest multilayer module in wakefulness vs. deep sleep and propofol-induced loss of consciousness (LOC). Figure   6A presents a comparison of significant changes in this probability between wakefulness, N3 and LOC. While changes were more widespread and significant during N3 sleep, LOC was also associated with decreases in sensorimotor regions, and increases in frontal regions. A scatter plot of the change in the probability of belonging to the largest multilayer module for LOC vs. N3 shows that even though less regions were significant for LOC, the pattern of changes was similar to that measured during N3 sleep (Fig. 6B). Also, both N3 and LOC were characterized by smaller sizes of their largest multilayer modules. As shown in Fig.   6C, the normalized size of the largest modules was significantly smaller for N3 and LOC compared to wakefulness. strated the usefulness of adopting a dynamic approach towards brain activity [33], shedding light on the large-scale brain dynamics of physiological and pathological brain states [21], the coordination between distributed networks of regions that gives rise to human cognition and behavior [17][18][19][20], and the relationship between multimodal brain activity [34], anatomy [35], and genetics [36]. Global states of consciousness such sleep, general anesthesia, and disorders of consciousness have also been explored from this perspective, leading to the general observation of a reduced and more stable repertoire of brain state associated with unconsciousness [22,[37][38][39][40][41].
These results depend on the development of new methods capable of extracting temporally recurrent brain activity patterns in data sets comprising hundreds of thousands of voxels. At the core of dynamic analyses of fMRI data is the need for data-driven identification of recurrent brain states, and their mutual transition probabilities. This can be achieved by different means, such as clustering a point process representation of the data obtained by thresholding the BOLD signals [42][43][44], or applying more sophisticated methods based on hidden Markov chains [45]. One of the first and most frequently applied methods to unveil the dynamics of whole-brain fMRI activity consists of adopting a multi-layer network representation, where each layer corresponds to a different time window and intralayer connections indicate BOLD signal correlation (i.e. functional connectivity) [16]. By introducing interlayer connectivity, modularity optimization algorithms can then be applied to this representation for the data-driven discovery of temporal modules, i.e. non-overlapping groups of anatomical regions presenting high and low intra-and inter-group functional connectivity, respectively [16,46].
Module detection in temporal networks remains one of the most used methods to map the dynamic evolution of whole-brain activity, yet presents shortcomings related to the need to adjust free parameters. Our contribution advances our understanding of whole-brain dynamics from the two aforementioned perspectives. First, we developed a systematic framework to evaluate the performance of modularity optimization algorithms by constructing a dynamic benchmark with heterogeneous and realistic node distributions, allowing the justified choice of model parameters. Next, we started from these parameters to show that sleep progressively stabilizes module assignments in sensory and associative cortical regions, with the opposite behavior observed for frontal, parietal and temporal regions spanning key default mode network hubs. Sleep and propofol anesthesia reduced the regional probability of belonging to the largest multilayer module for most occipital, temporal and parietal brain regions, with these changes being correlated between two unconscious conditions, suggesting an overlapping mechanism implicated in loss of consciousness.
The size of the largest multilayer module represents a new metric to quantify integration simultaneously in the spatial and temporal domains. Since modules are computed from the multilayer representation of the network, this metric is different from the largest module at any given time point. The largest module can shrink to a few nodes or regions at a certain time, but grow at later times to encompass most of the brain; in other words, the identity of the module over time is taken into account. This is in line with the definition of the dynamic core by Edelman and Tononi, which corresponds to the same temporal process which involves different anatomical regions at different times, but nevertheless maintains its identity through time [2][3][4]. Accordingly, consciousness can vanish as consequence of a static or fragmented dynamic core (see Figure 1). The reduced size of the largest multilayer module during deep sleep and general anesthesia module seemingly supports the second possibility; furthermore, the regional correlation between the changes in the likelihood of belonging to the largest module during deep sleep and anesthesia suggests that this metric could capture a signature of fading consciousness present in both conditions.
In terms of flexibility, i.e. switching rate between module assignment, we did not observe significant changes at the global level, but detected two sets of regions with opposite changes during the descent to deep sleep. This dissociation became especially evident during N3 sleep, when flexibility increased in frontal and parietal regions, and decreased in primary sensory regions; these patterns persisted during N3 or slow wave sleep. We can speculate that the progressive fragmentation of largest module results in increased switching behavior in fronto-parietal regions, while the stabilization of module dynamics in sensory regions could be the results of neurophysiological processes preventing arousals from sleep, such as K-complexes and slow waves [47][48][49]. These processes seem specific to sleep, since general anesthesia did not results in local flexibility changes.
While this is the first study applying a multilayer network approach to understand the temporal dynamics of fMRI data during unconsciousness, there is one antecedent following similar methods applied to EEG data measured in patients suffering from disorders of consciousness [50]. This work demonstrated decreased flexibility in relatively slow frequency bands (mainly in the alpha range) and the opposite behavior over longer time scales. Our work might represent the counterpart of these findings in the hemodynamic activity measured using fMRI during sleep [51], where certain networks maintain high levels of network flexibility to revert towards a state of wakefulness upon arousal (as opposed to brain-injured patients who might not spontaneously revert towards conscious wakefulness).
Future research using multilayer network modularity to describe temporal fMRI dynamics should base parameter selection on an adequate benchmark. Using model networks with known behavior is a widely adopted practice in many applications of network science; ideally, this model should represent as many features of the empirical networks as possible [25].
Our methodological development filled this gap by introducing a benchmark for temporal modular networks with heterogeneous connectivity giving rise to realistic node and module size distributions. This benchmark can also find future applications in the comparison of different methods for modularity detection, and in the construction of synthetic datasets with known modular structure to serve as null models.
In summary, we investigated for the first time modular brain network dynamics during states of unconsciousness, finding converging evidence of a reconfiguration of the largest multilayer module during deep sleep and general anesthesia. We interpreted these changes in the light of the dynamic core theory, concluding that unconsciousness results in its fragmentation in spite of preserved stability. Future studies should assess whole-brain dynamics simultaneously with different methods to understand whether the dynamic core fluctuates over scales inaccessible to fMRI, and whether these fluctuations are manifest at the behavioral and cognitive levels.