Parcels and particles: Markov blankets in the brain

A persistent theme in systems neuroscience, especially neuroimaging, is the search for principles that underlie the functional anatomy of distributed neuronal processes. These principles are usually articulated in terms of functional segregation (or differentiation) and integration, which inherit from centuries of neuroanatomical, neurophysiological, and neuropsychological study (Zeki & Shipp, 1988). In recent thinking about functional integration, people have turned to formal accounts of (predictive) processing in the brain (e.g., Bastos et al., 2012; Keller & Mrsic-Flogel, 2018; Parr & Friston, 2018; Rao & Ballard, 1999; Spratling, 2008) to understand the nature of (neuronal) message passing on graphs, where edges correspond to connectivity and nodes correspond to neuronal populations. Crucially, this characterization rests upon the asymmetric and directed connectivity that defines cortical and subcortical hierarchies (e.g., Bastos et al., 2012; Crick & Koch, 1998; Felleman & Van Essen, 1991; K. J. Friston, Parr, & de Vries, 2017; Keller & Mrsic-Flogel, 2018; N. Markov et al., 2013; N. T. Markov et al., 2014; Mesulam, 1998; Stachenfeld, Botvinick, & Gershman, 2017; Zeki & Shipp, 1988). Usually, these asymmetries are expressed in terms of things like laminar specificity that distinguish between forward and backward connections (Buffalo, Fries, Landman, Buschman, & Desimone, 2011; Grossberg, 2007; Haeusler & Maass, 2007; Hilgetag, O’Neill, & Young, 2000; Thomson & Bannister, 2003; Trojanowski & Jacobson, 1977). More recently, asymmetries in spectral content have become an emerging theme (Arnal & Giraud, 2012; Bastos et al., 2015; Buffalo et al., 2011; Giraud & Poeppel, 2012; Hovsepyan, Olasagasti, & Giraud, 2018; Self, van Kerkoerle, Goebel, & Roelfsema, 2019; Singer, Sejnowski, & Rakic, 2019; van Kerkoerle et al., 2014).

In contrast, analyses of functional connectivity have focused on distributed patterns of coherent fluctuations in neuronal activity and phenomenological descriptions of the implicit dynamics (Bassett & Sporns, 2017; Biswal, Van Kylen, & Hyde, 1997; Bullmore & Sporns, 2009; Gilson, Moreno-Bote, Ponce-Alvarez, Ritter, & Deco, 2016; Gu et al., 2018; Lynall et al., 2010; van den Heuvel & Sporns, 2013). This phenomenology ranges from intrinsic brain networks—which are conserved over subjects in resting-state functional magnetic resonance imaging—to the dependence of neuronal dynamics on cortical excitability (Freyer, Roberts, Ritter, & Breakspear, 2012; Roy et al., 2014). The principles that are brought to bear on this kind of characterization could be seen as ascribing neuronal dynamics to various universality classes, such as self-organized criticality (Bak, Tang, & Wiesenfeld, 1988; Breakspear, Heitmann, & Daffertshofer, 2010; Cocchi, Gollo, Zalesky, & Breakspear, 2017; Deco & Jirsa, 2012; Haimovici, Tagliazucchi, Balenzuela, & Chialvo, 2013; Kitzbichler, Smith, Christensen, & Bullmore, 2009; Lopez, Litvak, Espinosa, Friston, & Barnes, 2014; Shin & Kim, 2006). (Note: Although we have subsumed criticality and dynamic instability under phenomenological approaches, criticality can refer to the dynamics of neurons and neural assemblies—as opposed to the statistical properties of data leveraged by functional connectivity. A simple example of criticality is a branching process, an inherently directed process. Neurobiological models of these kinds of processes have been derived from causal neural field models with directed (corticothalamic) interactions [Freyer et al., 2011] and with coupled oscillators [Deco & Jirsa, 2012]. In this setting, criticality acquires a more mechanistic aspect, as we will see below.) This dual-pronged approach to functional integration invites an obvious question: Is there a way of linking the two?

Practically, the study of context-sensitive, directed coupling between the nodes of neuronal networks calls for an estimate of effective connectivity, under some model of how measured brain signals are generated. One then has to resolve the ill-posed problem of recovering the underlying (connectivity) parameters of the model, usually using Bayesian inference. The best example here is dynamic causal modeling (K. J. Friston, Harrison, & Penny, 2003; K. J. Friston, Kahan, Biswal, & Razi, 2014; Razi & Friston, 2016). The complementary approach—based upon functional connectivity—borrows ideas from network science and graph theory. This entails specifying an adjacency matrix, usually formed by thresholding a functional connectivity matrix summarizing dependencies among nodes, where the nodes are generally defined in terms of some parcellation scheme (Bassett & Sporns, 2017; Bullmore & Sporns, 2009).

In what follows, we will consider a particular parcellation scheme based upon effective connectivity and ask whether it leads to the same phenomenology seen in network neuroscience. In doing so, we can, in principle, explain and quantify the emergence of large-scale intrinsic brain networks and their characteristic dynamics. A crucial aspect of the particular parcellation or partition—employed in this work—means that it can be applied recursively in the spirit of the renormalization group (Schwabl, 2002). This means that there is a formal way of quantifying the dynamics at various spatiotemporal scales. Our hypothesis was that the spatiotemporal dynamics of coarser scales would evince both the functional anatomy of intrinsic brain networks and the emergence of (self-organized) criticality—as assessed in terms of dynamical instability.

Although this work is framed as addressing issues in network neuroscience (Bassett & Sporns, 2017), it was originally conceived as a parcellation scheme for multiscale analyses of neuroimaging time series. In other words, it was intended as a first principle approach to dimension reduction and decomposition, as a prelude for subsequent graph theoretic or dynamic causal modeling (K. J. Friston, Kahan, et al., 2014; Razi, Kahan, Rees, & Friston, 2015; Razi et al., 2017; Zhou et al., 2018). However, the theoretical foundations—and uniqueness of the partition—proved too involved to support a simple and practical procedure. Instead, what follows is offered as a case study of emergence in coupled dynamical systems, using the brain as a paradigm example.

This paper comprises five sections. In the first, we review the notion of Markov blankets and how recursive applications of a partition or parcellation of states into Markov blankets allows one to express dynamics at increasing scales. We will use the notion of the renormalization group (RG) to motivate this recursive parcellation because there are some formal constructs (in terms of RG scaling) that furnish an insight into how dynamics change as we move from one scale to the next. The second section describes a simple (dynamic causal modeling) analysis of directed effective connectivity at the finest spatial scale, as summarized with a Jacobian. This plays the role of a directed adjacency matrix, which is all that is needed for successive renormalization to higher scales. The renormalization group is illustrated with an exemplar dataset, to show what the ensuing parcellation scheme looks like. This section concludes with a brief consideration of sparse coupling at the finest scale, in terms of excitatory and inhibitory connections. The subsequent sections consider dynamics at different scales of parcellation, in terms of intrinsic (within parcel) and extrinsic (between parcel) connectivity. Our focus here is on the progressive slowing of intrinsic dynamics as we move from one scale to the next—a slowing that organizes the dynamics at coarser (higher) scales towards critical regimes of instability and slowly fluctuating dynamical modes. The third section illustrates the emergence of autonomous dynamics, in terms of characteristic frequencies associated with intrinsic connectivity, and in terms of positive Lyapunov exponents that speak to transcritical bifurcations at, and only at, coarser scales. The fourth section focuses on extrinsic connectivity and the coupling between (complex) modes or patterns of activity and how this relates to functional connectivity and intrinsic brain networks (Fox et al., 2005). The final section reviews the dynamical phenomenology at hand from the point of view of statistical physics, with a special focus on dissipative dynamics and detailed balance at nonequilibrium steady state. We conclude with a brief discussion and qualification of this particular (sic) approach to functional integration.

The last section concluded with reference to a particular partition. The use of the word “particular” has a double entendre here. It is predicated on a more fundamental (or perhaps foundational) analysis of coupled dynamical systems that consider the emergence of “particles.” Full details of this treatment can be found in K. Friston (2019). From the current perspective, we just need to know how to define Markov blankets (Clark, 2017; Kirchhoff, Parr, Palacios, Friston, & Kiverstein, 2018; Pearl, 1988; Pellet & Elisseeff, 2008) and how Markov blankets engender particles and particular partitions (K. Friston, 2019). For readers interested in Markov blankets for dynamical systems, fairly comprehensive discussions can be found in K. Friston, Da Costa, and Parr (2020) and Parr, Da Costa, and Friston (2020).

In brief, a Markov blank et allows one to distinguish a collection of vector states (hereafter, simply states) that belong to a particle from states that do not. This provides an operational definition of a particle that, in the present setting, can be regarded as a region of interest or parcel of brain states. This means that a particular partition becomes a parcellation scheme, in terms of functional anatomy. The particular partition refers to a partition of a (potentially large) set of states into a smaller number of particles, where each particle is distinguished from other particles, in virtue of possessing a Markov blanket. A Markov blanket is simply a set of states that separate or insulate—in a statistical sense—states that are internal to the blanket and states that are on the outside; namely, external states. Technically, this means that internal states are conditionally independent of external states, when conditioned upon their blanket states (Pearl, 2009).

In a particular partition, all external states are assigned to particles, to create an ensemble of particles that are constituted by their blanket states and the internal states within or beneath the blanket. The crucial aspect of this partition is that we only need the blanket states to understand coupling between particles. This follows from the conditional independence between internal and external states, where the external states “that matter” are the blanket states of other particles. In short, the particular partition is a principled way of dividing states into particles or parcels that is defined in terms of statistical dependencies or coupling among states. In more complete treatments, one can divide the blanket states into active states and sensory states, according to the following rules: Sensory states are not influenced by internal states, while active states are not influenced by external states. Indeed, it is the absence of these influences that enables us to identify the Markov blanket of any given set of internal states. Please see the Supporting Information for a formal definition of Markov blankets in this dynamical context.

As noted above, we are dealing with vector states (not scalar variables). So, what is a vector state? A vector state is the multidimensional state of a particle, for example, the principal eigenstates of its Markov blanket. However, we have just said that a particle arises from a partition of states—and now we are saying that a state is an eigenstate (i.e., a linear mixture) of the blanket states of a particle. So, is a particle a collection of states or is a state the attribute of a particle (i.e., its blanket states)? The answer is both, because we have particles at multiple levels.

In our case, the random variables are states; the coarse-graining operation corresponds to the grouping into a particular partition (G) and a dimension reduction (R) inherent in retaining the principal eigenstates of particular blanket states. The dimension-reduction operator (R) has two parts. First, we can eliminate the internal states because they do not contribute to coupling between particles. Second, we can eliminate the eigenstates that dissipate very quickly; namely, those with large negative eigenvalues. These are the fast or stable modes of a dynamical system (Carr, 1981; Haken, 1983). This leaves us with the slow, unstable eigenstates picked out by the dimension reduction, which we can now see as an adiabatic approximation. Please note that in quantum mechanics, the adiabatic approximation refers to those solutions to the Schrödinger equation that make use of a timescale separation between fast and slow degrees of freedom.

In this (RG) setting, a relevant variable is said to describe the macroscopic behavior of the system. From our perspective, the relevant variables in question correspond to the slow eigenstates. In short, we can reduce many states to a small number of eigenstates that summarize the dynamics “that matter.” These eigenstates are the relevant variables that underwrite critical slowing. Figures 1 and 2 provide a graphical illustration of this recursive partitioning and reduction based upon an adiabatic approximation (i.e., eliminating fast eigenstates and approximating dynamics with the remaining slow eigenstates). This adiabatic reduction is commonplace in physics, where it plays a central role in synergetics through the enslaving principle (Haken, 1983) and, in a related form, in the center manifold theorem (Carr, 1981).

We now have at hand a principled procedure to repeatedly coarse-grain a system of loosely coupled particles (e.g., nonlinear neuronal oscillators) at successively coarser spatiotemporal scales. One can see that, by construction, as we ascend scales, things will get coarser and slower. It is this progressive slowing towards criticality that is the primary focus of the examples pursued below. However, before we can apply the particular partition to some empirical data, we first need to quantify the coupling among particles at a suitably fine or small scale. Having characterized this coupling in terms of some dynamical system or state space model, we can then use the Jacobian to identify a particular partition, compute the Jacobian of the blanket states, and then take the ensuing eigenstates to the next level, as described above. This furnishes a description of dynamics in terms of the intrinsic (within particle) coupling (i.e., eigenvalues) of any particle $λnn(i)$ and their extrinsic (between particle) coupling $λnm(i)$⁠. We will unpack the meaning of these terms later using numerical examples and analysis. At present, we will focus on estimating the coupling among a large number of particles at the smallest scale.

This linearized system can now be solved using standard (variational Laplace) schemes for parametric empirical Bayesian (PEB) models, to provide (approximate) Gaussian posteriors over the unknown elements of the Jacobian—and the unknown covariance parameters encoding the amplitude of various random effects (K. Friston, Mattout, Trujillo-Barreto, Ashburner, & Penny, 2007). This Bayesian model inversion requires priors on the parameters and hyperparameters (i.e., covariance component parameters), specified as Gaussian shrinkage priors. For nonnegative hyperparameters, Gaussian shrinkage priors are generally applied to log-transformed hyperparameters (i.e., a lognormal prior over nonnegative scale parameters).

Equipped with posterior densities over the coupling parameters—or elements of the Jacobian—we can now use Bayesian model reduction to eliminate redundant parameters (K. J. Friston et al., 2016); namely, parameters whose shrinkage to zero increases model evidence by removing redundancy or complexity. As described elsewhere (K. Friston & Penny, 2011), this can be done very efficiently, because we know the form of the posteriors, before and after reducing the model. Furthermore, we can apply other prior constraints to eliminate redundant coupling parameters.

In the examples below, we performed Bayesian model reduction to enforce reciprocal coupling among states, given that extrinsic connections in the brain are almost universally recurrent (Felleman & Van Essen, 1991; N. Markov et al., 2013; N. T. Markov et al., 2014). This was implemented by combining the changes in variational free energy—or log model evidence—when removing connections between two states in both directions. If model evidence increased by three natural units (i.e., a log odds ratio of exp(3):1 or 20:1), both connections were removed but not otherwise. In addition, we precluded long-range coupling (beyond 32 mm) and used Bayesian model reduction to identify the most likely upper bound on the spatial reach of coupling between nonhomologous particles (i.e., particles that did not occupy homologous positions in each hemisphere). These empirical connectivity priors were based upon a body of empirical work, suggesting that the density of axonal projections—from one area to another—declines exponentially as a function of anatomical separation (Ercsey-Ravasz et al., 2013; Finlay, 2016; Horvát et al., 2016; Wang & Kennedy, 2016). We will later examine the evidence for this kind of distance rule, based upon coupling among particles at the finest scale.

In summary, Equation 7 is used to evaluate the effective connectivity of a dynamic causal model based upon the Jacobian of a stochastic differential equation under some simplifying assumptions. In brief, we start with a linear state space model, in which the response variable y (the multivariate BOLD time series) is a linear convolution (K) of some hidden states x subject to observation and system noise ω_y and ω_x, respectively (cf. Equation 5). We can linearize this state space model (cf. Equation 6) such that it can be written in matrix form as a general linear model (cf. Equation 7). This linearized system is then solved using standard variational Laplace for parametric empirical Bayes (PEB) that provides the (Gaussian) posterior estimates for the system parameters (elements of the Jacobian) and hyperparameters (the unknown covariance of the observation and system noise). Since this scheme uses PEB for model inversion, it is automatically protected against becoming trapped in local minima. PEB uses the formal apparatus of variational Laplace, which optimizes a free energy lower bound, which optimizes the trade-off between model complexity and accuracy. Given the posterior distributions over the system parameters (and hyperparameters) and the model evidence, one can then use Bayesian model reduction procedures to prune away any redundant couplings. This constitutes a computationally efficient inversion scheme that can invert very large systems within minutes on a standard laptop. At the finest scale, the Jacobian had 1,024 by 1,024 elements taking about 45 min to infer the model parameters.

Computationally, the benefit of linearizing the system in this way means that one can evaluate the posterior coupling parameters or elements of the Jacobian region by region (cf. Frassle et al., 2017). This means that, provided one is prepared to wait long enough, one can invert large systems with thousands of regions or parcels. On a personal computer, it takes about an hour to evaluate the Jacobian for coupling among 1,024 states. To illustrate the renormalization group procedure practically, we applied it to the fMRI time series from a single subject. These time series are the same data used to illustrate previous developments in dynamic causal modeling. Time series data were acquired from a normal subject at 2 Tesla using a Magnetom VISION (Siemens, Erlangen) whole-body MRI system. Contiguous multislice images were acquired with a gradient echo-planar sequence (TE = 40 ms; TR = 3.22 s; matrix size = 64 × 64 × 32, voxel size 3 × 3 × 3 mm). Four consecutive 100-scan sessions were acquired, comprising a sequence of 10-scan blocks under five conditions. The first was a dummy condition to allow for magnetic saturation effects. In the second, Fixation, the subject viewed a fixation point at the center of the screen. In an Attention condition, the subject viewed 250 dots moving radially from the center at 4.7 degrees per second and was asked to detect changes in radial velocity. In No attention, the subject was asked to look at the moving dots. In the last condition, subject viewed stationary dots. The order of the conditions alternated between Fixation and photic stimulation. The subject fixated on the center of the screen in all conditions. No overt response was required in any condition, and there were no actual speed changes. In contradistinction to normal procedures in functional connectivity fMRI analyses, the time series were not smoothed (other than adjusting for ultraslow scanner drifts). This is because the random fluctuations at all timescales play a material role in the decomposition at hand. Informed consent from the subject was obtained and the study was approved by the Human Ethics Review Committee of University College London.

In the exemplar analyses below, we started at a scale where each particle can be plausibly summarized with a single state. This single state was the principal eigenstate following a principal components analysis of voxels that lay within about 4 mm of each other. This can be thought of as reducing the dynamics to a single mode of the Markov blanket of this small collection of voxels. Practically, this simply involved taking all voxels within a fixed radius of the voxel showing the largest variance, performing a singular value decomposition, and recording the first eigenvariate. These voxels were then removed, and the procedure repeated until the entire multivariate time series was reduced to 1,024 eigenstates, where each eigenstate corresponds to a simple particle. See Figure 3. Clearly, we could have summarized the dynamics of each collection of voxels with two or more eigenstates; however, for simplicity we will assume that the eigenstate with the greatest variance is a sufficient summary of the slow, non-dissipative, dynamics of this smallest scale. Interestingly, this variance is proportional to the characteristic time constant of systemic dynamics; namely, the negative inverse of the eigenvalues of the underlying Jacobian (see the final section). In other words, as the (negative) principal eigenvalue of effective connectivity approaches zero from below, the principal eigenvalue of functional connectivity (i.e., variance) increases; see Equation 9 in Lopez et al. (2014).

Following Bayesian model reduction, we now have a sparse Jacobian or directed, weighted adjacency matrix describing the dynamical coupling between univariate states of 1,024 particles (see Figure 4). This Jacobian can now be subject to a particular partition to identify the most connected internal states and their Markov blanket, following the procedures summarized in Figure 2. This grouping process (i.e., the G operator) was repeated until all 1,024 states are accounted for: in this example, grouping 1,024 states into 57 particles. For simplicity, and for consistency with the first level, each particle was assigned a single internal state. The ensuing partition was then subject to an adiabatic reduction (i.e., the R operator) by taking the eigenvectors of the Jacobian describing the intrinsic dynamics of each particle’s blanket states.

The eight principal eigenstates were retained if their eigenvalues were less than 1. This is the adiabatic approximation that dispenses with modes that dissipate quickly, here, within a second. This reduces the intrinsic coupling to a diagonal matrix $λnn(i)$⁠, corresponding to the eigenvalues of the intrinsic Jacobian ∂_xn$fn(i)$⁠. The extrinsic coupling $λnm(i)$ contains complex elements that couple the eigenstates of one particle to the eigenstates of another. We will return to how these Jacobians manifest in terms of connectivity later.

Note that it would have been possible to reevaluate the Jacobian using another dynamic causal model of the eigenstates at any particular level and then use Bayesian model reduction to eliminate redundant coupling parameters. This is an interesting alternative to using the estimates of the Jacobian based upon the first-order approximation at the smallest scale. We will explore the impact of reevaluating the Jacobian in subsequent work. For the purposes of the current illustration, we will retain the linear solutions at higher scales—based upon the finest scale—to illustrate that one can still reproduce the emergent phenomena of interest described below. These dynamical phenomena are therefore directly attributable to local linear coupling with a particular sparsity structure that is sufficient to produce interesting self-organized dynamics at higher scales. Before taking a closer look at dynamics over scales, this section concludes with a brief characterization of coupling at the smallest scale.

The Jacobian from the above analysis summarizes the effective connectivity at the smallest scale, where each node particle has a reasonably precise anatomical designation. This means that we can interpret the elements of the Jacobian in terms of directed (effective) connectivity. We had anticipated that this would mirror the exponential distance rule established through anatomical tracing studies (Finlay, 2016; Horvát et al., 2016; Wang & Kennedy, 2016). However, it did not. Instead, this (linear) characterization of effective connectivity was better explained with a power law that, interestingly, was quantitatively distinct for inhibitory (i.e., negative) and excitatory (i.e., positive) connections (i.e., elements of the Jacobian).

Figure 9 summarizes the statistical characteristics of coupling among particles at the first level. The upper left panel shows each connection in terms of the real part of the corresponding Jacobian in hertz, against the distance spanned by the connection (i.e., Euclidean distance between the centers of the two particles). Two things are evident from this scatterplot: first, positive (excitatory, red dots) connections dominate in the short range (around 8 mm), while negative (inhibitory, blue dots) dominate around 16 mm. Although there is variability, the dependency of the coupling strength on distance shows some lawful behavior that is disclosed by plotting the log-coupling (real part) against log-distance (upper right panel). The circles are the averages in bins (discrete ranges) of the dots in the upper left panel. A linear regression suggests a scaling exponent of −1.14 for excitatory coupling and a smaller scaling exponent of −0.52 for inhibitory coupling. This dissociation is consistent with a Mexican hat–like coupling kernel, with short-range excitation and an inhibitory penumbra. This kind of architecture predominates in neural field models of cortical and subcortical coupling (e.g., Coombes, 2005; Itskov, Curto, Pastalkova, & Buzsaki, 2011; Moran, Pinotsis, & Friston, 2013).

The lower panel plots the strength of reciprocal connections against each other, to illustrate the relative proportions of recurrent excitatory and inhibitory coupling, here 65% and 31%, respectively. There are only about 4% of connections that show an antisymmetry, that is, excitatory in one direction and inhibitory in the other. The rarefied region in the center of this scatterplot reflects the fact that connections with small coupling strengths have been eliminated during Bayesian model reduction (see Figure 4). In the next section, we will see how this sparse local coupling engenders progressively more structured and itinerant dynamics at increasing spatial and temporal scales.

This section focuses on the intrinsic dynamics of each particle at different scales by associating the Jacobian of each particle with Lyapunov exponents. For people not familiar with dynamical systems theory, the Lyapunov exponents score the average exponential rate of divergence or convergence of trajectories in state space (Lyapunov & Fuller, 1992; Pavlos et al., 2012; Yuan, Ma, Yuan, & Ao, 2011). Because we are dealing with a linearized system, the Lyapunov exponents are the same as the eigenvalues of the Jacobian describing intrinsic coupling. By construction, this is a leading diagonal matrix containing intrinsic eigenvalues whose real values are close to zero. In terms of a linear stability analysis, the real part of these eigenvalues (i.e., self-induced decay) corresponds to the rate of decay. This means that as the eigenvalue approaches zero from below, the pattern of activity encoded by this eigenstate decays more and more slowly. This is the essence of critical slowing and means that, from the point of view of dynamical stability, this eigenstate has become unstable (Haken, 1983; Jirsa, Friedrich, Haken, & Kelso, 1994; Mantegna & Stanley, 1995; Pavlos et al., 2012). The complement of critical instability is a stable fast eigenstate that decays very quickly, that is, an eigenstate whose eigenvalue has a large negative real part.

The imaginary part of the eigenvalue describes the characteristic frequency at which this decay is manifest. If the imaginary part is zero, the system decays monotonically. However, complex values mean that the intrinsic dynamics acquire a sinusoidal aspect. Because each particle has a number of eigenstates, an ensemble of particles can be construed as loosely coupled phase oscillators (Breakspear & Stam, 2005; De Monte, d’Ovidio, & Mosekilde, 2003; Kaluza & Meyer-Ortmanns, 2010; Kayser & Ermentrout, 2010), featuring multiple frequencies. The associated dynamics of a single particle can be visualized by plotting its eigenvalues in the complex number plane. The closer the eigenvalue to the vertical axis, the slower the dynamics, such that as the real eigenvalue approaches zero (i.e., from the left), the particle approaches a transcritical bifurcation (at zero) and displays a simple form of critical slowing.

This characterization of intrinsic dynamics—at different scales—is illustrated in the right panels of Figure 10. Note that the complex values are symmetrically paired (dots of the same color). The key thing to observe here is that when we look at the eigenvalues of particles at higher scales, there are some eigenstates that approach criticality and start to show intrinsic oscillatory behavior. This is one of the key observations from the current renormalization scheme; namely, there is a necessary slowing as one progresses from one scale to the scale above. Furthermore, at higher scales intrinsic dynamics start to appear with progressively slower frequencies.

Another way of characterizing temporal dynamics is to use linear systems theory to map the eigenvalues to the spectral density of the time series that would be measured empirically. This rests upon standard transforms and the convolution theorem that enables us to express the systems’ first-order kernels as a function of the Jacobian (Lopez et al., 2014). In frequency space, these kernels correspond to transfer functions and describe the spectral power that is transferred from the random fluctuations to the macroscopic dynamics of each eigenstate. The left panels of Figure 10 show the transfer functions of the eigenstates of each particle at different scales. At the finest scale, power is spread over a large range of frequencies. At progressively higher scales, the power becomes more concentrated in the lower frequencies with a transfer function that has a characteristic Lorentzian form. Crucially, the frequencies at the highest scale correspond to the characteristic ultraslow frequencies studied in resting-state fMRI; namely, < 0.01 Hz. This is an interesting observation, which suggests that one can explain ultraslow fluctuations in resting-state fMRI purely in terms of local directed coupling among small particles of brain tissue. Note that this explanation does not involve any hemodynamics because the Jacobian that gives rise to these slow oscillations pertains to the neuronal states (prior to hemodynamic convolution). In other words, this is not an artefact of removing fast frequencies from the measured fMRI signals.

One might ask whether there is any systematic variation of these ultraslow frequencies across the brain. Figure 11 reports the implicit intrinsic timescales at intermediate scales (second scale, upper rows; third scale, lower rows). The left column shows the eigenmodes in terms of their principal frequency, that is, the largest complex eigenvalue (divided by 2π). The right column shows the corresponding eigenmodes in terms of their principal time constants, that is, the reciprocal of the largest negative real part. These two characterizations—principal frequency and time constant—speak to different aspects of intrinsic timescales, both of which contribute to the shape of an eigenstate’s transfer function. The first quantifies the frequency of solenoidal flow, while the second reflects the rate of decay associated with the dissipative flow. We will unpack solenoidal and dissipative flows in the last section, in terms of density dynamics. In brief, the flows of a random dynamical system, at nonequilibrium steady state, can be decomposed into two orthogonal components. The first (dissipative) component is a gradient flow, in the direction of increasing density. This flow counters the dispersive effect of random fluctuations. The solenoidal (non-dissipative) component circulates on the isoprobability contours, thereby having no effect on the nonequilibrium steady-state density. Heuristically, if one imagines water flowing down a plughole, the dissipative part is the flow out of the basin, while the solenoidal part accounts for the circular motion of water. Generally speaking, when the random fluctuations are small, the gradient flows disappear (at steady state), leaving the solenoidal flow, that is, the classical mechanics that would apply to heavenly bodies. However, when random fluctuations are in play, both solenoidal and dissipative flows characterize the dynamics at steady state.

It is clear from these results that caudal (i.e., posterior) regions have faster intrinsic frequencies, relative to rostral (i.e., anterior) regions. Interestingly, in this example, the inferotemporal and ventral eigenmodes also show a relatively high frequency. At the third scale, this caudal-rostral gradient is more evident, suggesting that faster solenoidal dynamics dominate in posterior parts of the brain. This is consistent with both theoretical and empirical findings suggestive of a gradient of timescales—as one moves from the back to the front of the brain and, implicitly, from hierarchically lower areas to higher areas (Cocchi et al., 2016; Hasson, Yang, Vallines, Heeger, & Rubin, 2008; Kiebel, Daunizeau, & Friston, 2008; Liegeois et al., 2019; Murray et al., 2014; Wang & Kennedy, 2016). Note that the frequencies in question here are very slow; namely, about 0.01 Hz or below. These are the ultraslow frequencies typically characterized in resting-state fMRI (Liegeois et al., 2019). In the present setting, these ultraslow frequencies are an emergent property of scale-invariant behavior, when one moves from spatial temporal scales suitable for describing lobar dynamics or large intrinsic brain networks.

The eigenvalues in Figure 10 take positive real values at higher (second and third) scales. This means that they have crossed the zero threshold to engender a transcritical bifurcation. Strictly speaking, these produce solutions that cannot be realized because of an exponential divergence of trajectories. This reflects the first-order approximation that we are using to summarize the dynamics. Although this linear approximation precludes stochastic chaos, positive real values speak to the notion that some particles at higher scales become excitable for short periods of time. This means that we are moving away from a loosely coupled oscillator model—that has a fixed point or limit cycle attractor—towards what physicists call active or excitable matter (Keber et al., 2014; Ramaswamy, 2010). This is a nice metaphor for the brain and means that if the particles that constitute active (gray) matter are considered in isolation, they can show autonomous dynamics that can be characterized as stochastic chaos or itinerancy.

This scale-free behavior means that we can evaluate the time constants at scales that we have not characterized empirically. Table 1 lists these extrapolated or projected timescales right down to the nanoscale and up to higher scales that would be appropriate to talk about networks of brains or social communities or institutions. Note that the extrapolation of scaling to conferences (Table 1) should not be taken too seriously. There would be good reason to estimate the scaling exponent from a larger database, since small errors in the estimation of scaling exponents can amplify quickly in extrapolation.

It may help to distinguish between scalable and scale-free systems. A scalable system is one in which performing a scale transformation to a system’s state gives a new state, but crucially, one that is also a solution of the governing equation of motion. This is what Landau refers to as “mechanical similarity” in the context of Lagrangian mechanics. For example, orbital trajectories are scalable, seeing as increasing the scale of the system gives new (larger) orbits with different (slower) periods. Crucially however, these new orbits also follow Newton’s second law regardless of size—hence qualifying as “scalable.” On the other hand, a scale-free system has no characteristic-length scale and images taken at different resolutions are statistically unchanging.

This completes our discussion of scale invariance and associated dynamics, where we have taken a special interest in the temporal scaling behavior that emerges from local connectivity at smaller scales of analysis. In the next section, we turn to the coupling between particles and see what this has to say in terms of how intrinsic brain networks influence each other.

The implications of complex extrinsic coupling can be understood in terms of cross-covariance functions of time that characterize delayed or lagged dependencies. Please note that the cross-covariance functions can be evaluated in a straightforward manner from the complex transfer functions, shown in Figure 10. In other words, they can be derived directly from the Jacobian, under first-order assumptions. For example, Figure 14 characterizes these dependencies between the two eigenstates with the strongest coupling at this (third) scale. The implicit coupling is mediated by the corresponding element of the (complex) Jacobian—circled in red in the upper middle panel. The flanking panels on the left and right show the associated eigenmodes in voxel space. The middle row shows the autocovariance functions of the two eigenstates, illustrating serial correlations that can last for many seconds. The interesting part of this figure is in the lower panels: These report the cross-covariance function between the two eigenstates, over 256 s (lower left panel) and 32 s (lower right panel), respectively. The key thing to observe here is that the peak cross-covariance emerges at an 8-s lag from the 24th to the third eigenstate. This asymmetrical cross-covariance (and implicitly cross-correlation) function reflects the solenoidal coupling and implicit breaking of detailed balance accommodated by the particular decomposition (see the next section). Note that the (zero lag) correlation is almost zero. This speaks to the potential importance of using cross-covariance functions (or complex cross spectral in frequency space), when characterizing functional connectivity in distributed brain responses (K. J. Friston, Bastos, et al., 2014; Mohanty, Sethares, Nair, & Prabhakaran, 2020). This brief treatment of extrinsic coupling has made much of the complex nature of dynamical coupling and how it manifests in terms of functional connectivity. In the final section, we revisit this kind of coupling in terms of nonequilibrium steady states.

Above, we have referred to solenoidal and dissipative flows, in relation to the complex and real parts of intrinsic eigenvalues—and how they manifest in terms of intrinsic brain networks. This section unpacks this relationship by applying the statistical physics of nonequilibrium steady states to neuronal fluctuations. Our focus will be on the relationship between conventional characterizations of functional connectivity and the more general formulation afforded by a particular decomposition. Specifically, we will see that conventional formulations assume a special case that discounts solenoidal flow—and implicitly assumes neuronal dynamics attain steady state at statistical equilibrium.

In the previous section, we examined extrinsic coupling among particles in terms of their covariance. Here, we return to coupling and dynamics that are intrinsic to a particle; namely, the final particle at the last level. In this example, the particle has eight eigenstates, whose complex eigenvalues imply a loss of detailed balance and implicit steady state that is far from equilibrium. To understand the link between detailed balance and equilibria versus nonequilibrium steady states, it is useful to consider the eigen-decomposition of the final particle in relation to standard analyses of functional connectivity (e.g., singular value decomposition or principal component analysis of covariance matrices). In what follows, we first rehearse the relationship between flow and steady-state distributions over states afforded by the Helmholtz decomposition. We then look at what this implies under the assumption of symmetric coupling—and how this leads to equilibrium mechanics and a simple relationship between the Jacobian and covariances among their respective eigenstates. We then revisit these relationships but replacing solenoidal flow, to clarify the differences between summarizing dynamics in terms of the eigenvectors of the Jacobian and the eigenvectors of the functional connectivity matrix.

In summary, if we assume detailed balance (i.e., discount solenoidal flow), we are assuming an equilibrium steady state of the sort studied in quantum and statistical mechanics. In this special case, there is a direct relationship between the (eigenvectors of) the Jacobian and the Hessian matrix (i.e., precision, or inverse functional connectivity) matrix. Furthermore, there is also a direct relationship between the Jacobian and the variance of the expected flow or dynamics. The assumption of detailed balance is licensed in many situations; in particular, if we are dealing with ensembles of states or particles that are exchangeable (e.g., an idealized gas). This renders the Jacobian symmetrical and ensures detailed balance. The Jacobian is symmetrical because the influence of one particle on a second is the same as the influence of the second on the first. However, this symmetry cannot be assumed in biological systems that break detailed balance, especially the brain. We now rehearse the above analysis by retaining the symmetry breaking, solenoidal flow that underwrites nonequilibrium steady-state dynamics.

In the presence of solenoidal flow, the eigenvectors of the Jacobian and Hessian are no longer the same. So, which is the best summary of dynamics? Clearly, there is no definitive answer to this question; however, if we are interested in relevant quantities “that matter,” we are specifically interested in non-dissipative, slow, unstable dynamics. By construction, this is what the particular decomposition “picks out,” by discarding fast fluctuations at each successive scale. This means that the eigenstates of the final particle should have identified slow, unstable, or critical dynamics. In contrast, had we just taken the principal components of the covariance matrix of the data (i.e., functional connectivity), we may not have identified the slow modes.

Note that when working with eigenstates, the solenoidal terms are encoded by Q, which is a leading diagonal matrix of imaginary values. Similarly, the dissipative terms are encoded by Γ, which is a leading diagonal matrix of real values. In other words, nonequilibrium steady state—as defined by the prevalence of solenoidal flow—manifests as the imaginary parts of the Helmholtz decomposition, when the system is projected onto the eigenvectors of the Jacobian.

Figure 15 shows the dissipative and solenoidal (kinetic) energy of the eigenstates at the final scale. (Note: Kinetic energy is normally positive because it involves a squared quantity—momentum—usually assumed to be purely real; see Equation 20. However, when we relax this assumption and allow complex momenta-like quantities, squaring no longer guarantees positivity. Hence kinetic energy can take negative values [see Figure 15, top left panel] in this complex state space.) The corresponding eigenmodes are shown in the subsequent panels as maximum intensity projections (of their absolute values). In terms of dissipative dynamics, the first eigenmode has the smallest dissipative energy. In other words, it features the slowest, most unstable mode macroscopic (intrinsic) dynamics. Eigenmodes 2 and 3 are a conjugate pair, with complex parts—and, implicitly, a solenoidal contribution to their (kinetic) energy. These nodes are most pronounced in dorsal mid-prefrontal regions. Note that the kinetic energy of the first eigenstate is negative. This may seem counterintuitive; however, it is a simple reflection of the fact that the principal eigenvalue has a real part that is greater than zero. Clearly, the implicit exponential divergence of trajectories cannot persist globally. In a more complete analysis, the (stochastic) trajectory would quickly enter regimes of dissipation, such that the average real part (cf. Lyapunov exponent) was less than zero. One might ask where the dissipative energy comes from. It is effectively driven by intrinsic fluctuations that, at the finest level, include the fluctuations in active states, which play the role of experimental or sensory inputs. This raises an interesting question: At what scale do experimental inputs manifest?

An intuitive way of thinking about the distinction between dissipative and solenoidal dynamics is in terms of the fluctuations in bath water when perturbed (e.g., when the tap or faucet is running), as opposed to the ripples and waves that persist after the perturbation is removed (e.g., when the tap or faucet is turned off). In one case, the water is trying to find its free energy minimum, while in the second case solenoidal, divergence free flow is more like the complicated swinging of a frictionless pendulum that neither consumes nor creates energy. In this view, it becomes interesting to characterize the response of the system to perturbation—here, the exogenous inputs provided by the experimental design. Conceptually, we can regard experimental inputs (such as visual afferents to the lateral geniculate) as (active states of) external particles that influence (but are not influenced by) the sensory states of particles at the finest level. Practically, these experimental inputs were included in the estimation of the coupling parameters that subtend the Jacobian at the finest scale.

Figure 16 characterizes the influence of exogenous, condition-specific effects at different scales in terms of correlations between fluctuations in the eigenstates that can be explained by any of the three experimental inputs (i.e., visual, motion, and attention). This analysis suggests that the effect of exogenous inputs can be detected at all scales. For example, at the third scale, eigenmode 17 shows an extremely significant effect of sensory perturbation, dominated by visual motion. The associated eigenmode picks out primary visual cortex and extrastriate areas, encroaching upon motion-sensitive regions. Note that eigenmode 17 has a negative time constant, and hence the real part of its eigenvalue is positive. This suggests that during part of its orbit, it becomes transiently unstable, while following the temporal structure of exogenous fluctuations.

Figure 17 provides a complementary and revealing perspective on the effects of sensory perturbation. This figure characterizes induced responses in terms of first-order Volterra kernels—that is, impulse response functions—of particles at the first (finest) scale of coarse-graining. These kernels are based upon the Jacobian and quantify the effects of changing an input on each eigenmode over time (based on the parameters mediating the influence of experimental inputs on motion of states at the first scale). The maximum intensity projections on the left report the variance attributable to each input, based upon the sum of squared kernels over time (i.e., the autocovariance function at zero lag under each input). Note that this is a fundamentally different characterization of brain “activation” because it is modeling the variance induced by an input that is distributed in space and time through recurrent coupling among brain regions.

In this example, motion induces responses in visual and extrastriate—presumably, motion-sensitive eigenmodes—while attention has protracted influences on parietal, prefrontal, and medial temporal regions, including the frontal eye fields and intraparietal sulcus. Visual input per se seems to be expressed preferentially in subcortical systems, including the lateral geniculate but also other subcortical and medial temporal regions. In addition, it appears to selectively engage posterior superior temporal regions in the left hemisphere, often associated with biological motion processing. The interesting aspect of this characterization is the protracted nature of the kernels, which decay to small values after 100 s or so. In effect, this means that although induced responses may be expressed in a regionally specific way almost instantaneously, there are enduring effects that can last for a minute or so, following any exogenous perturbation. Clearly, these effects will be overwritten by ongoing sensory input; however, this suggests that brain systems—and accompanying distributed neuronal responses—have a greater memory than might have been anticipated. Heuristically, this means that I should be able to tell you whether you have “seen something” in the past minute or so by examining your brain activity at this moment in time.

In summary, we have introduced a particular partition that plays the role of a functional parcellation scheme for a system of loosely coupled nonlinear oscillators, such as neuronal populations in the brain. The key aspect of this parcellation scheme is that it can be applied recursively in the spirit of the renormalization group. This enables one to examine the scale-invariant behavior of the ensuing spatiotemporal dynamics in a principled way. The numerical analyses above confirm the analytic intuitions that as we move from one scale to the next, there is a progressive slowing and loss of stability of eigenstates associated with each parcel or particle. This manifests as a form of self-organized criticality, in the sense that slow unstable (non-dissipative) eigenmodes supervene on lower scales. Quantitatively speaking, the spatial scale of a particle, its characteristic frequencies and Lyapunov exponents, all fit nicely with empirical observations of (dynamic) functional connectivity within and among large-scale intrinsic brain networks (Liegeois et al., 2019; Northoff, Wainio-Theberge, & Evers, 2019). The use of an eigen-decomposition of this sort is particularly interesting in relation to recent eigenmode-decompositions of structural connectivity (Atasoy, Donnelly, & Pearson, 2016), cortical geometry (Tokariev et al., 2019), and spatially embedded neural fields (Robinson et al., 2016). These related applications are slightly simpler than the current analysis because they can assume symmetric coupling—of one sort or another—and therefore need only deal with real variables and symmetric modes. For a complementary approach to coarse-graining fMRI data, and then plotting cross-scale covariance functions (including time-asymmetric decompositions), please see Breakspear, Bullmore, Aquino, Das, and Williams (2006), especially Figures 10 and 11. It is also worth noting that our use of the graph Laplacian (G) to define neighboring (i.e., coupled) internal states bears a similarity to the graph signal processing (GSP) scheme (Huang et al., 2018). GSP is used to analyze and integrate structural and functional connectomes, where the (functional) brain activity at the graph nodes is studied on the underlying graph’s (anatomical) structure. In the GSP framework, the eigen-decomposition of graph Laplacian (constructed from structural connectome) is used to identify eigenmodes of low and high frequency (together they define graph Fourier transform).

Because this paper is a technical (foundational) description of the procedures entailed by the existence of Markov blankets, we have focused on the simplest implementation. This means that we started off with linearization assumptions, and propagated this approximation to higher levels. Clearly, it would be nice to revisit the particular partition using higher order approximations that retain nonlinearity in the equations of motion; for example, Equation 14. This would require a more careful analysis of the Lyapunov exponents, which would involve integrating the system and averaging the eigenvalues over the ensuing state-dependent Jacobian: The Jacobian becomes a function of states when one includes nonlinearities in the equations of motion. This raises the interesting issue of how to identify the adjacency matrix used to define the Markov blankets. In other words, we need to establish the conditional independences in terms of a zero entry in the Jacobian. However, if the Jacobian is fluctuating over time, over an orbit in state space, then there may be times when the Jacobian element is zero (i.e., zero coupling) and nonzero at other times. Related numerical analyses of nonlinear systems (K. Friston, 2013) usually require that the Jacobian is zero over a suitably long period of time, when forming the adjacency matrix in Figure 2. Clearly, this would involve evaluating the Jacobian over all the solutions to the trajectory in state space. This may be a time-consuming but otherwise interesting exercise.

We have already mentioned some limitations and extensions. These include starting off with multivariate characterizations of intrinsic dynamics at the finest level. As noted above, this is easy to implement by using the first few principal eigenstates, following a singular decomposition of the smallest particles. Another extension is repeating the dynamic causal modeling at each scale, to reevaluate the Jacobian with suitable high-order (i.e., nonlinear) approximations to the equations of motion.

The nonuniqueness of the particular partition is a key practical issue. There is no pretense that there is any unique particular partition. There are a vast number of particular partitions for any given coupled dynamical system. In other words, by simply starting with different internal states —or indeed the number of internal states per particle—we would get a different particular partition. Furthermore, the thresholds used in the elimination of fast dissipative eigenmodes will also change the nature of the partition, leading to more or less inclusive dynamics at the scales above. This latter aspect is probably more defensible in terms of summarizing multiscale behavior, in the sense that we can easily motivate the adiabatic approximation in terms of the relative stability of eigenmodes at a particular level. However, the number of internal states to consider—and how to pick them—introduces a more severe form of nonuniqueness. In this paper, we used the state that was maximally coupled to other states as the internal state of the next particle. This was based upon the graph Laplacian of the adjacency matrix at the appropriate scale. This is a sensible but somewhat arbitrary definition of an internal state and speaks to the point that there are a multitude of particular partitions—and implicit Markov blankets—that could be used. There are two ways that one could handle this nonuniqueness. (Note: As noted by one of our reviewers, it might be useful to fix the [initial] centroids of the particles to the centroids of brain atlases [for example, Glasser et al., 2016; Schaefer et al., 2018]. This might ease the interpretation of the findings, as well as comparisons between the current and complementary brain network studies.) One would be to embrace it and focus on the statistics of characterizations over different particular partitions and look for scaling behaviors that are conserved over partitions. The alternative is to think about a unique particular partition and how this would be identified. This as an outstanding issue; namely, what is the “best” particular partition and, indeed, is the notion of the best partition appropriate?

Another important caveat is the fact that we have predicated the illustrative analyses in this paper on a single-subject dataset acquired under an experimental activation paradigm. We chose this dataset because it has been used to illustrate previous developments of dynamic causal modeling. Conceptually, this means that the particular partition is specific to this subject and the subject’s responses to the attentional paradigm (summarized in Figure 17). Because this paradigm introduced context-sensitive or condition-specific changes in effective connectivity, it was designed to change the Jacobian over different periods of stimulation (e.g., attentional modulation of coupling between visual motion areas and early visual cortex). We did not attempt to model these effects here; this would require the nonlinear modeling mentioned above. If this modeling was to second order, we would end up with a bilinear form for Equation 2, which is the basis of most DCM analyses of fMRI data. This speaks to the fact that the parcellation scheme may not produce the same results when applied to a different paradigm. In turn, this means that there is further work to be done in terms of finding a particular partition that accommodates variability in functional anatomy. In theory, this would probably be best addressed using a generative model; in other words, assuming one underlying sparse Jacobian at any given scale and then adding random effects, so that it could be used to explain multiple paradigms or subjects. Having said this, the current analyses can be taken as proof of principle that this sort of multiscale decomposition can be applied to empirical neuroimaging time series and leads to the same phenomenology reported in functional connectivity literature.

The software producing the figures in this figure are available as part of the academic software SPM. They can be accessed by invoking DEM graphical user interface and selecting the DCM and blankets button (DEMO_DCM_MB.m). Please see https://www.fil.ion.ucl.ac.uk/spm/.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00175.

Karl J. Friston: Conceptualization; Methodology; Visualization; Writing – original draft. Erik D. Fagerholm: Methodology; Writing – review & editing. Tahereh S. Zarghami: Writing – review & editing. Thomas Parr: Methodology; Writing – review & editing. Inês Hipólito: Writing – review & editing. Loïc Magrou: Writing – review & editing. Adeel Razi: Conceptualization; Methodology; Writing – review & editing.

Karl J. Friston, Wellcome Trust (http://dx.doi.org/10.13039/100004440), Award ID: 088130/Z/09/Z. Erik D. Fagerholm, King’s College London Prize Fellowship. Tahereh S. Zarghami, Cognitive Sciences and Technologies Council of Iran international research visit. Adeel Razi, Australian Research Council, Award ID: DE170100128. Adeel Razi, Australian Research Council, Award ID: DP200100757. Adeel Razi, National Health and Medical Research Council, Award ID: 1194910.

Functional connectivity:

A (undirected) measure of the statistical dependencies between spatially remote neurophysiological events.

Effective connectivity:

A measure of the directed (causal) influence of one neural system over another using a model of neuronal interactions

Dynamic causal modeling:

A Bayesian framework that is used to infer causal interaction between coupled or distributed neuronal systems.

Markov blanket:

A Markov blanket allows one to distinguish a collection of states that belong to a particle from states that do not.

Jacobian:

A matrix that contains a first-order partial derivative for a vector function.

Lyapunov exponent:

It gives the rate of exponential divergence from perturbed initial conditions.

Bayesian model reduction:

A Bayesian inversion and comparison of models that are reduced (or sparsed) forms of a full (or parent) model.

Eigenmode:

A stable state (i.e., mode) of a dynamical system in which all parts of the system oscillate at the same frequency.

Graph Laplacian:

A matrix representation of graph that combines node adjacency and node degree in mathematical formulation and belongs to spectral graph theory.

Generative model:

A model for randomly generating observable data values, typically given some hidden parameters.