Abstract
The rise of mathematical developments in the theories of consciousness has led to new measures to detect consciousness in a system. The Integrated Information Theory (IIT) is one of the best mathematical rooted attempts to quantify the level of consciousness in a system with Φ as the effective information generated in a system above its parts. Recently, the IIT has inspired the Perturbational Complexity Index (PCI) to detect conscious states in patients with disorders of consciousness, and it has shown to have almost perfect classification accuracy. In this study, we explore the statistical correspondence between the theoretical Φ and the experimental PCI through a neurocomputational model of coupled oscillators that can be artificially perturbed, which mainly focuses on the dynamics of collective synchronization between subsets of brain areas. Our results reveal that both measures are statistically related but, in principle, this relationship is far to be perfect. These results are discussed in the context of the model of coupled oscillators, which mainly focuses on the dynamics of collective synchronization between subsets of brain areas.
1. Introduction
One of the most challenging and still in progress tasks in science is to objectively quantify to what extent a person is conscious. One possible reason is that consciousness itself cannot be operatively defined in an easy way since it is a subjective phenomenon, and hence, it cannot be directly observed [1]. There are different theoretical approaches to consciousness [2]. For example, from a philosophical perspective, it has been proposed that there is a high order thought associated with consciousness. A conscious thought would be composed of simple percepts associated with other thoughts that provide further semantic value. Within neuroscience fields, theories can be divided into biological and functionalistic theories. Biological theories state that consciousness is a biological state of the brain; consciousness is studied through the association of different cognitive states with particular brain regions. On the other hand, from a functionalistic perspective, the existence of consciousness only requires an abstract structure to exist. For example, in silico structures could support conscious experiences as long as they obey the necessary conditions provided by the theory. A good example of functionalism is the Global Workspace account, which considers that conscious experience is the result of competition between functional networks in the brain [3] being conscious the winning network. However, within functionalistic accounts, the Integrated Information Theory (IIT) is probably the most solid and mathematically rooted attempt to define what consciousness might be and how it can be quantified [4–6]. Specifically, Tononi, in an early version of his theory [4], proposed that consciousness arises as integrated information in a system and provided Φ as a computable measure of how conscious any physical system might be. The main concept behind the calculation of Φ is to measure to what extent a system as a whole cannot be explained as a sum of its parts. As we will explain in the next section, in order to compute Φ, it is necessary to compare the information generated by the entire system with the information of the system considered as two subsystems (bipartition).
The main limitation in Φ computation is the necessary condition to find and calculate the information contained in all possible bipartitions of the system. This is a problem of complexity class NP, which cannot be exhaustively computed for systems with a large number of elements. Given that in the human brain the number of neurons is in the order of 1011 [7], even if we have the entire information from every single neuron, it would not be possible to obtain a value for Φ. The number of possible bipartitions for a brain would be in the order of a Stirling number of second kind that can be computed with and would give a number of combinations in the order of , which interestingly is much higher than the classic Eddington’s number that estimates the amount of protons in the known universe (≈1080) [8]. To solve this and other practical problems, a number of versions and estimators have been developed in different fields of research. One remarkable example was provided by Barrett and Seth [9] where they proposed a version of Φ for time series data. This version noted here as (TS stands for time series) can be applied to time series from a generative dynamical model, and it perfectly agrees with the IIT in its early version (and in the main concepts of the further versions).
In the work we present here, we will focus on this measure because it is theoretically well funded and easy to apply to time series data. However, shows the same limitation as the original Φ, which, as stated before, includes the computation of information from all possible bipartitions of the system, and hence, it is only possible to be applied to time series from systems with limited number of elements.
In a different line of research, one of the most important empirical estimators of Φ has been developed in the field of clinical assessment of unresponsive patients with disorder of consciousness. Casali et al. [10] presented a Perturbational Complexity Index (PCI) of integrated information to decide if a given patient is conscious or not. The PCI was designed to capture information and integration in the system. It quantifies the richness of information in the process of propagation of activity across the EEG channels right after discrete transcranial magnetic stimulation (TMS). The theoretical rationale of this measure is that a system with high Φ, when stimulated with TMS pulses, needs to show cortical propagation (reflecting integration of its elements) as well as diverse functional reactions at different areas of the cortex (differentiation between the elements ≈ information). Given that Φ is related to information and integration in a system, PCI is proposed as an estimator of Φ for data collected from real patients.
Although PCI has been demonstrated to classify patients into conscious or unconscious with an almost perfect accuracy rate [10], it is an indirect measure of Φ and it is not clear whether it reflects the theoretical Φ calculated from the brain as a physical system. In this study, we explore the relationship between the experimental PCI and the theoretical Φ. To investigate this question, it would be necessary to obtain both measures from the same system, and this is not easy since the latter () is designed for simple dynamical systems and the former (PCI) for real brains. The approach we take here is to obtain from an accepted neurocomputational model of whole brain resting-state activity, that is, a variant of the Kuramoto of coupled oscillators [10–13], which can be artificially perturbed to simulate TMS pulses in order to obtain a PCI for the same model.
Henceforth, in this study, we will focus on two versions of Φ: for simple dynamical systems and time series and PCI estimator developed to measure the level of consciousness in patients with different disorders of consciousness. With this work we want to explore the possible relationship between these two measures, , theoretically well founded, and the PCI, with indisputable clinical results. To the best of our knowledge, this is the first work that directly addresses this problem.
2. Methods
To investigate the potential relationship between the experimental estimator of consciousness PCI and the more theoretical index calculated over the same system of Kuramoto oscillators, we followed the next steps.
We first designed a Kuramoto model to simulate resting dynamics of the cortex. Second, we calculated several realizations of the model and obtained from the model. Third, we perturbed the system to simulate TMS pulses and compute the PCI of the model. In the next three sections, we explain these steps in detail.
2.1. The Kuramoto Model to Simulate Resting-State Dynamics
A Kuramoto model can be defined as a set of coupled oscillators modelled as the evolution of its phases according to the following set of coupled delay differential equations:where is the phase of the th oscillator on its limit cycle and is its natural frequency in radians. The control parameter is the global excitatory coupling strength, a parameter that scales all coupling strengths. is the total number of oscillators. Importantly, is the connectivity strength between each pair of oscillators and represents the time delays between these oscillators. Both adjacency matrices (connectivity strength and time delays) were obtained by Hagmann et al. [11] that defined the structure of a network coupled together according to human white matter tractography (in the work of Hagmann et al. following diffusion spectrum and MRI acquisitions, the segmented grey matter was partitioned into 66 anatomical regions according to anatomical landmarks. White matter tractography was used to determine which regions pairs were connected by putative white matter fiber tracts and to estimate their density and corresponding length, from which the structural connectivity and delays were obtained). in the network of oscillators we use here (see Figure 1). In short, each oscillator represents a cortical region of the brain located in a three-dimensional space with different connections to all other oscillators. The term allows us to dynamically modify the structural connectivity. It is an important term in this study because it can generate perturbations in the model that aimed to simulate TMS. For a Kuramoto model, the degree of synchrony between oscillators is conveniently measured by an order parameter, , that satisfieswhere measures the phase coherence or synchrony of the oscillators population; is the symbol for the imaginary operator; and is the average phase of the collective [12, 13]. indicates how coherent oscillators are in a given time and it qualifies if the phases of the collective are tightly clustered or widely distributed. We have chosen instead of other possible order parameters, as the average phase , because it does not describe with accuracy the oscillators’ collective behavior (for a given , for example, there are many possible phase distributions of oscillators).
According to several studies, the Kuramoto model shares dynamical similarities with resting-state brain functioning when it shows high metastability [14, 15]. This concept refers to high variance of which in other words can be defined as the tendency of a system of oscillators to continuously migrate between a variety of transient synchronous states, allowing a dynamical organization between the elements of the network. The system continuously goes from ordered to disordered states [16]. Then, the values of the parameters in the model were selected so that the global dynamics showed high metastability. For Hz (gamma rhythm), metastability was evident with = 3 ms, and 0.5 < < 6.5, with as the mean value of . Note that any change in can be considered a change in the mean velocity of the conduction delays between oscillators, and other values of with a different range of produced similar behaviors. For example, we found approximately the same effects for 2.5 < < 5.
The Kuramoto model was simulated for a wide range of values. As indicated before, represents the strength in the global connectivity of the model, and from a biological point of view it could be seen as a parameter to characterize integration between oscillators. Each simulation consisted of a baseline of 65 × 103 ms. As in Cabral et al. [17], we used an Euler scheme in which the time step of numerical integrations was set to .1 ms.
It would be important to state that since no experimental data are provided here, our results are obtained for parameters of the model (, , and metastability) that have been shown, in previous works [18–20], to parallel key characteristics of brain functioning. Hence, excluding the parameter , we did not manipulate the parameters in the coupled equations of the Kuramoto system.
2.2. Integrated Information in the Kuramoto Model
Integrated information was measured with the version and it is mainly based on the concept of effective information () [9].
Let be a multivariate random variable that takes values in the space . It is evident that the dimensions of are the number of elements in the system that generates . The effective information generated by a system in its current state about the state with respect to a bipartition of it is defined by the mutual information generated by the entire system minus the sum of the mutual information of its parts in the bipartition:
Mutual information in bits can be calculated with the expressiona measure that gives the average bits that can be predicted in given the state [21]. The calculation of mutual information includes the calculation of probabilities and joint probabilities of any estate and . Integrated information is the effective information with respect to the minimum information bipartition (MIB):wherewith stating for “the minimum number in the set” andis a normalization factor to correct excessive unbalanced bipartitions. here stands for Shannon entropy.
To apply to the Kuramoto model described in the previous section, we divided the 66 regions of the original network into 6 clusters proposed by Hagmann et al. in the original work [11]. This simplification allowed exploring all possible bipartitions in the system (the Stirling number for this case gives = 31 possible bipartitions). In addition, and following [22], time series for each cluster were calculated as the synchrony between the oscillators belonging to that cluster (). Then, we characterized these series as synchronized or not synchronized by constructing new binary time series from each . In our study, we considered a synchronization threshold of = 0.8 and 1 was assigned for all > γ. We selected this value because it was the median of , and using the median for thresholding eliminates the possible influence of extreme values, due to its robust properties. In addition, there was a theoretical reason for this election. The value of γ = 0.8; it was the one used in [22], which would allow us to compare the results we obtained.
Finally, was calculated to these binary series and the result was taken as the complexity value of the Kuramoto model. We set = 150 ms. Values of < 150 ms gave negative estimations of Φ and for > 150 ms the pattern of results for different values did not change.
2.3. Estimation of PCI in the Kuramoto Model
In order to calculate an estimation of the PCI in the Kuramoto model, it was necessary to solve two problems in the simulation process. The first one was to perturb or stimulate the system from an external source (to emulate the effects of TMS), and the second difficulty we found was to calculate reliable ERPs for the final PCI calculation. The problem of the stimulation was easily solved since it has been done in other studies. For example, Hellyer et al. [23] simulated external stimulation to a similar Kuramoto model we present here by increasing the connectivity between some nodes of the network. Similarly, Ibáñez-Molina and Iglesias-Parro [24] stimulated another Kuramoto version by transiently increasing the connectivity of key oscillators in the network. Hence, in our study we followed these studies and perturbed the system by transient increases in the connectivity between six oscillators located in the parietal cortex (see Figure 1). We arranged this procedure by introducing for these oscillators during short periods of time of 5 ms. We randomly repeated this stimulation 15 times for each numerical integration of the model.
The next step was to build reliable ERPs with the resulting phases from the oscillators. To achieve this goal we simulated EEG series, and then, ERPs were calculated for each model by averaging the segments associated with each period of stimulation. We explain this procedure in detail in the next two sections.
2.3.1. EEG and ERPs Simulation
The EEG activity from 32 sensors was simulated for each condition in agreement with the following weighted sum of the activity in each oscillator:where is the time series from sensor th and is the weighted contribution of source th in sensor th. Each was calculated using a standard forward model algorithm [25] applied according to the Talairach coordinates of the oscillators. After that, each oscillator was considered a cortical source. Second, the weights of these sources were normalized to a maximum value of 1. EEG signals obtained with this procedure gave a set of signals that changed in amplitude and frequency variations around 60 Hz (natural frequency of the oscillators). Because in the ERPs the interesting information is in the amplitude, we calculated the envelope of with the Hilbert transform. Envelopes of the signal were then used to construct the ERPs with the average of all segments in each realization of the model. Formally, ERPs were built with the analytical signal (in the complex plane) of which iswhere is the original signal and represents the real part of the new complex series and is the imaginary part from the Hilbert transform with as the imaginary operator. In the right side of the expression, and stand for the angle and modulus representing the complex values in Euler’s notation. The modulus is the amplitude or the analytic power of the signal and can be easily calculated with
These new series were considered the activity from each sensor and the ERPs were built with segments extracted from them. For each sensor ,where is the number of segments in a single realization of the model and is a vector with values from −600 to 600 ms after stimulation. = 15 for all numerical integrations (see Figure 2 for a visual inspection of ERPs).
2.3.2. PCI Calculation
ERPs were the input signal for PCI calculation. PCI was obtained following the original algorithm in Casali et al. [10]. For each condition of the Kuramoto model, first, we built a binary source matrix with dimensions corresponding to sensors and time steps from 0 to 1200 ms after stimulus presentation (SS). The signals were downsampled ten times to obtain a sampling rate similar to real data. SS = 1 if the absolute value of the poststimulus simulated signal was higher than the absolute value of the maximum entry encountered in any sensor and any time step from the prestimulus baseline, and SS = 0 otherwise. SS was used as input for the Lempel-Ziv measure [26] to estimate the algorithmic complexity (). gives the number of chains with nonredundant information contained in SS. The algorithm seeks for the minimal number of patterns necessary to describe the sequence. For random sequences, the asymptotic behavior of the measure is log2, where is the binary entropy for length where is the probability to find a “1” in the binary sequence of length . PCI is defined as the normalized value of :
3. Results and Discussion
The results for the Kuramoto simulations are characterized in the first place to understand the basic dynamics of the model. In addition, we include graphical descriptions of the ERPs to visualize the structure of the averaged waves at each sensor from simulated perturbations. Finally, we describe and compare Φ and PCI taking into account the values of with high and low metastability, which as mentioned above is a necessary condition for brain dynamics at resting state.
3.1. Kuramoto Simulations in the Baseline Condition
In Figure 3(a) we show a diagram with the behavior in the baseline condition at several values of its coupling parameter . The most important property in the evolution of is the metastability that can be estimated by the variability of . As can be seen, there is a bifurcation for that indicates the end of high metastability and hence the dynamics of the model for > 6.5 should be taken with caution since in principle, there is no functional correspondence with real cortical dynamics.
(a)
(b)
In addition, it is important to note that the frequency structure in is not fixed for all s. One can observe in Figure 3(a) that the frequency of seems to increase with the increase of . To better understand this phenomenon, we include a spectral decomposition of the evolution of in Figure 3(b). Surprisingly, we found a complex landscape in the oscillatory structure of . The general structure of the spectral diagram showed a resemblance with the bifurcation diagram of the classical logistic map. This similarity appeared because it showed a bifurcation-like proliferation of frequency components as the parameter increased. The nature of this spectral structure, however, was not explored and goes beyond the goals of this study. If we inspect Figure 3(b), it can be stated that the end of slow oscillatory properties of high metastability was evident at . Above this value, the critical coupling strength of certain clusters is achieved, so their synchronization becomes stable, while the order parameter of the whole system remains <1. As increases further, larger and larger synchronized clusters are formed, resulting in a reduced number of components, until approaches 1, with ultimately only a single component as tens to infinity [27]. Accordingly, in Figure 3(b), the main frequency of the signal slowly increased with , and the components of seemed to increase with as well following a complicated pattern. It is also noteworthy that in Figure 3(b) another bifurcation-like region for can be perceived that consists in a reduction in the number of components. Hence, by the end of the landscape, seems to be more simpler with less oscillatory properties and probably this could lead to low values of and PCI.
The shape of this diagrams led us to consider that and PCI could be sensitive to the bifurcation region in . If metastability is a necessary condition for brain functioning [16] it would be reasonable to think that should diminish in the low metastability region. The same should be true for PCI if this measure is closely related to . As we will show in the next sections, the reduction after metastability was only found for .
3.2. Comparisons between and PCI
The relation between and PCI was assessed using Pearson product-moment coefficient between the values of obtained for each of the -levels (from .5 to 15) and the corresponding values of the PCI obtained for those same levels of . The results showed a nonsignificant negative correlation ( = −.21, ). Thus, apparently, and PCI are linearly independent. However, taking into account metastability values, a bifurcation for = 6.5 is apparent. We therefore divided the series into two parts according to the bifurcation, before (from = .5 to 6) and after (from = 6.5 to 15) the bifurcation, and recalculated the correlation between and PCI in each of those two parts. In this case, results showed a significant positive relation between and PCI before ( = .64, ) and after (r = .51, ) the bifurcation.
A graphical description of the evolution of and PCI and the metastability of the model can be observed in Figure 4. Values were calculated with .5 -steps. A visual inspection of Figure 4 shows that both and metastability exhibited a big decrease around ≈ 6.5 indicating the dependence between and metastability. Actually, the correlation between metastability and over the whole range of was significant ( = .68; ). However, the PCI did not show a significant decrease in this region of . In fact, the PCI evolution seems to progressively increase with . The fact that there is no decrease in PCI does not mean that this measure is not related to ; as seen before, a closer exploration of both measures indicated a positive relation between them. Moreover, one can observe in Figure 4 that the PCI trend seems to stop after ≈ 11 which is in agreement with the bifurcation shown in Figure 3(b) showing that the synchrony dynamics are simpler for this region.
4. Conclusions
Under the assumption that conscious states come from integrated information in a system, various metrics have been proposed to try to quantify consciousness. In the present work, we tested two of them using a neurocomputational model. On the one hand, the theoretically well founded Φ has been proposed as a way to quantify the total amount of information that a conscious system can integrate [28]. On the other hand, the PCI distinguishes conscious versus unconscious states at a single patient precision [10]. Under the assumption that conscious states correspond to a distributed but nonuniform spatiotemporal pattern of current sources, Casali et al. applied a standard data compression scheme (the Lempel-Ziv algorithm) to distinguish between conscious and unconscious states. Despite the excellent results at applied level, the claim that the measure is theoretically grounded in a conceptual understanding of consciousness deserves a closer look. In the present work, we have tackled the possible relationship between these two measures of the degree of consciousness in a system.
As stated previously, according to the IIT, wakeful consciousness requires the ability to integrate information across multiple brain regions with a high degree of differentiated activity. Thus, loss of consciousness may result as a consequence of a loss of integration as well as a loss of differentiation (or both). Due to our manipulations in the present paper (i.e., increase of -values), the reduction in consciousness indicated by values would reflect such stereotypical behavior across different oscillators. Increasing the mean field connectivity of the model results in an increase of synchrony, as revealed by the Kuramoto order parameter () and in a reduction of its variability (metastability). From this point of view, results obtained from metastability and converge with the theoretical predictions that suggest that metastability is a necessary condition for healthy brain functioning and consciousness [16].
In general, as increases, the system as a whole is more coherent and hence, it is more integrated. When the system is above the bifurcation point (), synchrony is very high and the dynamics of the binary time series from the clusters tend to be 1 all the time (). Hence, information will tend to be low. So what happens is that information is much lower as the system crosses the bifurcation point. One can objectively see this by observing the entropy and mutual information (in bits) of the system (see Figure 5). This descent in the entropy after the bifurcation point could be responsible for the apparent inability of the PCI to capture the dynamic of the system after the bifurcation. In this respect, [29] have suggested that the inability of Lempel-Ziv to compress efficiently low entropy sequences is due to the inability to cope with long runs of identical symbols. In this respect, the Perturbative Integration Latency Index that characterizes the latency of extinction of a massive stimulation perturbing a basal state without drawing upon Lempel-Ziv algorithm, recently proposed by [30], is a promising option that future works could explore.
An important finding in this study is the significant positive correlations between and PCI before and after the bifurcation point. These correlations might indicate that there is a modulation in the PCI when changes due to connectivity manipulations. Hence, from the exploration we carried out in this study we can claim that the PCI is sensitive to modulations, and that this is true when the system is considered in a coherent region of metastability (high versus low).
One limitation in our study is that we have not found a critical value of at which both measures reached a maximum. An outstanding correspondence between and PCI would have led to an optimum parameter that characterizes the system in terms of integrated information for both theoretical and empirical estimators. It is evident that the reasons why we did not find this perfect convergence could be that our model is an oversimplification of brain functioning, or the procedure to calculate and PCI relies on many simplifications for the characterization of the system. Due to this oversimplification of the brain dynamics, the fact that the PCI after the bifurcation point tend to increase could be due to intrinsic characteristics of the model. However, to the best of our knowledge, Casali et al. did not test PCI when patients exhibit loss of consciousness due to global synchronization in the cortex (epileptic states, for example). Future works could explore the possibility that the loss of consciousness due to global synchronization cannot be fully captured by PCI. However, it is noteworthy that when we consider the metastable region of the model ( < 6.5) the maximum values for and PCI are found in the short range 2 < < 3, and the minimum values are found for when the connectivity of the system is relatively low. These two findings might indicate that and PCI could have a better agreement for high levels of metastability, and if this is true, it is not surprising that PCI is a good clinical indicator of conscious states.
Another potential limitation in our study could be the algorithm used to estimate originally proposed by Barrett and Seth [9]. This algorithm produces negative values that could hinder the interpretation of the obtained results. These limitations have given rise to new versions of Φ in which the disadvantage of negative values is solved [31, 32]. Although both measures are highly correlated, future studies could include both estimators.
Disclosure
This study was partially reported as oral presentation in the local workshop “Modelos Globales de Dinámica Cerebral,” held in Seville on December 12, 2017.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
Acknowledgments
This work has been supported by the Spanish Government and the European Union (via ERDF funds) through the Research Project MTM2014-61312-EXP and by Junta de Andalucía (Biomedical and Heath Science Research Project PI-0386-2016).