# Entropy production of Multivariate Ornstein-Uhlenbeck processes correlates with consciousness levels in the human brain

Matthieu Gilson, Enzo Tagliazucchi, Rodrigo Cofré

<figure><img src="image_1.png" /></figure>

## To cite this version:

Matthieu Gilson, Enzo Tagliazucchi, Rodrigo Cofré. Entropy production of Multivariate Ornstein-Uhlenbeck processes correlates with consciousness levels in the human brain. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2023, 107, pp.024121. 10.1103/Phys-RevE.107.024121. hal-04441253

HAL Id: hal-04441253

https://amu.hal.science/hal-04441253v1

Submitted on 6 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Entropy production of Multivariate Ornstein-Uhlenbeck processes correlates with consciousness levels in the human brain

Matthieu Gilson¹, Enzo Tagliazucchi², Rodrigo Cofré³,⁴

¹Institut de Neurosciences des Systèmes INSERM-AMU, Marseille, France

²Physics Department University of Buenos Aires and Buenos Aires Physics Institute Argentina

³ CIMFAV-Ingemat, Facultad de Ingeniería, Universidad de Valparaíso, Chile

⁴Institute of Neuroscience (NeuroPSI-CNRS) Paris-Saclay University Gif sur Yvette 91400, France

Consciousness is supported by complex patterns of brain activity which are indicative of irreversible non-equilibrium dynamics. While the framework of stochastic thermodynamics has facilitated the understanding of physical systems of this kind, its application to infer the level of consciousness from empirical data remains elusive. We faced this challenge by calculating entropy production in a multivariate Ornstein-Uhlenbeck process fitted to fMRI brain activity recordings. To test this approach, we focused on the transition from wakefulness to deep sleep, revealing a monotonous relationship between entropy production and the level of consciousness. Our results constitute robust signatures of consciousness while also advancing our understanding of the link between consciousness and complexity from the fundamental perspective of statistical physics.

Animal cognition is the most sophisticated example of information processing found in biological and technological systems [1]. Consciousness, understood as the capacity to sustain subjective experience, can be considered a property that emerges when a sufficiently high level of complex cognitive processing is achieved [2]. From the perspective of physics, consciousness and cognition seem unlikely to emerge from regular and predictable systems, such as those which are in thermodynamic equilibrium and obey the detailed balance equations [3]. Instead, recent research draws a close parallel between the level of consciousness and the entropy production rate of brain activity time series, highlighting temporal irreversibility as a landmark feature of conscious information processing [4–6]. These results suggest a close link between consciousness and non-equilibrium dynamics, prompting a rigorous evaluation from the perspective of stochastic thermodynamics.

In spite of these exciting results, the direct estimation of entropy production from neural activity recordings is undermined by insufficient spatio-temporal sampling, leading to the adoption of heuristics and approximations which lack rigorous justification [3, 4]. To circumvent these limitations, we adopted a framework based on Multivariate Ornstein-Uhlenbeck (MOU) processes, that are widely used for modeling the multivariate dynamics of time series. The importance of MOU derives from the fact that it is the only continuous stationary stochastic process that is simultaneously Gaussian and Markovian. The MOU process is at the heart of many models used to fit fMRI data and to interpret them in terms of whole-brain communication [7–9], in line with the present methodology. We first characterize the non-equilibrium steady state of a generic MOU process. The irreversibility of the process is encoded in the antisymmetric part of the Onsager matrix, while the linearity of the Langevin equations allows us to derive closed-form expression for the entropy production rate in terms of the matrices that define the MOU. As a result, we obtained a model-based

estimation of the entropy production rate for the MOU fitted to fMRI data of subjects transitioning different levels of consciousness during the descent from wakefulness to deep sleep.

## I. MULTIVARIATE ORNSTEIN-UHLENBECK PROCESS

We consider the MOU process closely following the notation in previous work [10]:

$$ \frac{d\boldsymbol{x}(t)}{dt} = -\boldsymbol{B}\boldsymbol{x}(t) + \boldsymbol{\eta}(t). \quad (1) $$

Boldfaced symbols denote vectors and matrices. The inputs $\boldsymbol{\eta}(t)$ correspond to Gaussian white noise with covariance

$$ \langle \boldsymbol{\eta}(t) \boldsymbol{\eta}^{\mathrm{T}}(t') \rangle_t = 2\boldsymbol{D}\delta(t-t'). \quad (2) $$

The angular brackets indicate the mathematical expectation over time and the superscript T the transpose for vectors or matrices. The N-dimensional MOU process is thus defined by two real $N \times N$ matrices, the input covariance matrix $\boldsymbol{D}$, which is symmetric with positive eigenvalues, and the friction matrix $\boldsymbol{B}$, which is not symmetric in general.

### A. Description of the state evolution

Knowing the initial condition $\boldsymbol{x}(0)$ and the realization of the stochastic input $\boldsymbol{\eta}$ over time, the trajectory of the solution of the Eq. (1) is given by:

$$ \boldsymbol{x}(t) = \boldsymbol{G}(t)\boldsymbol{x}(0) + \int_{0}^{t} \boldsymbol{G}(t-s)\boldsymbol{\eta}(s)\mathrm{d}s, \quad (3) $$

where $\boldsymbol{G}(t) = e^{-\boldsymbol{B}t}$ is the Green's function, also known as propagator. In addition to its mean value $\langle \boldsymbol{x}(t) \rangle =$
$G(t)x(0)$, the process is also characterized by its covariance matrix $S(t, t') = \langle x(t)x^T(t') \rangle$. The zero-lag covariance, denoted by $S(t, t)$, obeys the following deterministic differential equation:

$$ \frac{dS(t, t)}{dt} = -B S(t, t) - S(t, t) B^T + 2D. \quad (4) $$

Meanwhile, the lagged covariance with $t' > t$ exhibits an exponential decay as a function of the lag $t' - t$:

$$ S(t, t') = S(t, t) e^{-B^T(t'-t)}. \quad (5) $$

A standard method for analysing Eq. (1), consists in describing the evolution of the probability distribution $P(x, t)$ via the Fokker-Planck equation:

$$ \frac{\partial P(x, t)}{\partial t} = \nabla \cdot [Bx(t) P(x, t) + D \nabla P(x, t)], \quad (6) $$

where $\nabla$ denotes the spatial derivative with respect to $x$. Eq. (6) can be rewritten as a continuity equation of the form

$$ \frac{\partial P(x, t)}{\partial t} + \nabla \cdot \mathbf{J}(x, t) = 0. \quad (7) $$

with the following expression for the probability current (or flux)

$$ \mathbf{J}(x, t) = -\mathbf{D}\nabla P(x, t) - \mathbf{B}x(t)P(x, t) \quad (8) $$

## B. Stationary state and probability current

The Gauss-Markov property of the Ornstein-Uhlenbeck process ensures that the mean and covariances converge exponentially fast toward their respective fixed points, provided the eigenvalues of $B$ (which may be complex) have positive real part. The stationary state of the MOU process exhibits Gaussian fluctuations around a mean equal to zero. This corresponds to the time-independent multivariate probability density

$$ P(x) = \frac{1}{(2\pi)^{N/2}(\det S)^{1/2}} \exp\left(-\frac{1}{2}x^T S^{-1}x\right), \quad (9) $$

where $S$ denotes the fixed point of the zero-lag covariance matrix $S(t, t)$. From Eq. (9), the gradient of $P(x)$ simply reads

$$ \nabla P(x) = \frac{\partial P(x)}{\partial x} = -P(x) S^{-1} x \quad (10) $$

From Eq. (8), the stationary probability current $J(x)$ can thus be rewritten in a compact form

$$ \begin{aligned} \mathbf{J}(x) &= \mathbf{D} P(x) S^{-1} x - \mathbf{B} x P(x) \\ &= \mu x P(x), \end{aligned} \quad (11) $$

with

$$ \mu = D S^{-1} - B \quad (12) $$

## C. Entropy production rate

Going a step further, the (ir)reversibility can be described using thermodynamic variables evaluated for the dynamic process. Using the well-known definition for entropy for the probability distribution $P(x, t)$, now considering its time dependent version, we have

$$ e[P] = - \int_{\mathbb{R}^n} P(x, t) \log P(x, t) \, dx. \quad (13) $$

It can be shown that the rate of the increase of entropy over time can be decomposed into two factors, namely $\dot{e}[P] = EPR - HDR$, where EPR is the entropy production rate and HDR the heat-dissipation rate [10–12]. The EPR is the main quantity of interest here, which we denote by $\Phi$. Now calculating $\Phi$ for the time-independent distribution $P(x)$, we have

$$ \Phi = \int \frac{\mathbf{J}^{\mathrm{T}}(x) \mathbf{D}^{-1} \mathbf{J}(x)}{P(x)} dx = \langle \mathbf{\Pi}^{\mathrm{T}} \mathbf{D} \mathbf{\Pi} \rangle \quad (14) $$

where $\mathbf{\Pi}$ is called the thermodynamic force and is related to $\mathbf{J}$ by the Onsager’s reciprocal relations [11]:

$$ \mathbf{\Pi} = \frac{\mathbf{D}^{-1} \mathbf{J}}{P} \quad (15) $$

The heat-dissipation rate can be computed as follows:

$$ \mathrm{HDR} = \int_{\mathbb{R}^n} \mathbf{D}^{-1} \mathbf{B} x \cdot \mathbf{J} dx \quad (16) $$

In the context of the stationary MOU diffusion processes, a general expression for the entropy production rate per unit time in the stationary state is the following [10, 11, 13]

$$ \Phi = \int (\nabla \log P(x) - D B x)^T D (\nabla \log P(x) - D^{-1} B x) P(x) dx \quad (17) $$

which can be obtained from (12), (14) and (15) as follows:

$$ \begin{aligned} \mu &= \mathbf{D} \mathbf{S}^{-1} - \mathbf{B} \\ \mathbf{D}^{-1} \mu &= \mathbf{S}^{-1} - \mathbf{D}^{-1} \mathbf{B} & \mathbf{D}^{-1}. \\ \mathbf{D}^{-1} \mu x &= (\mathbf{S}^{-1} - \mathbf{D}^{-1} \mathbf{B}) x & \cdot x \\ \mathbf{D}^{-1} \mu x P &= (\mathbf{S}^{-1} - \mathbf{D}^{-1} \mathbf{B}) x P & \cdot P \quad (18) \\ \mathbf{D}^{-1} \mathbf{J} &= (\mathbf{S}^{-1} - \mathbf{D}^{-1} \mathbf{B}) x P & \text{from (11)} \\ \mathbf{\Pi} &= (\mathbf{S}^{-1} - \mathbf{D}^{-1} \mathbf{B}) x & \text{from (15)} \\ \mathbf{\Pi} &= \mathbf{S}^{-1} x - \mathbf{D}^{-1} \mathbf{B} x \end{aligned} $$

Now, as $\nabla \log P(x) = \mathbf{S}^{-1}x$, we obtain (17). From (14)

$$ \langle \mathbf{\Pi}^{\mathrm{T}} \mathbf{D} \mathbf{\Pi} \rangle = \langle x^{\mathrm{T}} (\mathbf{D}^{-1} \mathbf{B} - \mathbf{S}^{-1})^{\mathrm{T}} \mathbf{D} (\mathbf{D}^{-1} \mathbf{B} - \mathbf{S}^{-1}) x \rangle, $$
we obtain that

$$ \Phi = \langle \mathbf{x}^{\mathrm{T}} (\mathbf{D}^{-1}\mathbf{B} - \mathbf{S}^{-1})^{\mathrm{T}} \mathbf{D} (\mathbf{D}^{-1}\mathbf{B} - \mathbf{S}^{-1}) \mathbf{x} \rangle, (19) $$

where the average is taken over the stationary state of the process. From this equation we can verify that when $\mathbf{S} = \mathbf{B}^{-1}\mathbf{D}$, then $\Phi = 0$.

Following previous results [11, 14], a sufficient condition for the MOU process in Eq. (1) to be a time reversible stationary process corresponds to a specific relation between the matrices $\mathbf{B}$ and $\mathbf{D}$:

$$ \mathbf{B}\mathbf{D} = \mathbf{D}\mathbf{B}^{\mathrm{T}}. \qquad (20) $$

To quantify the time (ir)reversibility of the MOU process, it is advantageous to examine the Onsager matrix $\mathbf{L}$ reparameterized using the matrices $\mathbf{B}$, $\mathbf{D}$, and the pairwise zero-lag covariance $\mathbf{S} = \langle \mathbf{x}(t)\mathbf{x}^{\mathrm{T}}(t) \rangle_t$:

$$ \mathbf{L} = \mathbf{B}\mathbf{S} = \mathbf{D} + \mathbf{Q}, \qquad (21) \\ \mathbf{L}^{\mathrm{T}} = \mathbf{S}\mathbf{B}^{\mathrm{T}} = \mathbf{D} - \mathbf{Q}. $$

Here the antisymmetric part $\mathbf{Q}$ of $\mathbf{L}$ provides a measure for the irreversibility of the process. When the process is time reversible $\mathbf{Q} = 0$ and $\mathbf{L}$ is symmetric. The following expression for the entropy production rate $\Phi$ can then be derived from the differential entropy of a multivariate Gaussian, which is a well defined quantity.

From equations (21) and (12), we have $\mathbf{D}^{-1}\mathbf{B} - \mathbf{S}^{-1} = \mathbf{D}^{-1}\mathbf{Q}\mathbf{S}^{-1} = -\mathbf{D}^{-1}\mathbf{\mu}$. Thus, from Eq. (19) considering that $\mathbf{S}$ and $\mathbf{D}$ are symmetric and $\mathbf{Q}$ is anti-symmetric we obtain:

$$ \Phi = - \langle \mathbf{x}^{\mathrm{T}} \mathbf{S}^{-1} \mathbf{Q} \mathbf{D}^{-1} \mathbf{Q} \mathbf{S}^{-1} \mathbf{x} \rangle = \langle \mathbf{x}^{\mathrm{T}} \mathbf{\mu}^{\mathrm{T}} \mathbf{D}^{-1} \mathbf{\mu} \mathbf{x} \rangle \quad (22) $$

The entropy production rate $\Phi$ is non-negative. It is strictly positive if the process is irreversible, and it vanishes only if the process is reversible. Since the stationary state of the MOU is Gaussian with covariance matrix $\mathbf{S}$, we have the following property: $\langle \mathbf{x}^{\mathrm{T}} \mathbf{A} \mathbf{x} \rangle = \mathrm{tr}(\mathbf{S} \mathbf{A})$, and so

$$ \Phi = -\mathrm{tr}(\mathbf{S}^{-1}\mathbf{Q}\mathbf{D}^{-1}\mathbf{Q}) = \mathrm{tr}(\mathbf{S}\mathbf{\mu}^{\mathrm{T}}\mathbf{D}^{-1}\mathbf{\mu}), \quad (23) $$

which can be written into the following equivalent expressions, that does not involve the covariance matrix $\mathbf{S}$ nor its inverse explicitly:

$$ \Phi = \mathrm{tr}(\mathbf{B}^{\mathrm{T}}\mathbf{D}^{-1}\mathbf{Q}) = -\mathrm{tr}(\mathbf{D}^{-1}\mathbf{B}\mathbf{Q}), \quad (24) $$

The entropy production rate $\Phi$ provides a scalar measure for the (ir)reversibility of the whole network process, vanishing only if the process is reversible.

## II. METHODS

### A. Empirical covariance from fMRI data

The model is fitted to reproduce the two covariance matrices calculated from the empirical BOLD signals,

with zero lag and a lag equal to 1 TR:

$$ \hat{S}_{ij}(0) = \frac{1}{T-2} \sum_{1 \le t \le T-1} [x_i(t) - \bar{x}_i][x_j(t) - \bar{x}_j], \quad (25) $$

$$ \hat{S}_{ij}(1) = \frac{1}{T-2} \sum_{1 \le t \le T-1} [x_i(t) - \bar{x}_i][x_j(t+1) - \bar{x}_j](26) $$

Here $\bar{x}_i$ denotes the mean empirical signal: $\bar{x}_i = \frac{1}{T} \sum_t x_i(t)$ for all $i$, which is used to center the data as all variables $x_i$ have mean zero in the model. These are the empirical counterparts of the model covariances $S_{ij}(t, t)$ and $S_{ij}(t, t+1)$ averaged over time $t$.

### B. Parameter estimation of the MOU process

We fit the MOU process from the fMRI time series data for each subject in each sleep condition. We rely on a recent estimation method that tunes the MOU model such that its covariance structure reproduces the matrices in Eq. (25), optimizing its parameters the Jacobian matrix $-\mathbf{B}$ as well as the input covariance matrix $2\mathbf{D}$ [7]. Importantly, this optimization procedure incorporates topological constraints on $\mathbf{B}$, adjusting only existing anatomical connections, also keeping the input cross-covariances $D_{ij} = 0$ for $i \neq j$. Note that our current notation corresponds to a previous publication [7], using the following $-\mathbf{B} \leftrightarrow \mathbf{J}$ and $2\mathbf{D} \leftrightarrow \mathbf{\Sigma}$; note that $-\mathbf{B} \leftrightarrow \mathbf{J}^{\mathrm{T}}$ in the subsequent paper [15].

The model is first calibrated by calculating the time constant $\tau$ from the empirical signals.

$$ \tau = - \frac{N}{\sum_{1 \le i \le N} a(v_i | u)}, \quad (27) $$

where $a(v_i | u)$ is the slope of the linear regression of $v_i = [\log(\hat{S}_{ii}^0), \log(\hat{S}_{ii}^1)]$ by $u = [0, 1]$.

We rely on a gradient descent to iteratively adjust $\mathbf{B}$ and $\mathbf{D}$ until reaching the best fit [7]. At each optimization step, we calculate the model counterparts of the covariance matrices in Eq. (25) $\mathbf{S}(0)$ and $\mathbf{S}(1)$, assuming stationarity over each fMRI session. They can be calculated by solving the Lyapunov equation using e.g. the Bartels-Stewart algorithm, which yields here

$$ \mathbf{B}\mathbf{S}(0) + \mathbf{S}(0)\mathbf{B}^{\mathrm{T}} = 2\mathbf{D}, \quad (28) $$

once again equating the derivative with zero in Eq. (4), and the equation involving the propagator. We calculate the lagged covariance rewriting Eq. (5) for the time-lag equation here as

$$ \mathbf{S}(1) = \mathbf{S}(0) e^{-\mathbf{B}^{\mathrm{T}}}. \quad (29) $$

We then calculate the difference between the model and empirical covariances, $\Delta\mathbf{S}(t) = \hat{\mathbf{S}}(t) - \mathbf{S}(t)$ with $t \in$
{0, 1}. The parameter update is given by differentiating Eqs. (29) and (28):

$$ \Delta \mathbf{B} = \epsilon_B [\mathbf{S}(0)]^{-1} [\Delta \mathbf{S}(0) - \Delta \mathbf{S}(1) e^{\mathbf{B}^\mathrm{T}}] , \quad (30) $$

$$ \Delta \mathbf{D} = \epsilon_D \mathbf{B} \Delta \mathbf{S}(0) + \epsilon_D \Delta \mathbf{S}(0) \mathbf{B}^\mathrm{T}, $$

with $\epsilon_B$ and $\epsilon_D$ small learning rates. The best fit corresponds to minimising the squared norm of both $\Delta \mathbf{S}(0)$ and $\Delta \mathbf{S}(1)$.

### C. MOU-based anatomo-functional model to fit empirical fMRI data

We fitted a MOU process to the time series of blood oxygen level-dependent (BOLD) activity measured using fMRI for a whole-brain parcellation consisting of $N = 90$ regions of interest (ROIs). The BOLD signals were recorded from 15 healthy participants during wakefulness and three sleep stages of progressively deeper unconsciousness (N1, N2, N3). Further details about the data preprocessing like detrending and filtering can be found in [16]. Example BOLD time series are illustrated in Fig. 1A. Fig. 1B-C show two functional connectivity matrices, here calculated as covariances with zero lag $\hat{\mathbf{S}}(0)$ and lag of 1 timestep $\hat{\mathbf{S}}(1)$. These matrices are the empirical counterparts of the model pairwise covariance $\mathbf{S}(l) = \langle \mathbf{x}(t)\mathbf{x}^\mathrm{T}(t+l) \rangle_t$ with lag $l$, which is symmetric for $l = 0$ and was denoted above by $\mathbf{S} = \mathbf{S}(0)$.

In this application, the activity $x_i$ of the MOU process describes the BOLD activity of node $i$. Its friction matrix $\mathbf{B}$ quantifies the propagation of BOLD activity between ROIs, ignoring hemodynamics [18]. Specifically, the diagonal elements $B_{ii}$ are related to a time constant $\tau$ (identical for all ROIs) and the off-diagonal elements $C_{ij} = -B_{ij}$ correspond to the concept of effective connectivity from ROI $j$ to ROI $i$ (excitatory when $C_{ij} > 0$):

$$ -B_{ij} = -\frac{\delta_{ij}}{\tau} + C_{ij}, \quad (31) $$

where $\delta_{ij}$ is the Kronecker delta. The variance $D_{ii}$ reflects the fluctuation amplitude of ROI $i$.

For each subject and condition, the model was fitted to reproduce the two covariance matrices calculated from the empirical BOLD signals $\hat{\mathbf{S}}(0)$ and $\hat{\mathbf{S}}(1)$ (see Fig. 1B-C). We used a recent estimation method based on gradient descent to iteratively adjust $\mathbf{B}$ and $\mathbf{D}$ until reaching the best fit [7]. At each optimization step, we calculate the model counterparts of the covariance matrices $\mathbf{S}(0)$ and $\mathbf{S}(1)$, assuming stationarity over each fMRI session. Importantly, this optimization procedure incorporates topological constraints on $\mathbf{B}$, adjusting only existing anatomical connections (see Fig. 1D-E), also keeping the input cross-covariances $D_{ij} = 0$ for $i \neq j$. Model fit is quantified by two measures: model error, defined using the matrix distance and Pearson correlation between vectorized FC matrices (model versus data). All sleep states have Pearson correlation above 0.6, corresponding to an R2 of 0.36 (See fig S2 in the supplemental material).

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 1. A) Example of the filtered BOLD time series with 198 repetition times (TR) of 2 seconds, corresponding to the 90 ROIs of the AAL parcellation during wakefulness of one participant. B-C) Functional connectivity matrices calculated from the filtered BOLD signals in panel A, $\hat{\mathbf{S}}(0)$ with zero lag and $\hat{\mathbf{S}}(1)$ with a lag of one timestep (TR=2 s). These matrices are used in the objective functions used to fit the anatomo-functional model. D) Generic structural connectivity (SC) obtained from DTI data as described in [17]. E) Mask for existing directional connections to constrain the topology of the $\mathbf{B}$ matrix in the network model (symmetric here).</figcaption></figure>

### D. Robust decoding of sleep stages from MOU parameters

Following previous work [15, 19], we used the scikit-learn Python library for the implementations of multinomial logistic regression (MLR) classifier. The input features corresponded to the vectorized $\mathbf{C}/\mathbf{D}/\mathbf{S}$ matrices after discarding zero or redundant elements. We implemented a stratified cross-validation scheme with 80% of the samples for the train set and 20% for the test set, where the ratio of classes is the same in both sets. We also use the subject identity as “group information” to avoid mixing subject data between the train and test sets. In practice, we use 100 random splits of the data and report the distribution of the accuracies of the 100 splits.

As illustrated in Fig. 2B, both the empirical BOLD variances and the model estimates exhibit global differences across the four sleep stages, although they do not
<figure><img src="image_6.png" /><figcaption>FIG. 2. A) Our dynamic network model has two sets of optimized parameters: the matrix **C** (effective connectivity), which describes the causal interaction between brain regions, and the input variance **D**, which represents the spontaneous activity of each brain region. Note that the topology of the matrix **C** corresponds to the mask inferred from the SC data in Fig. 1D-E, but the weights are estimated from the empirical functional connectivity (FC) matrices Fig. 1B-C, resulting in an anatomo-functional model. B) Changes in total **C** and **D** weights across sleep stages (x-axis), pooled over the 15 subjects. The sleep stages are represented by the blue contrasts, from light for wake (W) to dark for the deepest sleep (N3). C) Classification accuracy based on the model estimates, **C** and **D**, and the empirical covariance matrices. The classifier is the multinomial logistic regression (MLR), which captures changes in individual features across sleep stages. The gray violin plots correspond to the chance-level accuracy calculated empirically by shuffling the labels of the sleep stages.</figcaption></figure>

exhibit a clear trend. These differences in global measures, which are averages over all ROIs, may hide more specific changes at the ROI level, as well as interactions between them. Supplementary Figure S1 shows the good fit of the anatomo-functional model to fMRI data obtained for all sleep stages, with mean correlation between simulated and empirical FC matrices exceeding 0.6 for all stages. Fig. 2C shows that the model estimates give good classification accuracy, both for **C** (in red) and **D** (in purple). This indicates that the model captures the differences in brain dynamics across the sleep stages. Notably, the matrix **C** gives a better classification accuracy than the empirical functional connectivity $\hat{S}(0)$ (in blue), meaning that the model inversion is robust and captures refined information about the sleep stages. Note that the MLR has better accuracy than the 1-nearest-neighbor (1NN) in Suppl Fig S1A, indicating that the changes across sleep stages concern specific features, i.e. connectivity weights (**C**) or nodal spontaneous activity (**D**), rather than their global profile.

<figure><img src="image_7.png" /><figcaption>FIG. 3. A) Violin plots comparing the entropy production rate across sleep stages. Same color coding used in previous plots. The average entropy production values across subjects are for the four sleep stages are 1.99, 1.65, 1.54, and 1.49, respectively. The stars indicate statistical significance for the Mann-Whitney test with *p* < 0.05. B) Comparison of the nodal irreversibility for each ROI (x-axis) between the W and N3 states (in light and dark blue, respectively). The plotted values correspond to the absolute value of sums over rows of **Q**, averaged for homotopic regions; error bars indicate the variability across subjects measured as the standard error of the mean. C) Heatmap plots of the nodal irreversibility on the cortical surface for the W and N3 sleep stages. Note the different color scales for the two stages, for the purpose of better readability. D) To gain insight into the effects of the matrices **C** and **D** on $\Phi$, we shuffle their values as a way to destroy their detailed structures (redistributing their values keeping the topology). We plot $\Phi$ for the estimated **C** and **D** matrices across subjects and sleep stages in blues (same color code used in previous plots), and the shuffled **C** matrices (light gray), **D** matrices (middle gray) and both (dark gray), as a function of the sum of **C** values (left panel), the time constant $\tau$ (middle panel) and the sum of **D** variances (right panel).</figcaption></figure>

E. Reduced entropy production in the transition from wakefulness to deep sleep

Using the condition-specific estimated parameters, we calculated the entropy production rate in the MOU model using Eq. (24). These results in Fig 3A show that entropy production decreases as a function of sleep depth, which in turn implies that dynamics become closer to equilibrium.

The model-based approach allows us to dissect this phenomenon. For all ROIs, we observe that the contri-
bution to Φ, as measured via the nodal irreversibility, defined as $\sum_j |Q_{ij}|$ for each ROI $i$, decreases, as illustrated in Fig 3B. This suggests that the reduction of Φ from W to N3 is a rather global phenomenon, but with a differentiated magnitude across brain regions. Notably, regions in the occipital lobes (cuneus, calcarine, lingual), as well as regions associated to hubs in the default-mode network (precuneus, post cingulate), and the thalamus, remain at a high level of nodal irreversibility in the deep sleep N3; these regions have been shown to exhibit sleep-related changes in previous studies [20–22]. See Suppl Fig S3 for a more detailed comparison across sleep stages.

Last, we examine how the model parameters **C** and **D** contribute to Φ and its reduction across sleep stages. Fig 3D shows a positive relationship between Φ and the sum of weights in **C**, as well as the sum of variances in **D**; conversely, a larger $\tau$ (directly calculated from the empirical BOLD signals) corresponds to a lower Φ. Then we assess the importance of the detailed structures in the **C** and **D** estimates by randomizing them spatially, namely redistributing the total weight/variances across non-zero elements while keeping the same topology and overall sum. We observe the same trends with respect to the **C** and **D** sums, but shifted up or down depending on the surrogates in Fig 3C: randomizing **C** (light gray) decreases slightly Φ, whereas randomizing **D** (middle gray) increases Φ; randomizing both (dark gray) decreases Φ. This indicates that Φ strongly depends on the detailed structures of the **C** and **D** estimates, being larger in the data than in the randomized surrogates. The opposing effects in randomizing **C** and **D** also suggest a balance implemented by the detailed brain dynamics, which results in a controlled level of Φ. Together, our results hint at a positive relationship between the measured Φ and the different levels of consciousness.

### III. DISCUSSION

We measured the entropy production using our anatomo-functional MOU process associated to resting-state fMRI activity recorded from human subjects in different sleep stages. The advantage of our model-based approach is that the entropy production has a closed-form expression from first principles of stochastic thermodynamics for the MOU process, which is numerically fitted to the fMRI data. Our results show high entropy production rate in conscious wakefulness, i.e. correlating positively with the presumed level of cognitive processing. This is consistent with converging theoretical accounts that identify consciousness with an emergent property of a highly complex physical system [2]. These results are also consistent with previous findings relating entropy production with states of consciousness [4–6], with the advantage that do not depend on heuristic approximations. Importantly, our approach allows for identifying the brain regions that contribute most to entropy production. The fulfillment of detailed balance in the

brain is scale-dependent [3]. At the large scale, its violation might relate to the large-scale circuit operations critical for healthy cognition and for the global broadcasting of information which is identified with the computational aspect of consciousness [23]. Because of this, metrics related to the departure from detailed balance (such as entropy production rate) might offer valuable tools to determine levels of consciousness in brain-injured patients and other neurological populations. In summary, assessing temporal irreversibility through entropy production of MOU processes derived from fMRI signals has the potential to highlight different states of consciousness and cognition. More generally, can bridge brain dynamics and thermodynamics, and ultimately help to understand fundamental questions about the brain and consciousness.

**Acknowledgements.** M.G. and R.C. would like to thank the organizers of the 2021 Spring School of the European Institute of Theoretical Neuroscience, the place where this project began to develop. MG was supported by the European Union’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement No. 945539 (Human Brain Project SGA3), as well as the Excellence Initiative of the German federal and state governments (ERS PF-JARA-SDS005). E.T. is supported by grants PICT2018-03103 and PICT-2019-02294 funded by Agencia I+D+I (Argentina), by a Mercator fellowship granted by the German Research Foundation and FONDECYT regular 1220995. R.C was supported by the Human Brain Project, H2020-945539.

### SUPPLEMENTAL MATERIAL

In this supplemental material, we present details on sleep stage decoding sleep stages, supplementary analysis, and a detailed description of the fMRI data.

#### A. Decoding of sleep stages

We use the same approach as in previous work [15, 19]. We rely on two usual classifiers: multinomial logistic regression (MLR) and the 1-nearest-neighbor (1NN). The features correspond to vectorized $C/D/S$ matrices after discarding zero or redundant elements. The MLR is a canonical tool for high-dimensional linear classification, which tunes a weight for each feature, thus selecting the important ones to discriminate the classes. In addition, we use L2-regularization for the MLR ($C = 1.0$ in the scikit-learn implementation). In contrast, the 1NN assigns to a new sample the class to which belongs its closest neighbors with respect to a similarity metric, here chosen as the Pearson correlation coefficient between the feature vectors. It thus relies on the global profile of features to cluster samples into classes.

Following standards, we use a stratified cross-validation scheme with 80% of the samples for the train
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

FIG. S1. A) Similar plot to Fig 2C in the main text with the classification accuracy for the 1-nearest-neighbor (1NN) classifier, which relies on a similarity measure (here the Pearson correlation) between the input features to predict the class of the test sample. The x-axis indicates the features: the model estimates **C** and **D**, as well as the empirical FC denoted by $\hat{\mathcal{S}}$. B) Similar plot to Fig 2B in the main text but for the model input variance summed over all ROIs.

set and 20% for the test set, where the ratio of classes are the same in both sets. We also use the subject identity as “group information” to avoid mixing subject data between the train and test sets. In practice, we use 50 random splits of the data and report the distribution of accuracies of the 50 splits (see the violin plots).

Fig S1 shows that the decoding of the sleep states by the 1NN classifier, which assigns to a new sample the class to which belongs its closest neighbors with respect to a similarity metric, here chosen as the Pearson correlation coefficient between the feature vectors. It thus relies on the global profile of features to cluster samples into classes.

## B. Model fitting and goodness of fit across sleep stages.

The model fit was quantified using two measures. The model error, defined by matrix distance, and the Pearson correlation between vectorized FC matrices (model versus data). As shown in Fig S2A-B, all sleep states have Pearson correlation above 0.6. Note that the changes in the goodness of fit of the model, as measured by the Pearson correlation across sleep stages in Fig S2B, are likely due to the stronger similarity between the corresponding empirical FC matrices in the deep sleep stages than in light sleep stages. In other words, we do not expect this trend to indicate a much better fit of the model for N3 than for W (which might have implications for the calculation of entropy production).

## C. Complementary analysis

The comparison of the nodal irreversibility across sleep states shows an overall decrease for all ROIs from W to N3 (Fig S3A), with a pronounced decrease from W to N1 (Fig S3B) and to a lesser extent from N2 to N3 (Fig S3D), although N1 and N2 are rather similar (Fig S3C). Importantly, we can see heterogeneity in the reduction of

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

FIG. S2. A: Model error for each sleep stage (x-axis) across the 15 subjects. The sleep stages are represented by the blue contrasts, from light for wake (W) to dark for the deepest sleep (N3). B: Goodness of fit as measured by the Pearson correlation between the vectorized model and empirical FC matrices: $S(0)$ with $\hat{S}(0)$, and $S(1)$ with $\hat{S}(1)$. The average of the two Pearson correlation values is reported. C: Isomap decomposition of the empirical zero lag functional connectivity matrices FC0 (or $\hat{S}(0)$) in two dimensions to illustrate their stronger clustering in deep sleep stages N3 (squares) and N2 (triangles), compared to N1 (dots) and W (crosses). D: Pearson similarity across $\hat{S}(0)$ matrices in each sleep stage (x-axis) for all pairs of subjects

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

FIG. S3. A) Plot of the sum of absolute values in *Q* over each row across the W and N1 states. The error bars indicate the s.e.m. across subjects. This plot is another view of the same data in Fig 3B (main text). B-C) Similar plots to panel A for N1 versus N2, N2 versus N3 and W versus N3.

irreversibility across the ROIs in the transition to deep sleep.

Fig S4 shows the correlation between the goodness of fit and the entropy production across all subjects for each brain state. We observe a lack of statistically significant Spearman correlations between variables except for the awake state, where the effect was close to the significance threshold. Fig S2 shows the isomap decomposition of the empirical zero lag functional connectivity matrices FC0 in two dimensions to illustrate their stronger clustering in deep sleep stages compared to N1 and W. We also show the Pearson similarity across FC0 matrices in each sleep stage (to be compared with fig 3B in main text.
<figure><img src="image_16.png" /><figcaption>FIG. S4. Correlation between the goodness of fit and the entropy production across all subjects for each brain state. In each plot, we show the Spearman correlation, denoted *s*, with the corresponding *p*-value. The correlations are not statistically significant, except for the Awake state.</figcaption></figure>

## D. Resting-State fMRI signals

### Participants

A total of 63 healthy subjects (36 females, mean ± SD, 23.4±3.3 years) were selected from a data set previously described in a sleep-related study by Tagliazucchi and Laufs [16]. Participants entered the scanner at 7 PM and were asked to relax, close their eyes, and not fight the sleep onset. A total of 52 minutes of resting state activity were measured with a simultaneous combi-

nation of EEG and fMRI. According to the rules of the American Academy of Sleep Medicine, the polysomnography signals (including the scalp potentials measured with EEG) determine the classification of data into four stages (wakefulness, N1, N2, and N3 sleep). We selected 15 subjects with contiguous resting-state time series of at least 200 volumes to perform our analysis. The local ethics committee approves the experimental protocol (Goethe-Universität Frankfurt, Germany, protocol number: 305/07), and written informed consent was asked to all participants before the experiment. The study was conducted according to the Helsinki Declaration on ethical research.

### MRI data acquisition

MRI images were acquired on a 3-T Siemens Trio scanner (Erlangen, Germany) and fMRI acquisition parameters were 1505 volumes of T2-weighted echo planar images, TR/TE = 2080 ms/30 ms, matrix 64 × 64, voxel size 3×3×3 mm³, distance factor 50%; FOV 192 mm². An optimized polysomnographic setting was employed (chin and tibial EMG, ECG, EOG recorded bipolarly [sampling rate 5 kHz, low pass filter 1 kHz] with 30 EEG channels recorded with FCz as the reference [sampling rate 5 kHz, low pass filter 250 Hz]. Pulse oximetry and respiration were recorded via sensors from the Trio [sampling rate 50 Hz]) and MR scanner-compatible devices (BrainAmp MR+, BrainAmpExG; Brain Products, Gilching, Germany), facilitating sleep scoring during fMRI acquisition.

### Brain parcellation AAL 90 to extract BOLD time series and filtering

To extract the time series of BOLD signals from each participant in a coarse parcellation, we used the AAL90 parcellation with 90 brain areas anatomically defined.BOLD signals (empirical or simulated) were filtered with a Butterworth (order 2) band-pass filter in the 0.01-0.1 Hz frequency range.

[1] G. Piccinini and A. Scarantino, Journal of biological physics **37**, 1 (2011).

[2] A. K. Seth and T. Bayne, Nature Reviews Neuroscience **23**, 7 (2022).

[3] C. W. Lynn, E. J. Cornblath, L. Papadopoulos, M. A. Bertolero, and D. S. Bassett, PNAS **118**, 47 (2021).

[4] Y. Sanz Perl, H. Bocaccio, C. Pallavicini, I. Pérez-Ipiña, S. Laureys, H. Laufs, M. Kringelbach, G. Deco, and E. Tagliazucchi, Physical Review E **104**, 014411 (2021).

[5] L. de la Fuente, F. Zamberlan, H. Bocaccio, M. Kringelbach, G. Deco, Y. S. Perl, and E. Tagliazucchi, Cerebral Cortex **bhac177** (2022).

[6] R. N. Muñoz, A. Leung, A. Zecevik, F. A. Pollock, D. Cohen, B. van Swinderen, N. Tsuchiya, and K. Modi, Phys-

ical Review Research **2**, 023219 (2020).

[7] M. Gilson, R. Moreno-Bote, A. Ponce-Alvarez, P. Ritter, and G. Deco, PLoS Computational Biology **12**, 3 (2016).

[8] K. J. Friston, K. H. Preller, C. Mathys, H. Cagnan, J. Heinzle, A. Razi, and P. Zeidman, NeuroImage **199** (2019).

[9] S. Frässle, Z. M. Manjaly, C. T. Do, L. Kasper, K. P. Pruessmann, and K. E. Stephan, NeuroImage **225** (2021).

[10] C. Godrèche and J. M. Luck, Journal of Physics A: Mathematical and Theoretical **52**, 035002 (2019).

[11] H. Qian, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences **457**, 10.1098/rspa.2001.0811 (2001).
[12] R. Cofré and C. Maldonado, Entropy 20, 10.3390/e20010034 (2018).

[13] H. Qian, M. Qian, and X. Tang, Journal of Statistical Physics (2002).

[14] H. Risken and T. Frank, Springer Series in Synergetics (1996).

[15] M. Gilson, G. Zamora-López, V. Pallarés, M. H. Adhikari, M. Senden, A. T. Campo, D. Mantini, M. Corbetta, G. Deco, and A. Insabato, Network Neuroscience 4, 2 (2020).

[16] E. Tagliazucchi and H. Laufs, Neuron 82, 3 (2014).

[17] I. P. Ipiña, P. D. Kehoe, M. Kringelbach, H. Laufs, A. Ibañez, G. Deco, Y. S. Perl, and E. Tagliazucchi, NeuroImage 215, 116833 (2020).

[18] M. Gilson, R. Moreno-Bote, A. Ponce-Alvarez, P. Ritter, and G. Deco, PLoS Computational Biology 12,

10.1371/journal.pcbi.1004762 (2016).

[19] V. Pallarés, A. Insabato, A. Sanjuán, S. Kühn, D. Mantini, G. Deco, and M. Gilson, NeuroImage 178 (2018).

[20] S. G. Horovitz, A. R. Braun, W. S. Carr, D. Picchioni, T. J. Balkin, M. Fukunaga, and J. H. Duyn, PNAS 106, 27 (2009).

[21] T. T. Dang-Vu, M. Schabus, M. Desseilles, V. Sterpenich, M. Bonjean, and P. Maquet, Sleep 33, 12 (2010).

[22] L. Mirandola, G. Cantalupo, A. E. Vaudano, P. Avanzini, A. Ruggieri, F. Pisani, G. Cossu, C. A. Tassinari, P. F. Nichelli, F. Benuzzi, and S. Meletti, Epilepsy and Behavior Case Reports 106 (2013).

[23] S. Dehaene, M. Kerszberg, and J. P. Changeux, PNAS 95, 24 (1998).