Distinct weak asymmetric interactions shape human brain functions as probability fluxes

TL;DR

This paper reveals that the human brain's functional computations are driven by subtle, task-dependent modifications in asymmetric interactions among brain regions, as evidenced by probability flux patterns and Ising spin models.

Contribution

It introduces a novel mechanism showing that small, asymmetric interaction changes enable brain functions without large energy increases, supported by data analysis and modeling.

Findings

Probability flux patterns vary with task

Symmetric interactions are strong and constant

Antisymmetric interactions are subtle and task-dependent

Abstract

The functional computation of the human brain arises from the collective behaviour of the underlying neural network. The emerging technology enables the recording of population activity in neurons, and the theory of neural networks is expected to explain and extract functional computations from the data. Thermodynamically, a large proportion of the whole-body energy is consumed by the brain, and functional computation of the human brain seems to involve high energy consumption. The human brain, however, does not increase its energy consumption with its function, and most of its energy consumption is not involved in specific brain function: how can the human brain perform its wide repertoire of functional computations without drastically changing its energy consumption? Here, we present a mechanism to perform functional computation by subtle modification of the interaction network among the brain regions. We first show that, by analyzing the data of spontaneous and task-induced whole-cerebral-cortex activity, the probability fluxes, which are the microscopic irreversible measure of state transitions, exhibit unique patterns depending on the task being performed, indicating that the human brain function is a distinct sequence of the brain state transitions. We then fit the parameters of Ising spin systems with asymmetric interactions, where we reveal that the symmetric interactions among the brain regions are strong and task-independent, but the antisymmetric interactions are subtle and task-dependent, and the inferred model reproduces most of the observed probability flux patterns. Our results indicate that the human brain performs its functional computation by subtly modifying the antisymmetric interaction among the brain regions, which might be possible with a small amount of energy.

Equations

$${\mathcal{J}}_{\mu ,\nu }:=w_{\mu ,\nu }{p}_{\nu }-w_{\nu ,\mu }{p}_{\mu },$$

$$\mathbf{J}={\mathbf{J}}^{\text{(s)}}+{\mathbf{J}}^{\text{(a)}},$$

$$w_{\mu ,\nu }{p}_{\nu }\approx \frac{1}{\tau }{n}_{\mu ,\nu },$$

$${\mathcal{J}}_{\mu ,\nu }:=w_{\mu ,\nu }{p}_{\nu }-w_{\nu ,\mu }{p}_{\mu }\approx \frac{1}{\tau }\left(n_{\mu ,\nu }-n_{\nu ,\mu }\right).$$

$$\mathbf{K}:=\left[\begin{array}{cccc}{\mu }_1& {\mu }_2& \cdots & {\mu }_{L-1}\\ {\mu }_2& {\mu }_3& \cdots & {\mu }_L\\ \Delta t_{2,1}& \Delta t_{3,2}& \cdots & \Delta t_{L,L-1}\end{array}\right],$$

$$w_{\mu ,\nu }{p}_{\nu }\approx \frac{1}{{\sum }_{j=1}^{L-1}{K}_{3,j}}\sum _{j=1}^{L-1}{\delta }_{\mu ,K_{2,j}}{\delta }_{\nu ,K_{1,j}},$$

$$p_{\nu }\approx \frac{1}{{\sum }_{j=1}^{L-1}{K}_{3,j}}\sum _{j=1}^{L-1}{K}_{3,j}{\delta }_{\nu ,K_{1,j}}.$$

$${\left\{t_j^{\prime }\right\}}_{j=1}^{M^{\prime }}=\underset{t}{\mathrm{argwhere}}\left[\frac{{\mathrm{d}}^{\alpha }}{\mathrm{d}{t}^{\alpha }}{f}_i\left(t\right)=0,\forall i\right],\alpha \in \left\{0,1,2\right\},$$

$$s_i\left(t\right)=\{\begin{array}{ll}+1& \text{if }\frac{{\mathrm{d}}^0}{\mathrm{d}{t}^0}{f}_i\left(t\right)>0\\ -1& \text{if }\frac{{\mathrm{d}}^0}{\mathrm{d}{t}^0}{f}_i\left(t\right)<0\\ -s_i\left(t-0\right)& \text{if }\frac{{\mathrm{d}}^0}{\mathrm{d}{t}^0}{f}_i\left(t\right)=0\end{array},$$

$$s_i\left(t\right)=\{\begin{array}{ll}+1& \text{if }\frac{{\mathrm{d}}^1}{\mathrm{d}{t}^1}{f}_i\left(t\right)>0\\ -1& \text{if }\frac{{\mathrm{d}}^1}{\mathrm{d}{t}^1}{f}_i\left(t\right)<0\\ -s_i\left(t-0\right)& \text{if }\frac{{\mathrm{d}}^1}{\mathrm{d}{t}^1}{f}_i\left(t\right)=0\end{array},$$

$$s_i\left(t\right)=\{\begin{array}{ll}+1& \text{if }\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}{f}_i\left(t\right)<0\\ -1& \text{if }\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}{f}_i\left(t\right)>0\\ -s_i\left(t-0\right)& \text{if }\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}{f}_i\left(t\right)=0\end{array}.$$

$$\mathbf{K}=\left[\begin{array}{ccccc}\mu \left(0\right)& \mu \left(t_1^{\prime }\right)& \cdots & \mu \left(t_{M^{\prime }-1}^{\prime }\right)& \mu \left(t_{M^{\prime }}^{\prime }\right)\\ \mu \left(t_1^{\prime }\right)& \mu \left(t_2^{\prime }\right)& \cdots & \mu \left(t_{M^{\prime }}^{\prime }\right)& \mu \left(\tau \right)\\ \Delta t_{1,0}^{\prime }& \Delta t_{2,1}^{\prime }& \cdots & \Delta t_{M^{\prime },M^{\prime }-1}^{\prime }& \tau -t_{M^{\prime }}^{\prime }\end{array}\right],$$

$$\Delta p_{\mu }:=\sum _{\nu =1}^{2^N}{\mathcal{J}}_{\mu ,\nu }.$$

$$\Sigma :=⟨\left(\mathbf{s}-⟨\mathbf{s}⟩\right){\left(\mathbf{s}-⟨\mathbf{s}⟩\right)}^{\top }⟩,$$

$${\tilde{\mathbf{e}}}_i:=2\left[\begin{array}{c}{v}_{i;j}\\ v_{i;k}\end{array}\right],$$

$$\mathcal{H}\left(\mathbf{s}\right):=-\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{s}_i{J}_{i,j}{s}_j-\sum _{i=1}^N{s}_i{h}_i=-\frac{1}{2}{\mathbf{s}}^{\top }\mathbf{J}\mathbf{s}-{\mathbf{s}}^{\top }\mathbf{h},$$

$$w_{\mu ,\nu }\approx \frac{1}{\tau }\frac{n_{\mu ,\nu }}{p_{\nu }}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{p}_{\mu }\left(t\right)$$

$$=\sum _{\nu =1}^{2^N}{\mathcal{J}}_{\mu ,\nu }\left(t\right),$$

$$w_{\mu ,\nu }{p}_{\nu }=w_{\nu ,\mu }{p}_{\mu },\forall \left(\mu ,\nu \right),$$

$$\frac{w_{\mu ,\nu }}{w_{\nu ,\mu }}=\frac{p_{\mu }}{p_{\nu }}.$$

$${\pi }_{\mu }=\frac{1}{\mathcal{Z}}\exp \left[-\beta \mathcal{H}\left({\mathbf{s}}_{\mu }\right)\right],$$

$$\frac{w_{\mu ,\nu }}{w_{\nu ,\mu }}=\frac{\exp \left(-\beta \frac{\Delta E_{\mu ,\nu }}{2}\right)}{\exp \left(-\beta \frac{\Delta E_{\nu ,\mu }}{2}\right)}=\frac{\frac{1}{1+\exp \left(\beta \Delta E_{\mu ,\nu }\right)}}{\frac{1}{1+\exp \left(\beta \Delta E_{\nu ,\mu }\right)}},$$

$$w_{\mu ,\nu }=A\exp \left(-\beta \frac{\Delta E_{\mu ,\nu }}{2}\right),$$

$$w_{\mu ,\nu }=A\frac{1}{1+\exp \left(\beta \Delta E_{\mu ,\nu }\right)},$$

$$\Delta E_{{\mu }^{\left(k\right)},\mu }$$

$$=-{\left({\mathbf{s}}_{{\mu }^{\left(k\right)}}-{\mathbf{s}}_{\mu }\right)}^{\top }\left(\mathbf{J}{\mathbf{s}}_{\mu }+\mathbf{h}\right),$$

$$w_{{\mu }^{\left(k\right)},\mu }$$

$$=A\exp \left[\beta \frac{1}{2}{\left({\mathbf{s}}_{{\mu }^{\left(k\right)}}-{\mathbf{s}}_{\mu }\right)}^{\top }\left(\mathbf{J}{\mathbf{s}}_{\mu }+\mathbf{h}\right)\right].$$

$$L\left(\beta \mathbf{J},\beta \mathbf{h},A\right)$$

$$:=$$

$$\displaystyle\bigg{.}~-2\ln\left(\frac{1}{\tau}\frac{n_{\mu^{\left(k\right)},\mu}}{p_{\mu}}\right)\bigg{]}^{2},$$

$$\beta \tilde{\mathbf{J}},\beta \tilde{\mathbf{h}},\tilde{A}=\underset{\left\{\beta \mathbf{J},\beta \mathbf{h},A\right\}}{\mathrm{argmin}}\left[L\left(\beta \mathbf{J},\beta \mathbf{h},A\right)\right].$$

Full text

Introduction

How the human brain works is mysterious even with the accumulation of detailed knowledge [28]. Computation seems to be a possible and probable analogy to explain the function of the human brain [29]. The neuron, the fundamental building block of the brain, exhibits the binary spiking activity, and such observation has led to assume the human brain as the computer. Except for the digital nature of the computation, there are various differences between the human brain and the computer, but their aims are the same: information processing.

Following the analogy of the logic circuit of computers, the population of neurons is modelled as a network of neurons, i.e., a neural network [30]. Interacting neurons exhibit the emergent collective behaviour [3], which is more than the sum of the individual neurons [31]. Such collective behaviour arising from the underlying interaction corresponds to the property of the system. For biological system like neural networks, it is the functional property [1]. Among the emergent collective behaviour, the state transition dynamics describes the computation [2]. Emerging function as state transition dynamics arise from the underlying network structure of the interacting neurons [1] and revealing and explaining them is the ultimate aim of neuroscience [4].

Over the decades, the technological advancements enable one to record the population activity of neurons [5, 6], and the brain research has been transformed into big science [32]. The collected big data has opened “the new century of the brain” [33] but such large, complex, and high-dimensional data is difficult to analyze, and explanatory rather than descriptive model is demanded [7]. The theory of neural network, rather than that of a single neuron, is expected to explain and extract emergent function of the human brain from data [8, 34, 35, 9]. Surge of data allow data-driven quantitative approach to human brain through the lens of network [10], and physics-rooted intuition reveals insights not only into the structure and dynamics of the human brain but also into controlling them [11].

In addition to digital circuit nature, the human brain has another characteristic aspect: its tremendous energy consumption rate. The human brain accounts only for 2% of the body weight, but it consumes 20% of the total metabolic rate at rest [12, 13, 14], and the consumption rate reaches even 66% in childhood [15]. That of other vertebrates account for $\lesssim$10% of the basal metabolic rate [12, 13], thus the energy consumption by human brain consume more than twice of that of the other vertebrates, including primates [12, 13]. Furthermore, among human organs, energy consumption rate of the brain is the highest [14], and it is the third-highest in energy consumption per weight after the heart and kidney [14]. Such intense energy consumption is due to the cerebral cortex, which accounts for the majority ($\gtrsim$60%) of energy consumption in the brain [13]. The living systems are alive by keeping themselves out of equilibrium [36, 37] through the constant consumption of energy [37]. The energy consumption of the human brain drives itself to exhibit nonequilibrium dynamics to perform its functions [16, 17].

Although the high energy consumption rate is the characteristic trait of the human brain, its relation with the functional computation is vague. Indeed, the brain increases its energy consumption rate as demanded, but the increase is as small fraction ($\sim$1%) of total energy consumption rate [19, 20, 38]. Furthermore, most (60 to 80%) of the energy consumption in human brain does not involve in specific brain function [19, 20, 38]. Such function-unrelated energy consumption is referred to as “dark energy” [19, 38, 39]. The dark energy of human brain arise from interneuron communication [18, 19, 20, 22, 21] rather than the computation [22, 21]. Thus, the human brain consumes less energy to perform its computational function and a large portion of energy is used for task-unrelated interneuron communications. Then, how can the human brain perform its functions without drastically increasing its energy consumption rate from its baseline? Moreover, how can human brain operation switches from one function being performed to another—among the wide repertoires of functions—only with the subtle change of energy consumption rate? To begin with, why does the human brain consume a large portion of energy for network structure?

Here, we reveal that nonequilibrium state transition dynamics of the human brain exhibits the unique sequential patterns depending on tasks and such patterns emerge by a subtle change in the asymmetric part of the interaction network among the brain regions. Analysing the whole-cerebral-cortex activity data recorded through blood-oxygen-level-dependent (BOLD) functional magnetic resonance imaging (fMRI), we show that the state transition dynamics of the human brain exhibit unique pattern depending on tasks. We then show that the task-dependent pattern of dynamics arise from the spatio-temporal pattern of human brain activity—the brain regions collectively change its activity with time according to task. To understand the structural origin of these task-dependent dynamics, we investigate the underlying interaction network among the brain regions by developing a method to fit the Ising spin system with the asymmetric interaction—the prototypical model to study nonequilibrium collective behaviour on networks—to the data. We find that the symmetric part of the interaction matrix is similar across tasks, but the asymmetric part of the interaction matrix is not. Finally, we confirm that our model captures the task-dependent dynamics observed in the data by comparing the predicted dynamics from the Ising spin system with the empirical data.

Task-dependent irreversible dynamics of human brains

Firstly, we develop the procedure to analyze the irreversible dynamics of the human brain. Our targeted data (BOLD fMRI data in Methods) contain the whole-cerebral-cortex (divided into 100 cortical parcels [40]) activity recorded through BOLD fMRI as a part of the Human Connectome Project [23]. It consists of BOLD fMRI signal from 590 healthy adults at rest and during seven cognitive and motor tasks. Each of them are time series and approximately first three minutes are analyzed.

To examine the task-dependent irreversible dynamics of human brains, we perform the hypercubic probability flux analysis (Hypercubic probability flux analysis in Methods). The probability flux is a measure of the broken detailed balance—the condition of reversibility—or arrow of time. The probability flux characterize the nonequilibrium state transition [41, 42] and probability flux analysis reveal the non-trivial probability flux from data [43]. Formally, the probability flux from state $\nu$ to $\mu$, $\mathcal{J}_{\mu,\nu}$, is defined as the difference between the forward (from state $\nu$ to $\mu$) and backward (from state $\mu$ to $\nu$) joint transition rate,

[TABLE]

where $w_{\mu,\nu}$ is the transition rate from state $\nu$ to $\mu$, and $p_{\mu}$ is the probability of finding the system in state $\mu$. To perform the probability flux analysis (Probability flux analysis in Methods) from the empirical data, we first reduce the spatial-dimensionality of the data (Spatial coarse-graining of brain region through hierarchical clustering in Methods), then we temporally binarize the data (Temporal coarse graining through time series binarization in Methods).

As the first step of probability flux analysis [43] (Probability flux analysis in Methods), we define the discrete state space of data. The challenge of applying the probability flux analysis is that the high-dimensionality of the human brain dynamics, which necessitates the use of dimensionality reduction. In the spirit of the idea of renormalization group [44] of statistical physics, we seek the coarse-grained representation of the data while preserving the qualitative feature (Spatial coarse-graining of brain region through hierarchical clustering in Methods). We coarse-grain the brain regions following the correlation among them [45]. Based on the correlation matrix, we perform the hierarchical clustering [46], which identifies the hierarchical structure of the correlation among the brain regions (Fig. 1a). We then define the seven clusters of brain regions that exhibit similar dynamics by manually deciding the threshold for the dendrogram (Fig. 1b). The resulting cluster of brain regions is shown in Fig. 1c. To understand the functional meaning of the clusters, we compare that how our cluster structure aligns with the known functional clusters [47] (Fig. 1d). We find roughly seven clusters of brain regions that correspond to the known seven functional clusters with minor non-alignment (Fig. 1e). Our correlation based coarse-graining seems to capture the functional feature of human brain. As the neural activity within each cluster is highly correlated, we average the neural activity within each cluster to obtain a time series of coarse-grained activity for each cluster.

As the second step, we binarize the seven-dimensional time series into a series of seven-dimensional symbols (Temporal coarse graining through time series binarization in Methods). By applying a threshold to the time series itself or its differentiation with time, we obtain the binarized time series, where the brain region is either active ($+1$) or inactive ($-1$) at any time point. With this binarized representation, we can analyze the probability flux between the states following the procedure of probability flux analysis [43], where the probability flux is estimated as the number of observed transitions between states. The stationarity of the probability distribution is examined and most of the states are stationary (Examining the assumption of stationary distribution in Methods and Extended Data Fig. 1).

We then visually examine the probability flux of the human brain by projecting its state space onto the two-dimensional plane (Fig. 2 and Hypercubic probability flux diagram in Methods). Because the state space of binarized brain activity is high-dimensional cube or hypercube, the state transition corresponds to the edge of hypercube. Thus, we can visualize the probability flux as a hypercubic edge. We employ principal component analysis (PCA) [48] to project hypercube onto a two-dimensional plane [49]. Visualizing the probability flux allow one to identify the feature of nonequilibrium dynamics [43, 16, 17], and the such visualization of probability fluxes provides the theoretical foundation of nonequilibrium steady states [50, 51, 41, 42].

At first glance, the resulting probability flux diagrams (Fig. 2) reveal distinct patterns of brain activity depending on the task being performed. They exhibit cycles of probability flux, which is one of the characteristic feature of nonequilibrium steady state [50, 51, 41, 42]. The difference among the tasks is the size, strength, and number of cycles of probability flux: some tasks exhibit cycles involving more states than the others, the magnitude of the probability flux is different among the tasks, and the number of such structures varies depending on the tasks. As mentioned in ref. [16], our results indicates that the task-dependent unique pattern of probability flux is closely related to the cognitive processes involved in each task. Furthermore, those hypercubic probability flux diagrams also indicate that the probability flux or biased state transition pattern arise from the task-dependent sequential pattern of brain activity—the correlated activity among different brain regions depends on the task being performed. For example, as parts of state transition pathways—the chain of probability flux, the probability flux involving cluster 2 in Fig. 1b and 1c (related to frontparietal) exhibits strong magnitude in the social and language tasks, while the probability flux involving cluster 1 in Fig. 1b and 1c (related to default mode) exhibits strong magnitude in the motor task. The state transition pathways themselves are also very distinct. We conclude that the sequence of brain state transitions or order of activation (or inactivation) of brain regions is distinct across the task being performed—which suggest the such pattern represents the human brain functions. We note that the unique patterns are seen by different method to show the probability flux such as PCA projected state space (Extended Data Fig. 2). We confirm that clustering defined from known functional clusters [47] (Fig. 1d) does not change the results (Extended Data Fig. 3), number of clusters does not change the results (Extended Data Fig. 4), discretization method does not change the results (Extended Data Fig. 5). The visualization of the hypercubic probability flux diagram using other PCs are available in Extended Data Fig. 6.

Human brain function and the probability flux

We then compare the task-specific pattern of probability flux to investigate the relation among tasks. In Fig. 3a and 3b, we perform hierarchical clustering of the tasks as probability flux to reveal the similarity among tasks. We find that there is large correlated group consisting of rest and five tasks (social, gambling, relational, working memory, and motor) and two tasks (emotion and language) are not strongly correlated or negatively correlated with others. This result indicates that the emotion and language tasks may involve distinct information processing compared to the other tasks.

To characterize each state transition or directed hypercubic edge, we define the eight-dimensional vector, where each element corresponds to the probability flux in the task. We perform hierarchical clustering of the directed hypercubic edge as the eight-dimensional vector. In Fig. 3c, we show the results of the hierarchical clustering of the cosine similarity matrix of directed hypercubic edges. There are several correlated clusters in Fig. 3c, indicating that there are directed hypercubic edges sharing the probability flux pattern—which might be the fundamental state transition of the human brain.

To validate our findings, we visualize, in Fig. 3d the probability fluxes of each tasks as list sorted by the results of the hierarchical clustering of Fig. 3a and 3c. By comparing the columns of the matrix in Fig. 3d, we find that the probability fluxes exhibit patterns characteristic of the tasks as visualized in Fig. 2, but at the same time, there are some minor similarity between them. The comparison of the rows in the matrix of Fig. 3d reveal that there are indeed some group of probability flux which is similar across different tasks, as Fig. 3d indicates. Together, these results suggest that there are shared dynamics in the probability fluxes across tasks, although the major features of the probability fluxes are still task-dependent—which might reflect the separation of fundamental task-independent dynamics and specific task-dependent dynamics.

The structural origin of task-dependent irreversible dynamics

So far we reveal that the human brain exhibit distinct sequence of brain state depending on the task being performed. To reveal the mechanism of such task-dependent dynamics, or task-dependent emergent collective behaviour, we investigate the structural origin of the observed probability flux. To infer the underlying interaction network structure exhibiting nonequilibrium dynamics, we present a method (Infering Ising spin system from transition rates in Methods) to infer the (asymmetric) Ising spin system from the transition rates $\left\{w_{\mu,\nu}\right\}$. Because there are probability flux in the hypercubic state space of human brain, equilibrium based methods are incapable to capture the nonequilibrium feature of the system. Assuming the stationarity of probability distribution and asymmetric (nonconservative) pairwise interaction between the brain clusters, we build a stochastic model of state transition based on a (pseudo-)Hamiltonian. With our method, we can reconstruct the interaction network and external input among the brain clusters.

To examine the possible structure exhibiting the emerging probability flux, in Fig. 4, we infer the interaction networks $\left\{\bm{J}\right\}$ and external inputs $\left\{\bm{h}\right\}$ from the estimated transition rates of each task. We assume the Arrhenius type transition rate (Infering Ising spin system from transition rates in Methods), The inferred interaction matrix (Fig. 4a) shows the interaction between the brain clusters defined in Fig. 1b and 1c. The element $J_{i,j}$ represents the interaction strength from cluster $j$ to cluster $i$, and its sign is the type of interaction (positive means excitatory and negative means inhibitory). It seems there is no significant difference among the inferred interaction across tasks. Majority of the interaction is negative, indicating that the brain clusters tend to inhibit each other’s activity, which is consistent with the results of the hierarchical clustering of the original unbinarized data (Fig. 1a).

The task-dependent difference of the inferred interactions are revealed by the decomposition of the interaction matrix. In general, the interaction matrix $\bm{J}$ is decomposed into the symmetric and antisymmetric parts:

[TABLE]

where $\bm{J}^{\text{(s)}}\coloneqq\frac{1}{2}\left(\bm{J}+\bm{J}^{\top}\right)$ is the symmetric part of the interaction matrix, and $\bm{J}^{\text{(a)}}\coloneqq\frac{1}{2}\left(\bm{J}-\bm{J}^{\top}\right)$ is the antisymmetric part of the interaction matrix. Note that the element of the interaction matrix is generally asymmetric $J_{i,j}\neq J_{j,i}$, but that of the symmetric part is $J_{i,j}=J_{j,i}$, and that of the antisymmetric part is $J_{i,j}=-J_{j,i}$. In Figs. 4b and 4c, we show the symmetric and antisymmetric parts of the inferred interaction matrix, respectively. The symmetric part does not exhibit significant difference across tasks, while the antisymmetric part shows more variability across tasks. The asymmetric part of the interaction matrix is the source of the nonequilibrium steady state [16], and indeed the inferred model exhibits such interaction depending on the tasks. This indicates that the antisymmetric interactions are modified by the tasks being performed. Moreover, the magnitude of the elements of antisymmetric part seems negligibly smaller than the that of the symmetric part, indicating that the antisymmetric interactions may be easily modified for adaptation to the specific task demands. This suggests the answer to the question we raise in the beginning: the brain is expected to achieve efficient information processing by slightly modifying the asymmetric part of the interaction matrix based on the task being performed.

The interaction network structure of the brain clusters are shown in Fig. 4d to examine the structural features. As indicated from the interaction matrix (Fig. 4a), the majority of the interactions are negative, which leads the system to have the geometrical frustration [52]. Although our inferred interaction matrix is not strictly symmetric and ground states are ill-defined, the cycle structure of the probability fluxes (Fig. 2) seems to involve those ill-defined ground states. There are particularly strong negative interactions involving the cluster 1 (related to default mode): the brain region related to the default mode inhibits the activity of other brain regions, which is consistent with the known feature of default mode [20].

Moving onto the external input (Fig. 4e and 4f), we find the external input exhibits the large difference depending on the tasks. The task-dependent external stimulus from the out of cortex may be reflected in the inferred external input $\bm{h}$. Nevertheless, there is a shared feature of external input among the majority of tasks: the strong negative field to the brain cluster 1 (related to default mode). Considering that the negative interaction involving the cluster 1 as shown in Fig. 4a and 4d, this negative external input to cluster 1 indicates that the brain region related to default mode is suppressed by the external stimulus, and it leads to the activation of other brain regions to perform information processing. We note that transition rate of the Glauber type [53] can also infer the probability flux, but the rate constant largely differ among tasks (Extended Data Fig. 7 and Extended Data Fig. 8).

To validate our inferred model, we reconstruct the probability flux and compare it with the original data. We show, Fig. 5, the reconstructed probability flux. By comparing the original probability fluxes with the reconstructed ones, we find the quantitative agreement for most tasks, except for the working memory task. This is confirmed by calculating the correlation and indeed the two sets of probability fluxes are correlated (Pearson coefficient larger than 0.7) except for the working memory task. By introducing more complexity to the model, this may be improved, but we believe the overall feature of the probability flux is captured by the method we present in this study. Our result is not from the variability of the external input because if we fix the interaction network among tasks, the probability fluxes are not well reconstructed (Extended Data Figs. 9 and 10).

Conclusion

F. H. C. Crick once mentioned that when one think about the brain, “We sense there is something difficult to explain, but it seems almost impossible to state clearly and exactly what the difficulty is” [28], and theoretical models, particularly that of neural networks rather than the single neuron, have been the candidate to explain the human brain. The emergent functional computation of the neural network, at first glance, appears to consume energy, but the human brain operating in function does not significantly increase its energy consumption. It is unclear how the human brain performs its wide variety of functions without drastically increasing its energy consumption. In this study, through the analysis of fMRI data, we suggest that the human brain function is characterized by the sequential pattern of state transitions, and such pattern emerge from the asymmetric part of the interaction network between the brain regions. Rather than changing the energy consumption to perform its function, the human brain might slightly change its underlying interaction network to exhibit varying nonequilibrium state transition dynamics. The results indicate that the human brain function is defined by dynamic sequential pattern rather the static single pattern, which in common with the idea of wave-like motifs [54]. Our finding may provide a different view on the cognition as sequential and dynamic rather than single static computation and representation [55]. As a problem of physics, we show that the probability flux (how the system breaks the detailed balance) rather than the entropy production rate (which quantifies the breaking of detailed balance in the all state space, Entropy production rate in Methods) has a richer organization and information of the nonequilibrium dynamics. Perhaps the diagrammatic description of nonequilibrium steady state provide not only the alternative technique to calculate entropy production rate [50, 51] but also the unexplored aspect of nonequilibrium statistical physics, particularly that of stochastic thermodynamics [56, 57, 58, 59, 60]. On the methodological side of this work, unlike the method such as maximum entropy modelling [61] or energy landscape analysis [62], we present a conceptually different approach to analyze the time series data by focusing on the probability flux—the characteristic of the nonequilibrium steady state [41, 42]—rather than the static correlation. The method presented here is not restricted to analyze neuroscientific data and can be applied to reveal nonequilibrium aspects of other high-dimensional or many-body systems.

Although our study shed light on the mechanism of the human brain function from the perspective of nonequilibrium state transition dynamics, interaction network structure and energy consumption, there are several assumptions and limitations. The first assumption is that the coarse-graining of the brain regions does not change the quantitative properties of the data. We define the seven clusters of parcels based on the hierarchical clustering to reduce the dimensionality of the data. Practically this is for reducing the computational cost and finite data availability but biologically each of seven clusters is related to the specific functions of the brain [47]. As we show in Fig. 1e, our clusters is roughly aligned with the known functional clusters [47]. The second assumption is that the binarization captures the essential features of the data. Considering the binary nature of underlying neurons of human brain, grasping the trend of the activity by binarization seems reasonable, but further examination is needed to validate what features are not captured. The third assumption is that the probability distribution is independent of time, i.e., the stationary distribution. The timescale of the neuronal spiking and BOLD fMRI signal is several orders smaller than the duration of scanning, and we fairly assume that the underlying dynamics to be approximately stationary over the scanning period (Examining the assumption of stationary distribution in Methods and Extended Data Fig. 1). Assessing the non-stationarity of data is one of the future direction. Our fourth assumption, the Markov process, seems to be the reasonable first step to approach. Indeed, we successfully reconstruct the pattern of probability fluxes from the inferred Ising spin system, but the reconstruction is not perfect and in task of working memory it fails to capture the feature of dynamics (Fig. 5). Our method may be improved by introducing the higher-order Markov process, which reduces the error of the probability flux analysis [43]. Fifthly, the inferred Ising spin system is based on the assumption of pairwise interactions between brain regions. As we show in Fig. 5, the model struggles to reconstruct the probability flux of a task. This may be due to the simplification of the underlying network structure. It may be improved by considering the higher-order interactions beyond pairwise interactions [61, 63, 64]. Finally, the biological origin of the external input term of the inferred model is not clear. We expect some sensory or perceptional signals to contribute to the external input, but this needs further investigation. Nevertheless, as shown in Figs. 4e and 4f, the external input suppresses the brain regions which are active during the default mode, which is consistent with the known feature. Conquering those limitations might improve the present method and may even reveal a new aspect of human brain dynamics and underlying network structure.

From the engineering perspective, our work shows a possible mechanism of a new computational technology. In real, any computation is performed with spatio-temporal and energetic constraints and such constraints drive the system into far from equilibrium to perform the computation [26]. Although the human brain intensely consume the energy compared to the other vertebrates [12, 13], the energy consumption rate in the human brain is \sim$$10^{5} times more efficient than the digital simulation of them [65]. The constrained computation like human brain dynamics may provide insight into an energy efficient and, at the same time, advanced functional computing [27] utilizing the nonequilibrium dynamics. The neuromorphic system co-locating memory and processor—without von Neumann architecture [27]—-may be developed based on the nonequilibrium nature of the human brain. From the algorithmic or architectural side of the software, the neuroscience has been inspired the artificial intelligence research [25], and our study may lead to probability flux based models. The reverse is also true: to make the artificial intelligence more human [66], we need to obtain the deep understanding of “how matter becomes mind” [67]. Our work provides a new framework for understanding the mind as emergent patterns of probability flux.

Methods

BOLD fMRI data

We analyze the previously collected [24, 23] and preprocessed [16] BOLD fMRI data. The data consists collected from 590 healthy adults and for each participant, the recordings were performed during seven cognitive and motor tasks (rest, emotion, working memory, social, language, gambling, and motor). The cortex is parcelled into 100 regions following the previously published work [40]. The BOLD fMRI data contains different phase encoding directions, specifically left-to-right (LR) and right-to-left (RL). The data duration for each task is different, so we analyze the first 176 time points for each task following the shortest task duration (emotional task) to avoid the systematic difference due to varying task lengths [16].

Hypercubic probability flux analysis

By combining the Probability flux analysis, Spatial coarse-graining of brain region through hierarchical clustering, Temporal coarse graining through time series binarization, and Hypercubic probability flux diagram of this Methods, we establish the hypercubic probability flux analysis. By considering the dynamic correlation rather than the static correlation, our analysis reveals the nonequilibrium dynamics and underlying asymmetric interaction networks of data.

Probability flux analysis

Under the assumption of the stationary dynamics, the probability flux can be estimated from the time series of state (phase) of the target system [43]. Below is the procedure of the probability flux analysis [43]. First, the discrete state of the system is defined by coarse-graining, resulting in coarse-grained state space. Thus, the original trajectory data is converted to the discrete state time series. Next, the joint transition rate is estimated by assuming steady state, i.e., $\frac{\mathrm{d}}{\mathrm{d}t}p_{\mu}\left(t\right)=0$ for all $\mu$. The probability of finding system in state $\mu$ at time $t$ is denoted as $p_{\mu}\left(t\right)$, which satisfies the normalization $\sum_{\mu}p_{\mu}\left(t\right)=1$. The joint transition rate $w_{\mu,\nu}p_{\nu}\left(t\right)$ from state $\nu$ to $\mu$ at time $t$ is the multiplication of the probability of existence $p_{\nu}\left(t\right)\in\left[0,1\right]$ by transition rate from state $\nu$ to $\mu$, $w_{\mu,\nu}\in\mathbb{R}_{\geq 0}$. Thus, in steady state, the joint transition rate corresponds to the number of transition from one state to the other during the unit time, i.e.,

[TABLE]

where $n_{\mu,\nu}$ is the number of transition from state $\nu$ to $\mu$, and $\tau$ is the time duration of the observation. In coarse-grained state space, the original data is discretized into the coarse-grained state, and the number of transition from one state to the other can be counted. Finally, the probability flux in steady state is estimated as the difference between the forward joint transition rate and backward joint transition rate (equation (1)), i.e.,

[TABLE]

For finite observation time $\tau$, the probability flux can be affected by the finite sampling effects, and thus, it is vital to consider the statistical significance of the estimated flux values. The bootstrapping [68] is performed to assess the variability of the estimated flux values. Assuming Markov process, we perform the trajectory bootstrapping [43], which resamples the observed state transition. Consider the coarse-grained state time series with $L$ time points, $\left\{\mu_{1},\mu_{2},\ldots,\mu_{L}\right\}$. The state transition is recorded as a state transition matrix of size $3\times\left(L-1\right)$,

[TABLE]

where $\mu_{i}$ is the coarse-grained state at $i$th time point, and $\Delta t_{i,j}\coloneqq t_{i}-t_{j}$ is the time difference between the $i$th and $j$th time points, i.e., the time staying in the state $\mu_{j}$. Each column of the state transition matrix corresponds to the state transitions from a specific time point. From this matrix, the joint transition rate is calculated as

[TABLE]

where $\delta_{x,y}$ is the Kronecker delta function. The probability flux is calculated from equation (4). As another example, the stationary distribution is estimated from the state transition matrix

[TABLE]

Because of the Markov assumption, each column of the state transition matrix is independent of each other. Therefore, we can resample each column of the state transition matrix $\bm{K}$ randomly to create the bootstrapped state transition matrix of the same size. The value of interest is calculated from the bootstrapped state transition matrix, and we can estimate the error from the ensemble of realizations. The error of the probability flux is estimated as the standard deviation over the bootstrapped trajectories.

Spatial coarse-graining of brain region through

hierarchical clustering

The neural data has high-dimensionality and defining the coarse-grained state space is not straightforward. Previous work performs the PCA and $k$-means clustering [16] to reduce the dimensionality of data and that of state space. In the spirit of renormalization group of statistical physics [44], we seek the alternative description to reduce the dimensionality of data and state space at once. which is simple but capturing the essence of the system. The standard block spin transformation—coarse-graining the system based on the periodic regular structure of interaction network— is not valid for neural networks: we need to consider the highly irregular and heterogeneous structure of the interaction networks. The core idea is approaching from empirical correlation rather than actual interaction network which is unknown beforehand [45]. If the interaction of the system is local (as expected in cortex), the strongest correlation likely arise from the interaction.

For our purpose of coarse-graining of the human brain neural network system, we perform the hierarchical clustering [46] of the brain regions. Note that, as mentioned in ref. [45], there are many possible methods to perform coarse graining [69], and we emphasize the hierarchical clustering is just one of them. We begin with the correlation matrix between the brain regions averaged over all type of scans, tasks, participants, and time points. The data dimension is $N=100$ and its length is $M=2\times 8\times 590\times 176$, and we calculated the correlation matrix. Then, using the correlation matrix, we perform the hierarchical clustering of the brain regions through unweighted pair group method with arithmetic mean (UPGMA) with Euclidean distance metric. To see the clear leaf structure of dendrogram, we perform the optimal ordering [70] of the linkage matrix. We then threshold the dendrogram by its dissimilarity to define the cluster of brain regions. The neural signal is averaged over the brain regions within each cluster.

Temporal coarse graining through time series

binarization

To discretize and coarse-grain the state space, we apply a binarization technique to the time series data. We employ the symbolic string transformation [71] to convert the continuous time series into the discrete symbol sequence. We first approximate the original $N$-dimensional $M$ data points each-row-standardized data $\bm{X}\in\mathbb{R}^{N\times M}$ by the temporally continuous $N$ spline function written as vector $\bm{f}\left(t\right)\coloneqq\begin{bmatrix}f_{1}\left(t\right)&\cdots&f_{N}\left(t\right)\end{bmatrix}^{\top}\in\mathbb{R}^{N}$. The cubic spline interpolation, where the spline function can be differentiated twice, is employed for time series of each dimension. Using the spline function, we define the three type of binarization distinguished by differentiation. We obtain the event time points (intersection points, stationary points, or inflection points) of $i$th dimensional spline function through the differentiation,

[TABLE]

where the index of event time points is sorted in ascending order, $t_{1}^{\prime}<t_{2}^{\prime}<t_{3}^{\prime}<\cdots<t_{M^{\prime}}^{\prime}$, and $M^{\prime}$ is the number of event time points. Note that when $\alpha=0$, the points are the intersection points, when $\alpha=1$, the points are the stationary points, and when $\alpha=2$, the points are the inflection points. We then discretize the spline function into a time series of binary $s_{i}\left(t\right)$ by three types of transformation: static transformation

[TABLE]

dynamic transformation

[TABLE]

or curve transformation

[TABLE]

Here, $s_{i}\left(t-0\right)\coloneqq\lim_{\Delta t\to 0}s_{i}\left(t-\Delta t\right)$ is the limit approaching from the negative side of $t$, and $s_{i}\left(t^{\prime}\right)=\lim_{\Delta t\to 0}s_{i}\left(t^{\prime}+\Delta t\right)=-s_{i}\left(t^{\prime}-0\right)$. The static and dynamic transformation is suggested in ref. [71] and we add the curve transformation. Finally, we temporally discretize the $N$-dimensional binarized time series $\bm{s}\left(t\right)$ into a state transition matrix:

[TABLE]

where $\mu\left(t\right)=1+\sum_{i=1}^{2^{N}}2^{i-1}\frac{1+s_{i}\left(t\right)}{2}$ is the index of $N$-dimensional Ising state vector $\bm{s}\left(t\right)\coloneqq\begin{bmatrix}s_{1}\left(t\right)&\cdots&s_{N}\left(t\right)\end{bmatrix}^{\top}\in\left\{+1,-1\right\}^{N}$ at time $t$, $\mu\left(0\right)$ is the initial state index, and $\mu\left(\tau\right)$ is the final state index. With the state transition matrix, we perform the probability flux analysis. For analysing neural data, we create a state transition matrix of each task by combining the state transition matrix of each scan and subject.

Examining the assumption of stationary distribution

The assumption behind the probability flux analysis (Probability flux analysis in Methods) is that the probabity distribution is stationary. To validate this assumption, we calculate the probability change of state $\mu$, $\Delta p_{\mu}$, using the estimated probability flux [16]:

[TABLE]

This is zero in steady state thus we investigate the distribution of $\Delta p_{\mu}$ to validate the stationarity of the probability distribution $p_{\mu}$.

Hypercubic probability flux diagram

The $N$-dimensional Ising state vector $\bm{s}$ corresponds to vertex of the hypercube and the state transition involving single spin flip corresponds to the edges of the hypercube. To visualize the probability flux in such the hypercubic state space—or hypercubic probability flux—we need to project hypercubes onto two-dimensional plane. We employ PCA to obtain reproducible, interpretable and automatic projections of hypercubes [49]. With the empirical stationary distribution (equation (7)) of each hypercubic state $\bm{s}\coloneqq\begin{bmatrix}s_{1}&\cdots&s_{N}\end{bmatrix}^{\top}\in\left\{+1,-1\right\}^{N}$, we calculate the covariance matrix

[TABLE]

where $\left<\ast\right>\coloneqq\sum_{\mu=1}^{2^{N}}p_{\mu}\,\ast$ is the average over the stationary distribution. After the diagonalization of the covariance matrix, we obtain $i$th PC loadings $\left\{\bm{v}_{i}\right\}_{i=1}^{N}$ and PC scores $\left\{\bm{r}_{i}\right\}_{i=1}^{N}$. We then introduce the biplot vectors $\left\{\tilde{\bm{e}}_{i}\right\}_{i=1}^{N}$ of PC$j$ and PC$k$ as,

[TABLE]

where $v_{i;j}$ and $v_{i;k}$ are the $j$th and $k$th elements of the $i$th PC loading $\bm{v}_{i}$, respectively. The biplot vectors are the projection of the unit vectors of the original high-dimensional space onto the chosen PCs. Typically, we employ first two PCs. Using the biplot vectors, the probability fluxes are visualized as hypercubic arrows [49], where the width is proportional to the magnitude $\left|\mathcal{J}_{\mu,\nu}\right|$ and the direction corresponds to the sign of the probability flux $\operatorname*{sgn}\left(\mathcal{J}_{\mu,\nu}\right)=-\operatorname*{sgn}\left(\mathcal{J}_{\nu,\mu}\right)$. The direction of the arrow aligns with the direction of the probability flux.

Infering Ising spin system from transition rates

Here, we describe the procedure to infer the Ising spin system from the transition rate of master equation. To begin with, we derive the transition rate of the Ising spin system with symmetric interaction. We then apply the derived transition rate to infer the Ising spin system with asymmetric interaction.

The (pseudo-)Hamiltonian of the Ising spin system is given by the interaction and the external input:

[TABLE]

where $\bm{s}\coloneqq\begin{bmatrix}s_{1}&\cdots&s_{N}\end{bmatrix}^{\top}\in\left\{+1,-1\right\}^{N}$ is the Ising state of the system consisting of $N$ components, $\bm{J}\in\mathbb{R}^{N\times N}$ is the (not necessarily symmetric) interaction matrix where element $J_{i,j}$ is the interaction from component $j$ to $i$, and $\bm{h}\coloneqq\begin{bmatrix}h_{1}&\cdots&h_{N}\end{bmatrix}^{\top}\in\mathbb{R}^{N}$ is the external input where element $h_{i}$ is the external input to component $i$. The self-interaction is set to zero, i.e., $J_{i,i}\coloneqq 0$ for all $i$. The positive interaction, $J_{i,j}>0$, means it is excitatory (or ferromagnetic) interaction, and the negative interaction, $J_{i,j}<0$, means it is inhibitory (or antiferromagnetic) interaction. The same goes for the external input $\bm{h}$: the positive external input, $h_{i}>0$, means it is excitatory bias, and the negative external input, $h_{i}<0$, means it is inhibitory bias. If the interaction matrix is symmetric, i.e., $\bm{J}=\bm{J}^{\top}$, equation (16) is the Hamiltonian, but if the interaction matrix is asymmetric $\bm{J}\neq\bm{J}^{\top}$, equation (16) is the pseudo-Hamiltonian because the energy is ill-defined. Below we assume the symmetric interaction matrix, to derive the transition rate of the Ising spin system. Our goal is to infer the interaction matrix $\bm{J}$ and external input $\bm{h}$ from the estimated joint transition rates $\left\{w_{\mu,\nu}\right\}$ rather from the empirical probability distribution $\left\{p_{\mu}\right\}$.

From the probability flux analysis, we obtain the transition rate from state $\nu$ to $\mu$ by dividing the joint transition rate by the stationary distribution $p_{\nu}$:

[TABLE]

The stationary distribution is estimated from equation (7) Thus, we can estimate the transition rate from empirical data.

We then build a stochastic model to fit the observation. Assuming the continuous time Markov process, we introduce time-evolution of the probability distribution by employing the master equation [53],

[TABLE]

where $p_{\mu}\left(t\right)$ is the probability of finding the system in state $\mu$ at time $t$, and $w_{\mu,\nu}$ is the transition rate from state $\nu$ to $\mu$. The transition rate is determined by the detailed balance condition,

[TABLE]

where $p_{\nu}\coloneqq\lim_{t\to\infty}p_{\nu}\left(t\right)$ is the stationary distribution of the system. Unlike the balance condition, $\sum_{\nu=1}^{2^{N}}w_{\mu,\nu}p_{\nu}=\sum_{\nu=1}^{2^{N}}w_{\nu,\mu}p_{\mu}$, the detailed balance condition constraint the system being microscopically reversible. From the detailed balance condition, the ratio of forward and backward transition rates is given by the ratio of the stationary distributions:

[TABLE]

Assuming that the stationary distribution is the equilibrium distribution, the transition rate is determined. With the canonical ensemble, the equilibrium distribution $\pi_{\mu}$ is given by

[TABLE]

where $\beta\coloneqq\frac{1}{k_{\mathrm{B}}T}$ is the inverse temperature with $k_{\mathrm{B}}$ being the Boltzmann constant and $T$ the absolute temperature, $\bm{s}_{\mu}$ is the Ising state vector of state $\mu$, and $\mathcal{Z}\coloneqq\sum_{\mu=1}^{2^{N}}\exp\left[-\beta\mathcal{H}\left(\bm{s}_{\mu}\right)\right]$ is the normalization constant or partition function of statistical mechanics. Note that each Ising state vector $\bm{s}$ is uniquely identified as corresponding integer, $\mu=1+\sum_{i=1}^{N}2^{i-1}\frac{1+s_{i}}{2}$. Substituting the equilibrium distribution $\pi_{\mu}$ into the detailed balance condition, we obtain

[TABLE]

where $\Delta E_{\mu,\nu}\coloneqq\mathcal{H}\left(\bm{s}_{\mu}\right)-\mathcal{H}\left(\bm{s}_{\nu}\right)=-\Delta E_{\nu,\mu}$ is the energy difference from state $\nu$ to $\mu$. We note two transition rates: the Arrhenius transition rate

[TABLE]

and Glauber transition rate [53]

[TABLE]

where $A\in\mathbb{R}_{>0}$ is the rate constant. Below we use the Arrhenius transition rate for example but extension to the Glauber transition rate is straightforward. For Ising spin system, assuming single spin flip dynamics [53] and the symmetric interaction matrix, the energy difference from state $\mu$ to $\mu^{\left(k\right)}$ by flipping spin $k$ is given as

[TABLE]

where $s_{k;\mu}$ is the $k$th element of Ising state vector $\bm{s}_{\mu}$. The index $\mu^{\left(k\right)}$ is the index of state obtained by flipping the $k$th spin of state $\mu$: $\begin{bmatrix}s_{1;\mu}&\cdots&-s_{k;\mu}&\cdots&s_{N;\mu}\end{bmatrix}^{\top}=\bm{F}_{\left(k\right)}\bm{s}_{\mu}$. Here, $\bm{F}_{\left(k\right)}\coloneqq\bm{I}-2\bm{e}_{k}\bm{e}_{k}^{\top}$ is the spin flip matrix that flips the $k$th spin with $\bm{e}_{k}$ being the $k$th standard unit vector of $\mathbb{R}^{N}$. If the interaction matrix is symmetric, equation (26) is derived [49] from the Hamiltonian (16) but when the interaction matrix is asymmetric, the derivation is not possible. Nevertheless, one can interpret the equation (26) is a different way. The term $\sum_{i=1}^{N}J_{k,i}s_{i;\mu}+h_{k}$ is the effective field or input to spin $k$ of state $\mu$, hence the energy difference $\Delta E_{\mu^{\left(k\right)},\mu}$ drive the spin $k$ to align with the effective field. With asymmetric interaction matrix, this interpretation is valid: the spin $k$ flipping is determined by the incoming effective field. We assume the symmetric interaction matrix to derive equation (26) but we apply the result to the asymmetric interaction matrix [16]. From eqs. (24) and (26), the transition rate from state $\mu$ to $\mu^{\left(k\right)}$ by flipping spin $k$ is given by

[TABLE]

We set transition rate to zero between states where the difference of them is not a single spin flip. i.e., $w_{\mu,\nu}=0$ if $\mu\neq\nu^{\left(k\right)}$.

We then fit the model parameters of interaction $\bm{J}$, external input $\bm{h}$, and transition rate constant $A$ to the empirical observation. We define the loss function, which is the difference of the transition rates of models and that from data $\sum\left[\ln\left(w^{\text{(model)}}\right)-\ln\left(w^{\text{(data)}}\right)\right]^{2}$, as a function of these parameters,

[TABLE]

where the inverse temperature $\beta$ is the nuisance parameter and the sum $\sum_{\left(\mu^{\left(k\right)},\mu\right)}$ is over all pairs of states with single spin flip difference; there are $N2^{N}$ such pairs. The loss function is minimized when the logarithm of transition rates of model (equation (27)) match the that of empirical observations (equation (17)). The parameters minimizing the loss of function of equation (28) gives the fitted parameters: $\beta\tilde{\bm{J}}$, $\beta\tilde{\bm{h}}$, and $\tilde{A}$.

[TABLE]

The estimated interaction matrix is in general asymmetric, $\tilde{\bm{J}}\neq\tilde{\bm{J}}^{\top}$. Note that the number of inferred parameters is $\left(N^{2}-N\right)+N+1=N^{2}+1$ which is smaller than the number of constraints $N2^{N}$, and this difference constrains the system enough.

The steady state hypercubic probability flux diagrams are reconstructed from the inferred model as below. Using the estimated parameters of equation (29), we first calculate the transition rates of equation (27). Then we rewrite the master equation (equation (18)) as matrix–vector multiplication form:

[TABLE]

Here, $\bm{p}\left(t\right)\coloneqq\begin{bmatrix}p_{1}\left(t\right)&\cdots&p_{2^{N}}\left(t\right)\end{bmatrix}^{\top}\in\left[0,1\right]^{2^{N}}$ is the probability vector of the system at time $t$, which is the statistical state of the system. The element of the transition rate matrix $\bm{W}\in\mathbb{R}^{2^{N}\times 2^{N}}$ is defined by

[TABLE]

Because the transition rate matrix is the stochastic matrix, from Perron–Frobenius theorem, the largest eigenvalue of the transition rate matrix is zero and the corresponding eigenvector is the unnormalized stationary distribution vector $\bm{p}$. With the stationary distribution $p_{\mu}$, we calculate the probability flux $\mathcal{J}_{\mu,\nu}\coloneqq w_{\mu,\nu}p_{\nu}-w_{\nu,\mu}p_{\mu}$ to validate the inferred model with the empirical data.

Entropy production rate

The entropy production rate of the system governed by the master equation (equation (18)) is given by [50]

[TABLE]

When the detailed balance condition is satisfied, $\dot{\mathcal{S}}_{\text{tot}}\left(t\right)=0$. In general, the entropy production rate is non-negative, $\dot{\mathcal{S}}_{\text{tot}}\left(t\right)\geq 0$, which corresponds to the second law of thermodynamics. For more details, see the Supplementary Information.

Data availability

Previously collected [24, 23] and preprocessed [16] data are used in this study. All data are available online from Zenodo [72]. https://doi.org/10.5281/zenodo.0000000.

Code availability

Calculations and visualizations of this work were performed using open-source Python [73] libraries: BrainSpace [74], JAX [75], Matplotlib [76], NetworkX [77], neuromaps [78], nilearn [79], Numba [80], NumPy [81], pandas [82], scikit-learn [83], SciPy [84], and seaborn [85]. The colour map of some figures are generated by ColorCET [86]. All code are available online from Zenodo [72]. https://doi.org/10.5281/zenodo.0000000.

Acknowledgements This work was supported by KAKENHI Grant Number 22H00406 and 24H00061 of the Japan Society for the Promotion of Science. Y.H. is supported by JST SPRING, Grant Number JPMJSP2125 “THERS Make New Standards Program for the Next Generation Researchers”.

Author contributions Y.H. performed conceptualization, data curation, formal analysis, investigation, software development, validation, and visualization; Y.H. partially contributed to the funding acquisition; Y.H. and S.F. developed the methodology; S.F. provided computational resources; Y.H. wrote the original draft; and Y.H. and S.F. reviewed and edited the manuscript.

Competing interests All authors declare no competing interests.

**Additional information

Supplementary information** The online version contains supplementary material.

Correspondence and requests for materials should be addressed to Y.H.