Capturing the non-stationarity of whole-brain dynamics underlying human brain states

NeuroImage 244 (2021) 118551 Contents lists available at ScienceDirect NeuroImage journal homepage: www.elsevier.com/locate/neuroimage Capturing the non-stationarity of whole-brain dynamics underlying human brain states J.A. Galadía , b , ∗ , S. Silva Pereira b , Y. Sanz Perl c , d , M.L. Kringelbach e , f , g , I. Gayte a , H. Laufs h , E. Tagliazucchi c , J.A. Langa a , G. Deco b , i a Departamento de Ecuaciones Diferenciales y Análisis Numérico, Universidad de Sevilla, Spain b Center for Brain and Cognition, Universitat Pompeu Fabra, Barcelona, Spain c Buenos Aires Physics Institute and Physics Department, University of Buenos Aires, Argentina d Laboratory of Experimental Psychology and Neuroscience, Institute of Cognitive and Translational Neuroscience, INECO Foundation, Favaloro University, Argentina e Department of Psychiatry, University of Oxford, UK f Centre for Music in the Brain, Department of Clinical Medicine, Aarhus University, Denmark g Life and Health Sciences Research Institute (ICVS), School of Medicine, University of Minho, Portugal h Department of Neurology, Christian-Albrechts-University Kiel, Germany i Institució Catalana de la Recerca i Estudis Avançats (ICREA), Universitat Pompeu Fabra, Spain a r t i c l e i n f o Keywords: Brain dynamics Attractors Energy levels Sleep Consciousness Model transform a b s t r a c t Brain dynamics depicts an extremely complex energy landscape that changes over time, and its characterisation is a central unsolved problem in neuroscience. We approximate the non-stationary landscape sustained by the human brain through a novel mathematical formalism that allows us characterise the attractor structure, i.e. the stationary points and their connections. Due to its time-varying nature, the structure of the global attractor and the corresponding number of energy levels changes over time. We apply this formalism to distinguish quantitatively between the different human brain states of wakefulness and different stages of sleep, as a step towards future clinical applications. 1. Introduction High-dimensional complex systems usually exhibit a large number of attractors or stable steady states, and are thus far from thermodynamic equilibrium. The energy or attractor landscape (AL) characterizes all these states, limits the possible future behavior of the system, defines states with greater or smaller probability, and establishes the propensity of a system for more or less probable changes. In Engineering and Mathematics, the AL can be represented within the Dynamical System Theory (DST) by means of a Lyapunov function, while in Physics the so-called energy landscape is usually expressed by means of a Potential function. Progress has recently been made in connecting both points of view ( Yuan et al., 2014; Zhou et al., 2012 ). However, to date this has not been applied to brain dynamics. Here, we develop a novel mathematical framework that allows us to capture for the first time the non-stationarity character of the attractor landscape of brain dynamics in different brain states. A dynamical system is usually described by ordinary or partial differential equations -continuous time- or difference equations -discrete time-. Generally, the phase space of the system comprises the values of ∗ Corresponding author at: Departamento de Ecuaciones Diferenciales y Análisis Numérico, Universidad de Sevilla, Spain. E-mail address: [email protected] (J.A. Galadí). the dynamical variables. In DST, the global attractor is identified with a compact invariant subset of the phase space ( Hale, 1988; Robinson et al., 2001 ). The global attractor determines the global asymptotic behavior of a dynamical system , i.e. it describes the past and the future of the system. The structure of the global attractor can be given in terms of isolated invariant subsets (typically, stationary points or periodic orbits) and the trajectories in phase space connecting them ( Bortolan et al., 2020; Carvalho et al., 2012 ). According to DST, if the dynamical system supports a Lyapunov function, all isolated invariant subsets of the global attractor can be ordered by their level of attraction or stability. Then, a sequence of energy levels can then be defined ( Aragao-Costa et al., 2012 ), and, for instance, dynamical systems with billions of stationary points in which detailed descriptions of the global attractor and the AL would not be available, can be characterized by just a few tens of energy levels. Indeed, the number of energy levels (NoEL) allows us to characterize the attractor landscape without having to calculate the explicit expression of a Lyapunov function or to estimate the Potential function (see Materials and Methods for details). Traditional DST analysis focuses on the existence and local properties of a given steady state, but it does not explain the global dynamics https://doi.org/10.1016/j.neuroimage.2021.118551 . Received 24 March 2021; Received in revised form 22 August 2021; Accepted 1 September 2021 Available online 8 September 2021. 1053-8119/© 2021 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 on complex systems in which transitions between attractors constitute the characteristic behavior. Brain dynamics do not converge or stabilize around a fixed set of invariants and might be described as a continuous flow of quick and irregular oscillations ( Deco and Jirsa, 2012; Golos et al., 2016 ). Here, we propose to characterize the transitions between attractors by means of a non-stationary attractor landscape , that is, how the set of steady states and their connections evolve with time. Thus, a mathematical method to approximate the non-stationary attractor landscape is developed. To do this, we must rethink one of the central questions in computational neuroscience: the relationship between the human brain as a complex empirical system and the theoretical models expressed by means of equations. Usually, the parameters of the models are fitted in order to simulate general characteristics of the brain such as its functional connectivity or its metastability. However, whatever the values of the fitted parameters are, the attractor landscape offered by these models is static and the empirical transitions between attractors are still not well characterized. Here we propose a new relationship between the model and the empirical data by means of the model transform (MT), which constitutes a useful tool to approximate the non-stationary landscape of complex systems. The MT is a newly developed concept ( Galadí, 2020 ) which allows for the application of simple theoretical models to complex empirical systems in each short interval of time, and can be seen as a generalization of trivial concepts (such as the instantaneous velocity, the curvature or the torsion) when trivial models (the uniform motion in a straight line, the circumference or the helix) are applied in a small neighbourhood of the trajectory (see Materials and Methods for more details). In summary, we search for the simplest model that will help us find an answer to the questions “towards which state is the brain attracted to at a given time point? ”and “how do the corresponding attractor landscape surroundings look like? ”, and we do so with help of DST. In the section Materials and Methods, we show that the model obtained is well known and called the Lotka-Volterra (LV) model. Although it is not our intention to model neural activity, the LV equations have been previously used in neuroscience for this purpose: a LV equation for mean firing rate was derived from the conventional membrane dynamics of a neural network with lateral inhibition and self-inhibition in Fukai and Tanaka (1997) ; LV equations have been used to generate reproducible transient sequences in neural circuits ( Afraimovich et al., 2004 ). It has been also shown that LV equations are capable of representing switching dynamics between different states of neural networks ( Cardanobile and Rotter, 2011 ). The joint activity dynamics of excitatory and inhibitory populations has been analyzed employing a pair of mutually interacting nonlinear differential equations. In absence of a voltage leak for individual neurons and for negligible synaptic transmission delay, these equations take the form of LV equations (Lagzi et al., 2019; Lagzi and Rotter, 2015) . The Lotka-Volterra Transform (LVT) is the particular case of the MT when the harnessed model is the LV model. Thus, assuming time-varying growth rate parameters in the LV equations we obtain a time-varying AL (see Materials and Methods). The LVT can also be seen as a mathematical operator defined to exactly reproduce the empirical BOLD fMRI signals by finding the growth rate function that tracks their evolution over time. However, it is worth noticing that the equations are used to approximate the temporal evolution of the global attractor and the nonstationary attractor landscape, but not to model the temporal evolution of the system itself. In order to explain the difference between modeling and the model transform, it should be underlined that each specific MT uses a specific model, but it is not modeling the system. Therefore, we do not attempt to make predictions of the dynamics in the short or long term. MT renounces the predictive power in exchange for defining new variables. Here, the final objective of MT is to define new measures, in this case measures of consciousness. MT is a generalization of the fact that there are simple models that extract certain information from a trajectory. One trivial example is the circumference, the simplest model which allows to calculate the curvature (the inverse of the radius of curvature) of any smooth curve at any point. But this does not mean that the whole smooth curve is shaped by a circumference. In this work, we do not state that the BOLD signals follow a LV system of differential equations. We just use a transformation to get well characterized non-stationary attractors. Non-stationary attractors should distinguish between awake and deep sleep states as a first partial validation to create a new measure of consciousness. This innovative view relies on the fact that we are computing an asymptotic attractor at every time instant, i.e., the global attractor the system would achieve in the limit assuming a particular value of the growth rate parameter of the LV system. This is a more dynamical view of the problem not previously used in the literature, which let us adapt to the rich behaviour of brain dynamics. Precisely by characterizing the AL by means of the number of energy levels, we are able to identify the different brain states. Over the last couple of years, there has been an increasing interest in trying to identify the necessary and sufficient properties of different brain states ( Deco et al., 2019 ). Yet, thus far, there has been precious little progress in developing the necessary mathematical tools. Here, we propose the innovative nonstationary landscape approach to be assessed distinguishing between different human brain states as measured with neuroimaging. Developing such a framework would be a major step forward, potentially leading not only to a deeper understanding of the AL in different brain states but also to create sensitive and specific biomarkers characterising the brain states of individuals in wakefulness, deep sleep, anesthesia, or different levels of coma ( Deco and Kringelbach, 2014 ). Obviously, this formalism is not needed to discriminate wakefulness and deep sleep, but the AL is of interest in itself. Our main aim is not to provide a better method to discriminate or to describe states of awake and sleep than those already existing in the scientific literature, but to validate, for the first time, a theoretical proposal based on the AL for the study of states of consciousness. Our method is justified as a first attempt to approximate the state the brain is drawn to at each time point and to calculate the corresponding non-stationary AL. In this sense, we are also interested here in characterising human sleep as a previous step towards future clinical applications to coma patients. Human sleep is traditionally subdivided into different stages that alternate in the course of the night ( Nir et al., 2013 ), mainly non-rapid- eye-movement (NREM) and rapid-eye-movement (REM) sleep. NREM is further subdivided into light sleep (N1), NREM sleep stage (N2) and NREM deep sleep stage (N3). From N1 to N3, traveling brain waves become slower and more synchronized. In this study, we aim to 1) associate brain states with the average and variability of the number of energy levels (NoEL) of the AL across time and to 2) demonstrate that this mathematical formalism can be applied to distinguish quantitatively and rigorously between the different brain states of wakefulness and deep sleep using state-of-the-art neuroimaging data collected from healthy human participants. 1.1. Motivations The first motivation for this work is based on the thesis that the fundamental information of a dynamical system is expressed in the structure of its global attractor. One consequence of this thesis is that, in the case of the human brain as a dynamical system, the structure of the attractor could be related to the corresponding states of consciousness. Hence the importance of finding a method to approximate the attractor of brain dynamics and its testing in different states of consciousness. Note that we go from BOLD fMRI signals to a flow of ALs and its abstract structures. We consider of high interest to show that this flow of structures is informative about states of consciousness. Thus, another objective of our method is to define new measures of consciousness. Although we are not trying to improve other current methodologies for classification to distinguish among states of consciousness, we believe 2 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 Fig. 1. Flowchart illustrating the methods. For this study, we used an automated anatomical labelling (AAL) atlas. The Lotka-Volterra transform, constrained by the whole-brain model, is a mathematical operator that calculates the growth rate as a function that exactly reproduces the filtered empirical BOLD fMRI signals. Then, solving the linear complementarity problem (LCP) we calculate the globally asymptotically stable stationary solution (GASS) and the number of energy levels (NoEL) of the attractor landscape at each time instant. In order to assess how different the distributions of the NoEL for awake and for deep asleep are, we calculate a statistical hypothesis test (Wilcoxon test) and the 𝐽 ind (see Supplementary Materials for a more detailed flowchart) before classification. this new perspective opens the door to further detailed studies, so that it could serve as the seminal work for a full research plan. The novelty of this approach is essentially a first transformation of BOLD data into a non-stationary description of the brain attractor landscape by a dynamical system approach. Even beyond neuroscience and on the side of mathematical DST, this way to interpret attractors as an instantaneous object is also new. Another of the motivations for this work is the search for a method to calculate the empirical “ghost ” attractors of the resting (and sleep) state suggested in Deco and Jirsa (2012) , Vohryzek et al. (2020) . The resting state dynamics needs to be explained, and these attractors would correspond to distinct foci of high activity in particular brain areas. According to Deco and Jirsa (2012) , at the edge of the transition from resting state to task context the local attractors would not exist as stable fixed points yet, since they are either saddle points, or regimes with close to zero flow in the phase space. However, as a possible neurobiological explanation of the resting state dynamics, these states could be easily stabilized when needed in a given task context or for a given function. We believe that our new formalism could be used to detect empirical ghost attractors since, in a non-stationary AL, stable attractors are never reached. Thus, they are candidates to play the role proposed in Deco and Jirsa (2012) in the sense that brain dynamics seem to be attracted to them, which can further be interpreted as latently attracting the system and potentially stabilizable in transitions from resting state to specific tasks. 2. Materials and methods 2.1. Overview of the method We aim to characterize the different brain states of wakefulness and deep sleep by approximating the non-stationary attractor landscape (AL) as shown in Fig. 1 with the flowchart of steps. Briefly, first we enter in the Lotka-Volterra (LV) equations 𝑢 𝑖 = 𝑢 𝑖 ( 𝛼𝑖 − 𝑢 𝑖 + 𝑛 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑢 𝑗 ) , 𝑖 = 1 , …, 𝑛, the structural connectivity matrix (SCM) between 𝑛 = 90 standardized brain areas through an automated anatomical labelling (AAL) atlas ( Tzourio-Mazoyer et al., 2002 ) as the 𝛾𝑖𝑗 parameters, which are obtained by applying tractography algorithms to diffusion magnetic resonance imaging (dMRI). Second, we applied the Lotka-Volterra Transform (LVT) 𝛼𝑖 ( 𝑡 ) = 𝑢 𝑖 ( 𝑡 ) 𝑢 𝑖 ( 𝑡 ) + 𝑢 𝑖 ( 𝑡 ) − 𝑔 𝑛 ∑ 𝑗=1 𝛾𝑖𝑗 𝑢 𝑗 ( 𝑡 ) , 𝑖 = 1 , …, 𝑛. to calculate the growth rate 𝛼( 𝑡 ) that reproduces exactly the filtered empirical BOLD fMRI signals 𝑢 𝑖 ( 𝑡 ) . Both 𝑢 ( 𝑡 ) and 𝑢 ( 𝑡 ) are empirical values in the form of discrete time series. The global coupling strength parameter 𝑔will be fitted to balance information from the proper region 𝑖 , 𝑢 𝑖 ( 𝑡 ) 𝑢 𝑖 ( 𝑡 ) + 𝑢 𝑖 ( 𝑡 ) , and information from the rest of the brain ∑𝑛 𝑗=1 𝛾𝑖𝑗 𝑢 𝑗 ( 𝑡 ) . Third, we characterize the corresponding attractor landscape computing the globally asymptotically stable stationary solution (GASS) and the number of energy levels (NoEL) at each time instant 𝑡 ∗ solving 𝑢 𝑖 = 𝑢 𝑖 ( 𝛼𝑖 ( 𝑡 ∗ ) − 𝑢 𝑖 + 𝑛 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑢 𝑗 ) , 𝑖 = 1 , …, 𝑛. Finally, in order to assess differences in the ALs between wakefulness and deep sleep, we use the non-parametric Wilcoxon signed-rank test and the J index. We consider two quantitative descriptors of the ALs, the time average of NoEL ( 𝑞 ) and the standard deviation of NoEL ( 𝜎𝑞 ), both calculated for each participant in the two brain states. We obtained the 𝑝 -value of a paired and two-sided test for the null hypothesis that the distribution of average NoEL (or their standard deviation) presented the same median during awake vs. deep asleep subjects. In addition, we defined the following measure 𝐽 ind : 𝐽 ind = 𝑛𝑠 ∑ 𝑖 =1 ( 𝑞 𝑎𝑠,𝑖 − 𝑞 𝑎𝑤,𝑖 ) 𝑛𝑠 ∑ 𝑖 =1 ∣𝑞 𝑎𝑠,𝑖 − 𝑞 𝑎𝑤,𝑖 ∣ (1) where 𝑛𝑠 is the total number of subjects, 𝑞 𝑎𝑠,𝑖 is the value of 𝑞 for subject 𝑖 while asleep and 𝑞 𝑎𝑤,𝑖 is the value of 𝑞 for subject 𝑖 while awake. Then, 𝐽 ind =1 ( 𝐽 ind =−1 ) if the mean number of energy levels is larger (smaller) for deep sleep than wakefulness for all subjects and 𝐽 ind ≃0 if the null hypothesis holds. Please note that 𝐽 ind can also be defined for 𝜎𝑞 or for any other parameter (see Supplementary Materials for more details). 2.2. BOLD signal acquisition In order to compare different naturally occurring brain states we use data from 18 healthy participants: specifically, BOLD fMRI signals in resting state and asleep phases (N1, N2 and N3). Empirical data comes from a set of fifty-five subjects (thirty-six females, mean ± standard deviation age of 23.4 ± 3.3 years) who fell asleep during a simultaneous EEG- fMRI recording previously described in Tagliazucchi and Laufs (2014) , where 18 participants who reached stage N3 sleep (deep sleep) were selected. The mean duration ( ± standard deviation) of contiguous N3 sleep epochs for these participants was 11.67 ± 8.66 min. fMRI data was recorded at 3T (Siemens Trio, Erlangen, Germany) simultaneously with EEG data using an MR-compatible EEG cap (modified BrainCapMR, Easycap, Herrsching, Germany), MRI-compatible amplifiers (BrainAmp MR, BrainProducts, Garching, Germany), and sleep stages were scored manually by an expert according to the AASM criteria (AASM, 2007). 2.3. BOLD signal preprocessing fMRI data was realigned, normalized and spatially smoothed using SPM8 ( www.fil.ion.ucl.ac.uk ), (see Tagliazucchi and Laufs, 2014 for 3 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 full acquisition, pre-processing and sleep scoring details). Cardiac (with frequencies around 1 ∼2 Hz), respiratory (around 0 . 3 Hz), motioninduced and their temporal aliasing noises were regressed out from the fMRI BOLD signals, and data were band-pass filtered in different ranges. Empirical BOLD signals are generally band-pass filtered to remove the contribution of noise. In this study we take a particular look at this process to better illustrate how the attractor landscape found in resting state and in the deep sleep differentiate from each other. Different filters have been proposed in the literature to obtain the most reliable results. We consider all the possible filters below the Nyquist frequency to show how the non-stationary landscapes differ in awake and asleep conditions. 2.4. Structural connectivity matrix For the whole-brain network model, the interactions between the 90 brain areas were scaled in proportion to their white matter structural connectivity. For this study, we used the structural connectivity between the 90 automated anatomical labelling (AAL) regions obtained in a previous study ( Deco et al., 2017 ) averaged across 16 healthy young adults (5 females, mean ± SD age: 24.75 ± 2.54). Briefly, for each subject, a 90 ×90 structural connectivity matrix Γ= [𝛾𝑖𝑗 ]was obtained by applying tractography algorithms to Diffusion magnetic resonance imaging (dMRI) where the connectivity 𝛾𝑖𝑗 between regions 𝑖 and 𝑗was calculated as the proportion of sampled fibers in all voxels in region 𝑖 that reach any voxel in region 𝑗. Since dMRI does not capture fiber directionality, 𝛾𝑖𝑗 was defined as the average between 𝛾𝑖𝑗 and 𝛾𝑗𝑖 . Averaging across all 16 participants resulted in a structural connectivity matrix Γ representative of healthy young adults. We ordered the different brain areas in the neuroanatomical connectivity matrix in such a way that homotopic regions in the two cerebral hemispheres were arranged symmetrically with respect to the center of the matrix (see Supplementary Materials for details). 2.5. Global attractor and energy levels The global attractor (see Supplementary Materials for a formal definition) is usually comprised of isolated invariant subsets and trajectories connecting them ( Carvalho et al., 2012 ). The highest energy level is defined as formed by isolated invariant subsets that receive no solution, i.e. unstable or source sets. Then, each successive lower level is defined including those isolated invariant subsets that receive solutions only from the previously defined higher levels ( Aragao-Costa et al., 2012 ). These concepts can be illustrated using the Lotka-Volterra (LV) system as follows: 𝑢 𝑖 = 𝑢 𝑖 ( 𝛼𝑖 − 𝑢 𝑖 + 𝑛 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑢 𝑗 ) , 𝑖 = 1 , …, 𝑛, (2) where 𝛼𝑖 are the growth rates and 𝛾𝑖𝑗 the structural connectivity, both considered constant in time. This is an example of a dynamical system on a network where the conditions for existence and uniqueness of solutions are well-known ( Murray, 2013 ). An shown in the example of Fig. 2 A–C with 𝑛 = 2 , with suitable 𝛾𝑖𝑗 and 𝛼𝑖 there are four stationary points (i.e. four isolated invariant subsets): one located at the origin ( 𝑢 1 = 0 , 𝑢 2 = 0) , one at the x-axis ( 𝑢 1 ≠0 , 𝑢 2 = 0) , one at the y-axis ( 𝑢 1 = 0 , 𝑢 2 ≠0) , and one at the positive quadrant ( 𝑢 1 ≠0 , 𝑢 2 ≠0) . The stationary points of the LV system for 𝑛 = 2 can be described as two binary variables, null values indicating that the corresponding dynamical variable is zero (i.e. is located in one of the axes), and non-null values indicating that it is located outside that axis. In other words, the global attractor of the system (and thus its asymptotic behavior) could be encoded as a set of Boolean variables and connecting global solutions among them. Figure 2 C shows that the stationary points are partially ordered, i.e. they admit an order relationship in which not necessarily all the pairs of elements can be compared. In this ordering, the unstable and stable points are maximal and minimal respectively, with the two other stationary points (saddle points) not being ordered with respect to each other. Thus, there are three energy levels: the trivial solution (unstable point), the saddle points, and the stable solution. The number of energy levels (NoEL) 𝑞will be the most important parameter in our data analysis (see Supplementary Materials for more details and formal definitions). Figure 2 D–F show a LV system for 𝑛 = 4 where 𝛼𝑖 𝑖 = 1 , …, 4 , are now periodic functions of time. Changes in the 𝛼𝑖 produce changes in the AL so that the NoEL 𝑞also changes over time. For LV systems, each energy level is formed by stationary points with the same number of non-zero components, and the only stable point is a single stationary point in the lowest level with incoming solutions only, called globally asymptotically stable solution (GASS). The GASS is the point the system is attracted to, since any initial point of ℝ 𝑛 + will converge to the GASS ( Takeuchi and Adachi, 1980 ). It can be shown that there exists an equivalence between obtaining the GASS of LV systems and solving a linear complementarity problem (LCP, see Takeuchi, 1996 and Supplementary Materials). Furthermore, the NoEL 𝑞equals the number of non-zero entries in the GASS plus one (see Fig. 2 F). The LCP associated with the LV can be solved using the Complementary Pivot Algorithm ( Cottle et al., 1992 ). Although the LV model is relatively simple, with 𝑛 = 90 nodes it could include a complex AL with up to 2 90 ( 1 . 238 ×10 27 ) stationary points. The Informational Structure (IS) of the global attractor is defined as a directed graph composed of nodes associated with the isolated invariant subsets and links establishing their connections (see Fig. 2 C and F, and Supplementary Materials for a formal definition). 2.6. Model transform A formal definition of the Model Transform (MT) is provided in the Supplementary Materials. In short, the idea is the same as when a uniform motion in a straight line is fitted to any three-dimensional motion in an arbitrarily small segment of the trajectory. If we want to quantify how a motion 𝑥 ( 𝑡 ) “moves away from rest ”, we must look for the simplest possible motion that allows us to measure that distance. Simple motion in this case is the uniform motion in a straight line and, adjusting it to each small neighborhood of 𝑥 ( 𝑡 ) , the variation per unit of time of 𝑥 ( 𝑡 ) is a measure of “remoteness from rest ”. Therefore, we say that the uniform motion in a straight line is the minimum model that transforms a motion 𝑥 ( 𝑡 ) into 𝑣 ( 𝑡 ) , or conversely that 𝑣 ( 𝑡 ) is the MT of 𝑥 ( 𝑡 ) when the minimum model is uniform motion in a straight line. Note that first, we have an intuitive idea of what we want to measure and then, after choosing a minimal model, the intuitive idea is specified and formalized. Thus, the equations 𝑥 ′ 𝑖 ( 𝑡 ) = 𝑣 𝑖 , 𝑖 = 1 , 2 , 3 can be seen as a “model ”of uniform motion in a straight line when 𝑣 𝑖 has a fixed value for 𝑖 = 1 , 2 , 3 or, alternatively, they can be seen as a definition of the instantaneous velocity 𝑣 ( 𝑡 ) at each time point 𝑡 for any trajectory 𝑥 ( 𝑡 ) . In the latter case, the value of 𝑣 changes over time and the MT transforms 𝑥 ( 𝑡 ) into 𝑣 ( 𝑡 ) . The 𝑣 ( 𝑡 ) can also be seen as defined to exactly reproduce 𝑥 ( 𝑡 ) integrating 𝑥 ′ 𝑖 ( 𝑡 ) = 𝑣 𝑖 ( 𝑡 ) , 𝑖 = 1 , 2 , 3 when the initial point 𝑥 0 = 𝑥 ( 𝑡 0 ) is known. The uniform motion in a straight line is the simplest model that allows us to know “how far an object moves from rest at each time instant ”using the model parameters 𝑣 𝑖 , 𝑖 = 1 , 2 , 3 . In summary, we look for the simplest model that allows us to answer the questions “towards which state is the brain attracted to at one time point? ”and “how is the attractor landscape characterized at one time point? ”. These questions are answered with help of the DST. Mathematically speaking, the Lotka-Volterra model is the simplest nontrivial cooperative model which includes the trivial solution (0 , 0 , …, 0) and possesses a GASS. Models of global brain dynamics are usually cooperative in the sense that different brain areas act cooperatively (excitation), while there is inhibition only at local level (within each region). The inclusion of the trivial solution and the GASS are crucial in our study of the attractor structure, since they both mark the beginning and the end 4 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 Fig. 2. ( A ), Two-node structural network of a Lotka-Volterra (LV) system (2) , with 𝛼1 = 0.71, 𝛼2 = 2.85, 𝛾12 = 0.2 and 𝛾21 = 0.1. ( B ), Solutions of the LV system in A . The global attractor is the set of all bounded solutions (yellow area).( C ), The structure of the global attractor is a new network made by four nodes associated with the four stationary points and directed links associated with connecting solutions. The relation induced by the links is transitive such that only the minimal links are represented. Each node of this graph is represented as a subgraph of the original LV system where non-null components are shown in black and null in grey. In this example there are three energy levels: the trivial solution, the saddle points, and the stable solution. ( D ), Structural network of a 4- dimensional ( 𝑛 = 4 ) LV system (see Supplementary Materials for 𝛾values). ( E ), In this example 𝛼𝑖 are periodic functions of time. ( F ), Structures of the global attractor corresponding to time steps shown in E . The number of energy levels (NoEL) changes over time, and equals the number of non-zero entries in the lowest level plus one. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) of said structure; that is, both the first and the last energy levels, along with the “distance ” between them give us the NoEL. Furthermore, we show (see Supplementary Materials) that the Lotka- Volterra model: i) includes a rich attractor landscape with a large number of stationary points, which are a unique combination of active and inactive nodes of the dynamical system, ii) includes empirical information regarding the structural connectivity of the dynamical system, iii) has a well-known condition for existence and uniqueness of the state the system is attracted to, i.e. the GASS, iv) supports a Lyapunov function so that its isolated invariants could be ordered according to the energy levels, v) supports an algorithm of low computational complexity (the Complementary Pivot Algorithm for solving the LCP), so that vi) the global attractor has a well-known structure. The Lotka-Volterra Transform (LVT) is the particular case of MT when the model is given by the LV equations: 𝑢 𝑖 = 𝑢 𝑖 ( 𝛼𝑖 − 𝑢 𝑖 + 𝑔 𝑛 ∑ 𝑗=1 𝛾𝑖𝑗 𝑢 𝑗 ) , 𝑖 = 1 , …, 𝑛, (3) where 𝑔is a global coupling strength parameter. From these equations we can obtain time-dependent parameters 𝛼𝑖 ( 𝑡 ) as follows: 𝛼𝑖 ( 𝑡 ) = 𝑢 𝑖 ( 𝑡 ) 𝑢 𝑖 ( 𝑡 ) + 𝑢 𝑖 ( 𝑡 ) − 𝑔 𝑛 ∑ 𝑗=1 𝛾𝑖𝑗 𝑢 𝑗 ( 𝑡 ) , 𝑖 = 1 , …, 𝑛. (4) This expression defines the LVT. Finally, we can associate each parameter and variable of the LVT with its neurobiological counterpart: both 𝑢 ( 𝑡 ) and 𝑢 ( 𝑡 ) will be empirical values from BOLD fMRI signals in discrete form; structural connectivity information in 𝛾𝑖𝑗 will be obtained by applying tractography algorithms to diffusion magnetic resonance imaging (dMRI) where the connectivity 5 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 between regions 𝑖 and 𝑗is calculated as the proportion of sampled fibers in all voxels in region 𝑖 that reach any voxel in region 𝑗; the global coupling strength parameter 𝑔will be fitted to balance information from the proper region 𝑖 , 𝑢 𝑖 ( 𝑡 ) 𝑢 𝑖 ( 𝑡 ) + 𝑢 𝑖 ( 𝑡 ) , and information from the rest of the brain ∑𝑛 𝑗=1 𝛾𝑖𝑗 𝑢 𝑗 ( 𝑡 ) . LVT transforms 𝑢 ( 𝑡 ) into 𝛼( 𝑡 ) , i.e., LVT transforms empirical data into auxiliary functions 𝛼𝑖 ( 𝑡 ) . 2.7. Alternative model transforms Given a dataset like ours, with resting state and N3 deep sleep data, can any model transform differentiate between these states? In the search for biomarkers of consciousness, what does the LVT contribute with, compared to any other MT? We use four generative models possessing various levels of complexity: i) the Kuramoto model for coupled oscillators, ii) the simultaneous autoregressive (SAR) model, a purely spatial model with no dynamics that expresses BOLD fluctuations within one region as a linear combination of the fluctuations in other regions, iii) a rate fluctuations model, which is a simplified version of the Wilson–Cowan system considering exclusively the excitatory population, and iv) the Lotka-Volterra (LV) model. The Kuramoto model Kuramoto (1975) , Yeung and Strogatz (1999) consists of a set of coupled oscillators: 𝑑𝜙𝑖 ( 𝑡 ) 𝑑𝑡 = 𝜔 𝑖 + 𝑔 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑠𝑖𝑛 ( 𝜙𝑗 ( 𝑡 ) − 𝜙𝑖 ( 𝑡 )) where 𝜙𝑖 and 𝜔 𝑖 stand for the phase and intrinsic angular frequency of region 𝑖 . 𝑔and 𝛾𝑖𝑗 are defined as in the LV model. The simultaneous autoregressive (SAR) model is composed of a linear combination of the fluctuations within other regions ( Tononi et al., 1994 ): 𝑢 𝑖 = 𝑔 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑢 𝑗 + 𝜈𝑖 where 𝜈𝑖 stands for uncorrelated white Gaussian noise, and 𝑢 𝑖 , 𝑔and 𝛾𝑖𝑗 are defined as in the LV model. The Rate fluctuations model ( Galán, 2008 ) is a simplification of the Wilson–Cowan model where inhibitory neurons and saturation have been removed: 𝜏𝑑𝑢 𝑖 ( 𝑡 ) d 𝑡 = − 𝑢 𝑖 ( 𝑡 ) + 𝑔 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑢 𝑗 ( 𝑡 ) + 𝜈𝑖 . Here, 𝜏is the time scale of the excitatory population, and 𝜈𝑖 stands for uncorrelated white Gaussian noise. 𝑢 𝑖 , 𝑔and 𝛾𝑖𝑗 are defined as in the LV model. Specifically, to test which model stands out compared to the others when it comes to finding differences between brain states, we contrast the four different model transforms (MT) to validate the choice of LV. Always starting from the same dataset, we look for the model with greatest capacity to distinguish between states of consciousness. In order to test this capacity, we add the sign test and the paired t -test (for more details see Supplementary Materials) to the non-parametric Wilcoxon signed-rank test and the J index. The models that we compare with LV have shown high predictive power of the functional connectivity starting from the structural connectivity ( Messe et al., 2015 ). To calculate the Kuramoto Transform (KT), we previously perform a Hilbert transform ( Marple, 1999 ) of the BOLD signal to obtain the amplitude and phase of the signal. 𝜙𝑖 is now said empirical phase in each region i. 𝜔 𝑖 is the parameter of the model that becomes a function of time and that will be the actual result of the KT: 𝜔 𝑖 ( 𝑡 ) = 𝑑𝜙𝑖 ( 𝑡 ) 𝑑𝑡 − 𝑔 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑠𝑖𝑛 ( 𝜙𝑗 ( 𝑡 ) − 𝜙𝑖 ( 𝑡 )) . For the SAR Transform (SART), 𝜈𝑖 is the model parameter that is converted to a function of time. Therefore, instead of being modeled as a Gaussian noise, it is calculated exactly from the experimental data: 𝜈𝑖 ( 𝑡 ) = 𝑢 𝑖 − 𝑔 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑢 𝑗 In the Rate fluctuations Transform (RFT), 𝜈𝑖 is again the chosen parameter: 𝜈𝑖 ( 𝑡 ) = 𝜏𝑑𝑢 𝑖 ( 𝑡 ) d 𝑡 + 𝑢 𝑖 ( 𝑡 ) − 𝑔 ∑ 𝑗≠𝑖 𝛾𝑖𝑗 𝑢 𝑗 ( 𝑡 ) . For details about these model transform processes, see Supplementary Materials. 3. Results 3.1. Impact of free parameters We transform time series of BOLD signal data of 18 healthy participants into time series of attractor landscapes (ALs), and investigate the distribution of the number of energy levels (NoELs) of the AL during wakefulness and deep sleep. Although the 𝛼𝑖 become time-varying and are fitted in each small time interval, other parameters remain free and we study their impact on the classification of brain states. Impact of coupling strength: The strength of the connections between pairs of nodes is controlled by the coupling strength parameter 𝑔. Changes in coupling strength 𝑔lead to changes in 𝛼, and therefore have an impact on the distribution of NoEL. Figure 3 A shows the impact of 𝑔on the ability of AL to differentiate between the two brain states, i.e., showing how the gap between awake and deep asleep changes as a function of 𝑔. A vertical dashed line shows the maximum value ( 𝑔 < 0 . 36743 ) ensuring global stability. We compute the Pearson linear correlation coefficient between the two curves ( −0 . 9781 ) to corroborate that the 𝐽 ind basically measures the same as the Wilcoxon test (although 𝐽 ind is more intuitive and easier to calculate). The optimal value ( 𝑔 =0 . 29 ) is the same in both tests, and therefore we choose it for our computations. Note however, that the gap between awake and asleep is not restricted to this value. Indeed, the range of 𝑔in which the differences between the distributions of energy are manifest ( 𝑝 -values below 1% or 10 −2 ) is very broad. In Fig. 3 A we can see that there is a wide range of values for which the 𝑝 -value is small (from around 𝑔 =0 . 05 to the red dashed line in 𝑔 =0 . 36743 , where global stability of the system is guaranteed). In that same interval 𝐽 ind remain above 72% ( 𝐽 ind = 0 . 72 ). In other words, our results are robust, as they do not depend on a specific value of 𝑔. In addition we verified ( Fig. 3 A) that the suggested method does not produce good results when the information from the coupling term in the LVT is not considered, i.e. 𝑔tends to zero. At the other extreme, when 𝑔tends to its maximum value the LVT does not distinguish so well between brain states. We can state that our method’s ability to distinguish brain states only decreases when global connectivity is too small or too large. Impact of time series filtering: Fig. 3 B shows the results of the Wilcoxon test and 𝐽 ind index for the mean 𝑞 and the standard deviation 𝜎𝑞 of the NoEL assuming different filter ranges. At high frequencies, noise was not filtered and the ALs of the two brain states were very similar, the 𝑝 - value was high and the 𝐽 ind small. For 𝑞 , the filter range with minimum 𝑝 -value was 0 . 077 − 0 . 096 Hz ( 𝑝 = 0.00023 ). For 𝜎𝑞 , the filter range with minimum 𝑝 -value was 0 . 077 − 0 . 106 Hz ( 𝑝 = 0.00023 ) and the second best was 0 . 077 − 0 . 096 Hz ( 𝑝 -value = 0.00028 ). For 𝜎𝑞 , the filter range with maximum 𝐽 ind was 0 . 077 − 0 . 106 Hz ( 𝐽 ind = 0.9954 ) and the second best was 0 . 077 − 0 . 096 Hz ( 𝐽 ind = 0.9813 ). Therefore, we chose the filtering range 0 . 077 − 0 . 096 Hz. When the signal was filtered in this range, the difference between the distributions of NoEL and therefore, between the ALs of awake and sleeping participants, is maximal. The values obtained for most filters (wide filters, narrow filters or even filter absence, i.e. 0 . 0 − 0 . 24 Hz where 0.24 Hz is the Nyquist frequency) were good enough ( 𝐽 ind > 0 . 75 and 𝑝 -value < 0 . 01 ) to sustain that our results can not be attributed to a specific filter. Rather, we demonstrated that our 6 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 Fig. 3. ( A ), Impact of coupling strength 𝑔on the 𝑝 -value of the Wilcoxon test and the 𝐽 ind maximization, both for mean NoEL. The red dashed line ( 𝑔 < 0 . 36743 ) delimits the region in which global stability is guaranteed. In both cases the optimal value is 𝑔 =0 . 29 . ( B ), Impact of the filter used on the 𝑝 -value of the Wilcoxon test and the 𝐽 ind for the samples of both, 𝑞 and 𝜎𝑞 . Both, the Wilcoxon test 𝑝 -value and the 𝐽 ind are functions of the ends of the filtering range (see Supplementary Materials for a 2D version of the graphics).. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) hypotheses, in general, were true even for the unfiltered signal. Only high-pass filters provided negative results ( 𝐽 ind < 0 . 4 and 𝑝 > 0 . 2 ). Thus, our conclusions are again robust in this regard. 3.2. Statistical analysis of non-stationary attractor landscapes The Lotka-Volterra Transform (LVT) allows us to compute the timevarying vector 𝛼for every participant in both brain states. Using the linear complementarity problem (LCP) we computed the globally asymptotically stable solution (GASS). Figure 4 A shows a 208-seconds long sample of the NoEL as a function of time for one subject in awake state and deep sleep state. We studied the NoEL distributions for the 18 participants in both conditions. Generally, the NoEL increased when the components of alpha grow ( Fig. 4 B, where alpha components were computed using the LVT). In order to draw general conclusions beyond the inter-individual differences, we computed the average distribution of the 18 participants in the awake state and the average distribution of the 18 participants in the sleep state ( Fig. 4 C). The pattern typically observed was that the distribution for the awake state was more homogeneous along 𝑞 ∈[ 1 , 80 ] NoEL values. ALs with extreme NoEL values (that is, low: < 28 , or very high: 90 or 91), were more frequent in awake individuals. Figure 4 D shows the mean NoEL 𝑞 in awake vs. deep asleep (N3) with filtered signals in the range 0 . 077 − 0 . 096 Hz. The error bars were computed according to 𝑒𝑟𝑟𝑜𝑟 = 𝜎𝑞 ∕ √𝑇 where 𝜎𝑞 is the standard deviation of the NoEL and 𝑇 was the length of the data series. For all participants except one, the mean NoEL 𝑞 was higher when the participants were in deep sleep. In order to assess how different 𝑞 is for participants in awake and deep sleep, we use the signed-rank test Wilcoxon and we obtain 𝑝 = 0.00023 , i.e., differences between the means of both brain states are significant. In addition, we compute 𝐽 ind =0 . 9948 . Recall that 𝐽 ind =1 when 𝑞 is larger for all participants in the deep sleep state. Then, we establish patterns regarding the variability of the distributions. We compute the standard deviations of the distribution of NoEL in each participant for the two brain states. In Fig. 4 E each point represents a participant. The error bars were calculated according to the standard error of the standard deviation 1 2 𝜎√ 1 𝑇 ( 𝜇4 − 𝑇 − 3 𝑇 − 1 𝜎4 ) (5) where 𝜎is the standard deviation of the NoEL, 𝑇 is the length of the data series, and 𝜇4 = 𝐸( 𝑞 − 𝑞 ) 4 Rao (2009) . For 17 (out of 18) participants the standard deviation is larger for wakefulness than for the deep sleep state, 𝐽 ind = −0 . 9813 and 𝑝 = 0.00028 . Again, the differences between the standard deviations of both brain states are very significant. There is only one participant who shows a greater deviation when asleep. This may be more clearly understood recalling Fig. 4 C, where ALs with extreme NoEL values are more frequent in awake individuals. 3.3. Comparison including N1 and N2 Results of the analysis of data from N1 and N2 states are included in Fig. 5 A–C. where N1, N2 and N3 states are significantly ordered according to statistical tests ( Fig. 5 D) so that, as deeper levels of sleep are reached, the average NoEL increases while the standard deviation decreases. However, the differences between awake and N1 are not significant neither for the mean nor for the standard deviation ( Fig. 5 D-E). The same is true for the comparison between awake and N2 ( Fig. 5 D and F). Nevertheless, the NoEL distributions distinguish between these 7 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 Fig. 4. ( A ), Sample of 208 seconds of the number of energy levels (NoEL) for participant number 5 in wakefulness and deep sleep states (filter 0 . 01 − 0 . 1 Hz, 𝑔 =0 . 29 ). ( B ), Number of Energy Levels (NoEL) as a function of the average 90 alpha components (same 𝑔as before). Each dot represents a time point of a subject in a brain state. ( C ), Average distribution of the NoEL for 18 healthy participants in awake and in deep sleep states (see Supplementary Materials for individual histograms). ( D ), Average number of energy levels 𝑞 of 18 participants awake vs. deep asleep (N3) ( 𝐽 ind =0 . 9948 , 𝑝 = 0.00023 ). 𝑞 is the same for wakefulness and for deep sleep brain states along the red line. ( E ), Standard deviation of the NoEL for deep asleep state vs. wakefulness state. ( 𝐽 ind = −0 . 9813 , 𝑝 -value = 0.00028 ). In graphics B-E the NoEL was obtained from filtered data between 0 . 077 − 0 . 096 Hz (see Supplementary Materials for a similar figure using other filters).. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) states: Fig. 5 G shows the NoEL distribution averaged across the four states (awake, N1, N2 and N3) and the 18 subjects (inset), as well as the differences between the average distributions for each state with respect to the mean. We observe that the differences between states are more significant for certain NoELs. Thus, if we compare the frequency at which the AL reaches 15 levels in each of the 18 subjects for awake and N1 states, the differences are significant ( Fig. 5 H). The same is true when we compare the frequency at which the AL reaches 13 levels in the comparison between awake and N2 states ( Fig. 5 I). Further observations show that the average NoEL for five subjects shows a monotonic increase from awake to N3 sleep, whereas another four subjects shows the same monotonic increase but with lower average NoEL in N1 than in the awake state (see Supplementary Materials). 3.4. Comparison of model transforms We contrast different model transforms (MT) and show that LV stands out compared to other models when it comes to finding differences between brain states ( Fig. 6 A–E). The MTs for four different models (LV model, SAR model, Rate fluctuations model and Kuramoto model) are compared and validate the choice of LV. We use the same dataset for all MTs and start filtering the BOLD signals in the range 0 . 077 − 0 . 096 𝐻𝑧 . 𝑔was optimized in each MT to find the greatest differences between brain states. The results show that LV has a greater capacity to distinguish between states of consciousness than the other models. In our formalism, we transform the BOLD signals into alpha values using the LVT and then use these alphas to calculate the structure of the corresponding attractor, along with its energy levels. Figure 6 F compares the ability of alphas with the ability of NoEL to distinguish states of consciousness. The result is that the structure of the attractor through its energy levels is able to distinguish between states much better, so the simple calculation of the alphas would not be enough to generate a measure of consciousness. 3.5. Classification In addition, we have classified the subjects by training a KNN classifier with 𝑘 = 2 , obtaining high percentages of accuracy with minimal amounts of training data. Leave-one-out cross-validation and Receiver Operating Characteristics (ROC) analysis are included in this classification process. Figure 6 G–H shows the classification results training the classifier with the mean and the standard deviation of the NoEL of 36 subjects (18 in N3 and 18 in awake resting state). The area under the curve (AUC = 0 . 89 ) is a performance measure of our classification problem at various threshold values: it indicates how much our formalism is capable of distinguishing between brain states. The higher the AUC, the better the formalism to predict N3 subjects as being in N3 and awake subjects being awake. By analogy, the higher the AUC, the better to distinguish between brain states. Thus, subjects can be classified with an accuracy of 89% by training the classifier with only the standard deviation of the NoEL. If we add 8 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 Fig. 5. ( A–F ), Comparisons of NoEL standard deviation among awake, N1, N2 and N3 (deep asleep) states (see Fig. 4 for awake-N3 comparison). Error bars were calculated according to (5) ( D ), 𝐽 ind and 𝑝 -values of the Wilcoxon tests. The differences between N1, N2 and N3 were very significant ( 𝑝 -values < 5% ), but differences between awake and N1 or N2 were not significant (see Supplementary Materials for comparisons of average NoEL). ( G ), Deviation of the NoEL distribution of each state respect to average NoEL distribution (inset) of the four states. ( H ), Frequency of attractor landscapes (ALs) of 15 energy levels of 18 participants awake vs. N1 ( 𝐽 ind =0 . 6629 , 𝑝 = 0.0139 ). ( I ), Frequency of ALs of 13 energy levels of 18 participants awake vs. N2 ( 𝐽 ind =0 . 6418 , 𝑝 = 0.0139 ). In all graphics (A-I) the NoEL was obtained from filtered data between 0 . 077 − 0 . 096 Hz. the mean NoEL, the accuracy rises to 92% as shown in Fig. 6 G. It is noteworthy that the classifier can be trained with particular values of the NoEL distribution (instead of global characteristics of the distribution like mean and standard deviation) such as the frequency of a given NoEL level. For instance, with the frequency of level 59, an 89% accuracy is obtained. Moreover, if we add the frequency of level 3, the accuracy rises further to 92%. 3.6. Attracting brain areas Next, we investigated which areas were active (non-zero components) in the attractor (GASS). On average, the attractor of the deep sleep state was more populated (larger NoEL) than that of the wakefulness state. Hence, most of the brain areas tended to appear more frequently in the deep sleep attractor than in the wakefulness attractor. In other words, as 𝑞 was greater in the participants in deep sleep, the general tendency was that each area of the brain was more frequently present in the deep sleep attractors than in the wakefulness attractors (with some exceptions). As shown in Fig. 7 AB, we computed the frequency of appearance of an area in the GASS for each participant and for each brain state, averaged over all participants, and then we computed the difference between the frequencies in awake and in deep sleep states for each participant and area ( Fig. 7 C). As can be seen, most of the areas showed less presence in the wakefulness attractor than in the deep sleep attractor. But there were some exceptions (indicated by the warm colours), i.e., greater presence of the area in the wakefulness attractor (see Supplementary Materials for individual differences, where we found inter-individual variability but, in all cases, warm colors were exceptional). Considering the averaged results (see Fig. 7 C) for the 18 participants, 14 areas (out of 90) showed this opposite tendency. 12 out of these 14 are actually six pairs of homotopic zones, i.e., pairs of areas that occupy symmetrical zones in each hemisphere. These six pairs of areas are the olfactory, the calcarine (V1), the cuneus (basic visual), the lingual (visual letters), the occipital superior (visual), and the paracentral lobule (motor-sensory), all of them found in both left and right hemispheres. Furthermore, the fusiform in the left hemisphere (facial recognition) and the supplementary motor area in the right hemisphere (control movement). Due to their specific functions, in all cases, it makes sense that these areas were foci of attraction in the wakefulness state. 4. Discussion In this paper we have proposed a novel mathematical formalism for approximating the time-varying attractors landscape (AL) depicted by the dynamical activity of the human brain. In the rich and complex landscape provided by the Lotka-Volterra (LV) equations for collaborating species, the energy levels of the global attractor can be found by identifying the globally asymptotically stable solution (GASS). Our main aim is not to provide a better method to discriminate or describe states of awake and sleep, but to validate, for the first time, a theoretical proposal for the study of consciousness states. Although we do not pretend to find a better classification method than those that already exist in the scientific literature, our novel method allows us to clearly characterise and distinguish between different brain states by means of the variability of the associated energy levels, opening up for relevant clinical applications. We defined the Lotka-Volterra transform (LVT) as the time-varying growth rate parameter, 𝛼=𝛼( 𝑡 ) , that exactly reproduces the empirical BOLD signal. This allowed us to compare the corresponding non-stationary attractor landscapes in wakefulness and sleep conditions. Lotka-Volterra model was chosen as the simplest model that helped us find an answer to the questions “towards which state is the brain attracted to at a given time point? ”and “how do the corresponding attractor landscape surroundings look like? ”. Therefore, the choice of the model and the time-dependent parameters (rates of time-dependent growth and fixed structural connectivities) ensure maximum simplicity and coherence. Furthermore, the results show that the method is able to characterise and significantly distinguish ( 𝑝 -values < 5%) among the 9 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 63723-P, PGC2018-096540-B-I00, Proyecto US-1254251 del Fondo Europeo de Desarrollo Regional (FEDER) y la Consejería de Economía, Conocimiento, Empresas y Universidad de la Junta de Andalucía, dentro del Programa Operativo FEDER 2014-2020 and Proyecto PAIDI 2020 P20_00592. S.S.P. is supported by the Research Project PGC2018- 096641-B-I00 (Ministerio de Ciencia, Innovación y Universidades / Agencia Estatal de Investigación /FEDER, UE). G.D. is supported by the Spanish Research Project PSI2016-75688-P (Agencia Estatal de Investigacíon/Fondo Europeo de Desarrollo Regional, European Union); by the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreements 720270 (Human Brain Project [HBP] SGA1) and 785907 (HBP SGA2); and by the Catalan Agency for Management of University and Research Grants Programme 2017 SGR 1545. M.L.K. is supported by the European Research Council Consolidator Grant: CAREGIVING (615539) and Center for Music in the Brain, funded by the Danish National Research Foundation ( DNRF117 ). H.L. and data acquisition were supported by the Bundesministerium für Bildung und Forschung (grant no. 01 EV 0703); H.L. and E.T. by the LOEWE Neuronale Koordination Forschungsschwerpunkt Frankfurt (NeFF); H.L. also by Christian-Albrechts-University Kiel and the CRC1261 (Deutsche Fors-chungsge-meinschaft). We thank Astrid Morzelewski for sleep scoring of the EEGs. Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.neuroimage.2021.118551 . References Afraimovich, V., Zhigulin, V., Rabinovich, M., 2004. On the origin of reproducible sequential activity in neural circuits. Chaos 14 (4), 1123–1129. doi: 10.1063/1.1819625 . Andrillon, T., Nir, Y., Cirelli, C., Tononi, G., Fried, I., 2015. Single-neuron activity and eye movements during human rem sleep and awake vision. Nat. Commun. 6 (1), 7884. doi: 10.1038/ncomms8884 . Aragao-Costa, E. , Caraballo, T. , Carvalho, A. , Langa, J. , 2012. Continuity of Lyapunov functions and of energy level for a generalized gradient semigroup. Topol. Methods Nonlinear Anal. 39 (1), 57–82 . Bassett, D., Bullmore, E., 2006. Small-world brain networks. Neuroscientist 12 (6), 512– 523. doi: 10.1177/1073858406293182 . PMID: 17079517 Bortolan, M. , Carvalho, A. , Langa, J. , 2020. Attractors under autonomous and non-au- tonomous perturbations. Mathematical Surveys and Monographs. American Mathematical Society . Cardanobile, S., Rotter, S., 2011. Emergent properties of interacting populations of spiking neurons. Front. Comput. Neurosci. 5, 59. doi: 10.3389/fncom.2011.00059 . Carvalho, A. , Langa, J. , Robinson, J. , 2012. Attractors for infinite-dimensional non-au- tonomous dynamical systems. Applied Mathematical Sciences. Springer New York . Cottle, R. , Pang, J. , Stone, R. , 1992. The linear complementarity problem. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104) . Deco, G., Cabral, J., Woolrich, M., Stevner, A., van Hartevelt, T., Kringelbach, M., 2017. Single or multiple frequency generators in on-going brain activity: a mechanistic whole-brain model of empirical meg data. NeuroImage 152, 538–550. doi: 10.1016/j.neuroimage.2017.03.023 . Deco, G. , Cruzat, J. , Cabral, J. , Tagliazucchi, E. , Laufs, H. , Logothetis, N. , Kringelbach, M. , 2019. Awakening: predicting external stimulation forcing transitions between different brain states. PNAS 116 (36), 18088–18097 . Deco, G., Jirsa, V., 2012. Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J. Neurosci. 32 (10), 3366–3375. doi: 10.1523/JNEUROSCI.2523-11.2012 . http://www.jneurosci.org/content/32/10/3366.full.pdf Deco, G., Kringelbach, M., 2014. Great expectations: using whole-brain computational connectomics for understanding neuropsychiatric disorders. Neuron 84 (5), 892–905. doi: 10.1016/j.neuron.2014.08.034 . Esteban, F., Galadí, J., Langa, J., Portillo, J., Soler-Toscano, F., 2018. Informational structures: a dynamical system approach for integrated information. PLoS Comput. Biol. 14 (9), 1–33. doi: 10.1371/journal.pcbi.1006154 . Fukai, T., Tanaka, S., 1997. A simple neural network exhibiting selective activation of neuronal ensembles: from winner-take-all to winners-share-all. Neural Comput. 9 (1), 77–97. doi: 10.1162/neco.1997.9.1.77 . Galadí, J. , 2020. Dynamical systems applied to consciousness and brain rhythms in a neural network. Department of Differential Equations and Numerical Analysis. University of Seville Ph.D. thesis . Galadí, J., Torres, J., Marro, J., 2020. Emergence and interpretation of oscillatory behaviour similar to brain waves and rhythms. Commun. Nonlinear Sci. Numer. Simul. 83, 105093. doi: 10.1016/j.cnsns.2019.105093 . Galán, R.F., 2008. On how network architecture determines the dominant patterns of spontaneous neural activity. PLoS ONE 3 (5), 1–10. doi: 10.1371/journal.pone.0002148 . Golos, M., Jirsa, V., Dauce, E., 2016. Multistability in large scale models of brain activity. PLoS Comput. Biol. 11 (12), 1–32. doi: 10.1371/journal.pcbi.1004644 . Gugger, M., Bögershausen, S., Schäffler, L., 1993. Arousal responses to added inspiratory resistance during rem and non-rem sleep in normal subjects.. Thorax 48 (2), 125–129. doi: 10.1136/thx.48.2.125 . https://thorax.bmj.com/content/48/2/125.full.pdf Hahn, G., Zamora-López, G., Uhrig, L., Tagliazucchi, E., Laufs, H., Mantini, D., Kringelbach, M.L., Jarraya, B., Deco, G., 2021. Signature of consciousness in brain-wide synchronization patterns of monkey and human fMRI signals. NeuroImage 226, 117470. doi: 10.1016/j.neuroimage.2020.117470 . Hale, J. , 1988. Asymptotic behavior of dissipative systems. Mathematical Surveys and Monographs. American Mathematical Society . Hobson, J.A. , 2014. Introduction. Springer International Publishing, Cham, pp. 3–7 . Jobst, B., Hindriks, R., Laufs, H., Tagliazucchi, E., Hahn, G., Ponce-Alvarez, A., Stevner, A., Kringelbach, M., Deco, G., 2017. Increased stability and breakdown of brain effective connectivity during slow-wave sleep: mechanistic insights from whole-brain computational modelling. Sci. Rep. 7 (1), 4634. doi: 10.1038/s41598-017-04522-x . Kalita, P., Langa, J., Soler-Toscano, F., 2019. Informational structures and informational fields as a prototype for the description of postulates of the integrated information theory. Entropy 21 (5). doi: 10.3390/e21050493 . Kuramoto, Y. , 1975. Self-entrainment of a population of coupled non-linear oscillators. In: Araki, H. (Ed.), International Symposium on Mathematical Problems in Theoretical Physics. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 420–422 . Maki-Marttunen, V., Diez, I., Cortes, J., Chialvo, D., Villarreal, M., 2013. Disruption of transfer entropy and inter-hemispheric brain functional connectivity in patients with disorder of consciousness. Front. Neuroinform. 7, 24. doi: 10.3389/fninf.2013.00024 . Markov, N., Ercsey-Ravasz, M., Essen, D.V., Knoblauch, K., Toroczkai, Z., Kennedy, H., 2013. Cortical high-density counterstream architectures. Science 342 (6158), 1238406. doi: 10.1126/science.1238406 . Marple, L., 1999. Computing the discrete-time “analytic ”signal via fft. IEEE Trans. Signal Process. 47 (9), 2600–2603. doi: 10.1109/78.782222 . Mediano, P.A., Seth, A.K., Barrett, A.B., 2019. Measuring integrated information: comparison of candidate measures in theory and simulation. Entropy 21 (1). doi: 10.3390/e21010017 . Messe, A., Rudrauf, D., Giron, A., Marrelec, G., 2015. Predicting functional connectivity from structural connectivity via computational models using MRI: an extensive comparison study. NeuroImage 111, 65–75. doi: 10.1016/j.neuroimage.2015.02.001 . Murray, J. , 2013. Mathematical biology. Biomathematics. Springer Berlin Heidelberg . Nir, Y. , Massimini, M. , Boly, M. , Tononi, G. , 2013. Sleep and Consciousness. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 133–182 . Nobili, L., Ferrara, M., Moroni, F., Gennaro, L.D., Russo, G.L., Campus, C., Cardinale, F., Carli, F.D., 2011. Dissociated wake-like and sleep-like electro-cortical activity during sleep. NeuroImage 58 (2), 612–619. doi: 10.1016/j.neuroimage.2011.06.032 . Oizumi, M., Albantakis, L., Tononi, G., 2014. From the phenomenology to the mechanisms of consciousness: integrated information theory 3.0. PLoS Comput. Biol. 10 (5), 1–25. doi: 10.1371/journal.pcbi.1003588 . Rao, C. , 2009. Linear statistical inference and its applications, 2nd ed. (With Cd). Wiley Series in Probability and Statistics. Wiley India Pvt. Limited . Robinson, J. , Crighton, D. , Ablowitz, M. , 2001. Infinite-dimensional dynamical systems: an introduction to dissipative parabolic PDEs and the theory of global attractors. Cambridge Texts in Applied Mathematics. Cambridge University Press . Searle, J. , Dennett, D. , Chalmers, D. , 1998. The Mystery of Consciousness. Granta Books . Shanahan, M., 2010. Metastable chimera states in community-structured oscillator networks. Chaos 20 (1), 013108. doi: 10.1063/1.3305451 . Song, C., Boly, M., Tagliazucchi, E., Laufs, H., Tononi, G., 2019. Bold signatures of sleep. bioRxiv doi: 10.1101/531186 . https://www.biorxiv.org/content/early/ 2019/01/26/531186.full.pdf Tagliazucchi, E., Laufs, H., 2014. Decoding wakefulness levels from typical fMRI restingstate data reveals reliable drifts between wakefulness and sleep. Neuron 82 (3), 695– 708. doi: 10.1016/j.neuron.2014.03.020 . Tagliazucchi, E., von Wegner, F., Morzelewski, A., Brodbeck, V., Jahnke, K., Laufs, H., 2013. Breakdown of long-range temporal dependence in default mode and attention networks during deep sleep. Proc. Natl. Acad. Sci. 110 (38), 15419–15424. doi: 10.1073/pnas.1312848110 . https://www.pnas.org/content/ 110/38/15419.full.pdf Takeuchi, Y. , 1996. Global Dynamical Properties of Lotka-Volterra Systems. World Scientific . Takeuchi, Y., Adachi, N., 1980. The existence of globally stable equilibria of ecosystems of the generalized volterra type. J. Math. Biol. 10 (4), 401–415. doi: 10.1007/BF00276098 . Tononi, G. , 2014. How Does Your PHI Formula Deal with the Evidence that Consciousness Is State Dependent? More Specifically, if PHI Were Higher in REM Sleep Than in Waking, Would You Conclude That Dreaming Was More Conscious Even Than Waking?. Springer International Publishing, Cham, pp. 215–217 . Tononi, G., Sporns, O., Edelman, G.M., 1994. A measure for brain complexity: relating functional segregation and integration in the nervous system. Proc. Natl. Acad. Sci. 91 (11), 5033–5037. doi: 10.1073/pnas.91.11.5033 . https://www.pnas.org/content/91/11/5033.full.pdf Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., Mazoyer, B., Joliot, M., 2002. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage 15 (1), 273–289. doi: 10.1006/nimg.2001.0978 . Vohryzek, J., Deco, G., Cessac, B., Kringelbach, M.L., Cabral, J., 2020. Ghost attractors in spontaneous brain activity: recurrent excursions into functionally-relevant bold phase-locking states. Front. Syst. Neurosci. 14, 20. doi: 10.3389/fnsys.2020.00020 . Vyazovskiy, V. , Olcese, U. , Hanlon, E. , Nir, Y. , Cirelli, C. , Tononi, G. , 2011. Local sleep in awake rats. Nature 472, 443–447 . 16 J.A. Galadí, S. Silva Pereira, Y. Sanz Perl et al. NeuroImage 244 (2021) 118551 Yeung, M., Strogatz, S., 1999. Time delay in the Kuramoto model of coupled oscillators. Phys. Rev. Lett. 82 (3), 648–651. doi: 10.1103/PhysRevLett.82.648 . Cited By 405 Yuan, R.-S., Ma, Y.-A., Yuan, B., Ao, P., 2014. Lyapunov function as potential function: a dynamical equivalence. Chin. Phys. B 23 (1), 010505. doi: 10.1088/1674-1056/23/1/010505 . Zamora-López, G. , Chen, Y. , Deco, G. , Kringelbach, M. , Zhou, C. , 2016. Functional complexity emerging from anatomical constraints in the brain: the significance of network modularity and rich-clubs. Sci. Rep. 6, 38424 . Zhou, J.X., Aliyu, M.D.S., Aurell, E., Huang, S., 2012. Quasi-potential landscape in complex multi-stable systems. J. R. Soc. Interface 9 (77), 3539–3553. doi: 10.1098/rsif.2012.0434 . https://royalsocietypublishing.org/doi/pdf/10.1098/ rsif.2012.0434 Further readings Lagzi, F., Atay, F.M., Rotter, S., 2019. Bifurcation analysis of the dynamics of interacting subnetworks of a spiking network. Sci. Rep. 9, 11397. doi: 10.1038/s41598-019-47190-9 . Lagzi, F., Rotter, S., 2015. Dynamics of Competition between Subnetworks of Spiking Neuronal Networks in the Balanced State. PLoS ONE 10(9): e0138947. doi: 10.1371/journal.pone.0138947 . 17