Entropy, free energy, symmetry and dynamics in the brain

Predictive coding is essentially based on three equations central to its discussion and associated concepts [9, 12–14].

$\begin{equation*}p(x,y)=p(y\vert x)p(x)\end{equation*}$

$\begin{equation*}\dot{Q}=f\left(Q,k\right)+\mathrm{v}\end{equation*}$

$\begin{equation*}Z=h\left(Q\right)+\mathrm{w}.\end{equation*}$

The first equation establishes a reduced form of the Bayesian theorem, expressed in terms of probability distribution functions p, where p(x, y) is the joint probability of the two state variables x and y, and p(y|x) is the conditional probability of the variable y given the state of the variable x. In the Bayesian framework, parameters and state variables have a similar statue in the sense that they can each be described by distributions and enter as arguments into the probability function p. For instance, the probability to jointly find states x and y for a given parameter k is written as p(x, y|k), establishing the likelihood of obtaining a set of data x, y, given a set of parameter values k, which have a prior distribution of p(k). The prior represents our knowledge about the model and the initial values.

The second equation is known as the Langevin equation ¹ and establishes the generative model, in which brain activations at the neural source level are represented by the N-dimensional state vector$Q=(x,y,\dots )\in {\mathrm{R}}^{N}$, $f\left(Q,k\right)$ are the deterministic influences expressed as an M-dimensional flow vector f depending on the state Q and the parameter k or a set of multiple parameters $\left\{k\right\}$. $\mathrm{v}\in {\mathrm{R}}^{N}$ establishes the fluctuating forces, typically assumed to be Gaussian white noise with ⟨v_i (t)v_j (t')⟩ = cδ_ij δ(t − t'), where δ_ij is the Kronecker-delta and δ(t − t') the Dirac function. More general formulations of the influence of noise including multiplicative or colored noise are possible and we here refer the reader to the relevant literature.

The third equation establishes the observer model and links the source activity Q(t) to the experimentally accessible sensor signals Z(t) via the forward model h(Q) and measurement noise w. For electro-encephalographic measurements, h is the gain matrix established by the Maxwell equations; for functional magnetic resonance imaging measurements, h is given by the neurovascular coupling and the hemodynamic Ballon–Windkessel model. In the present context, the observer model is of no concern, although it is of immense importance in real world applications and plays often the role of a major contaminating factor in problems of model inversion and parameter estimation. We mention these engineering issues here for completeness, but will now for simplicity assume h to be the identity operation with zero measurement noise, thus Z = Q.

Predictive coding relates to a large field of research in behavioral neurosciences [21, 22], in particular ecological psychology [23, 24] devoted to the scientific study of perception-action and dynamical systems. Here James J Gibson stressed the importance of the environment [25], in particular, the perception of how the environment of an organism affords various actions to the organism. This loop of perception-action closely matches the loop of predictions from and update of the internal generative model in predictive coding. A particular nuance emphasized in ecological psychology has been the ecologically available information—as opposed to peripheral or internal sensations—leading to the emergence of the perception-action dynamics. Scott Kelso and colleagues have made major contributions to the formalization of this framework and developed experimental paradigms to test the dynamical properties of the internal models underlying the coordination between the organism and the environment [26–28]. These paradigms were theoretically inspired by Hermann Haken's synergetics [4], which has been foundational for the theories of self-organization. It generated intense research efforts with a strong focus on transitions between states in perception and action, including modeling [29–32] and systematic experimental testing in bimanual and multisensorimotor coordination [33–44]. The approaches were then generalized to larger range of paradigms [45–53] with the goal to extract the principle features of the organism's behavior in interaction with the environment [22]. This large body of work generated substantial evidence for the benefit of dynamical descriptions of behavior, which can be subsumed in the framework of SFMs for functional perception-action variables in behavior (see [48] for a recent review) and the brain [54–56].

Our understanding of the deterministic and stochastic processes determining the discrete value x_i of the state variable x is captured by the computation of the corresponding probability distributions p_i . Shannon demonstrated the remarkable fact that there exists an information quantity H(p₁, ..., p_n ), which uniquely measures the amount of uncertainty represented by these probability distributions [8]. In his original proof, he showed that the requirement of three basic conditions, most notably the composition law linking events and probabilities, leads naturally to the following expression

$\begin{equation*}H({p}_{1},\dots ,{p}_{n})=-K\sum\limits _{i}\enspace {p}_{i}\enspace \mathrm{ln}\enspace {p}_{i},\end{equation*}$

where K is a positive constant. As the measure H corresponds directly to the expression for entropy in statistical mechanics, it is referred to as information entropy [7]. Commonly the discussion of its properties evolves from the perspective of specification of probabilities linked to the amount of information available to an observer. Laplace's principle of insufficient reason assigns equal probabilities to two events in absence of any distinguishing information. This subjective school of thought regards probabilities as expressions of human ignorance, expressing formally our expectation if an event occurs or not. Such thinking is fundamental to theory of predictive coding and its interpretation of cognitive processes generated by the brain. The objective school of thought is rooted in physics and regards probability as an objective property of the event, which can be, in principle, always measured by frequency ratios of events in a random experiment. Here we wish to take an intermediate position in this regard by investigating the deterministic and stochastic processes underlying predictive coding. The reference to deterministic and stochastic forces typically pertains to the language of the objective school, where the subsequent interpretation within the predictive coding framework is part of the subjective school. The generative model in predictive coding represents both types of forces through the Langevin equation and shapes the probability functions, which then provide us with access to this information through the empirical measurement of a function ⟨g(x)⟩. The triangular brackets denote the expectation value

$\begin{equation*}\langle g(x)\rangle =\sum\limits _{i}\enspace {p}_{i}g({x}_{i}).\end{equation*}$

These correlations of g(x), together with the normalization requirement, ${\sum }_{i}\enspace {p}_{i}=1$, express the constraints provided by deterministic and stochastic influences, within which the information entropy has a maximum. Any other assignment than the maximum of information entropy would introduce another deterministic bias or arbitrary assumption, which by hypothesis we do not have. This insight is the essence of Edwin T Jaynes' maximum information principle [7]. Consequently, we maximize the information entropy H under these constraints by introducing Lagrange parameters λ_i in the usual way known from classical mechanics and obtain

$\begin{equation*}{p}_{i}=\mathrm{exp}\left(-{\lambda }_{0}-\sum\limits _{j}\enspace {\lambda }_{j}{g}_{j}({x}_{i})\right),\end{equation*}$

where the expression for the probability distribution is generalized for multiple measured functions g_j (x_i ) of the event x_i . The entropy S of this distribution is maximal and reads

$\begin{equation*}{S}_{\mathrm{max}}={\lambda }_{0}+\sum\limits _{j}\enspace {\lambda }_{j}\langle {g}_{j}(x)\rangle .\end{equation*}$

The link between empirical measurements and the maximum information entropy is established by the partition function Z (as in Zustandssumme)

$\begin{equation*}Z({\lambda }_{1}\dots {\lambda }_{n})=\sum\limits _{i}\mathrm{exp}\left(-\sum\limits _{j}\enspace {\lambda }_{j}{g}_{j}({x}_{i})\right)\end{equation*}$

and the relations

$\begin{equation*}\langle {g}_{j}\left(x\right)\rangle =-\frac{\partial }{\partial {\lambda }_{j}}\enspace \mathrm{ln}\enspace Z\end{equation*}$

$\begin{equation*}{\lambda }_{0}=\mathrm{ln}\enspace Z.\end{equation*}$

The correlations captured by the empirical measurements ⟨g_j (x)⟩ link through these equations to the stationary probability distribution functions and consequently to the generative model (aka deterministic forces) expressed through the Langevin equation. In the following sections, we will provide explicit examples with clear neuroscience relevance.

Emergence in self-organizing systems requires the change of stability of some typically low-dimensional attractor. The system under consideration is nonlinear and high dimensional with N degrees of freedom. In the space spanned by these degrees of freedom, each point is a state vector and represents a potential state of the full system. As time evolves, the state of the system changes and thus traces out a trajectory in state space. The rules that the system follows can be understood as forces that cause the changes of the state vector and define a flow. In order to allow this system to generate low-dimensional behavior, that is, M dimensions withM ≪ N, there must be a mechanism in place that is capable of directing the trajectories in the high-dimensional space toward the lower M-dimensional sub-space. Mathematically, this translates into two flow components that are associated with different time scales: first, the low-dimensional attractor space contains a manifold M, which attracts all trajectories on a fast time scale; second, on the manifold a structured flow F(.) prescribes the dynamics on a slow time scale, where here slow is relative to the collapse of the fast dynamics toward and onto the attracting manifold, as illustrated in figure 2. For compactness and clarity, imagine that the state of the system is described by the N-dimensional state vector Q(t) at any given moment in time t. Then we split the full set of state variables into the components q and s, where the state variables in q define the M task-specific variables linked to emergent behavior in a low-dimensional subspace (the functional network) and theN − M variables of s define the remaining recruited degrees of freedom. Naturally, N is much greater than M and the manifold in the subspace of the variables q has to satisfy certain constraints to be locally stable, in which case all the dynamics is attracted thereto.

Zoom In Zoom Out Reset image size

We briefly outline some of the concepts of synergetics first, which will contextualize and aid in better understanding the concepts of SFM. Synergetics is the theory of self-organized pattern formation in open systems (i.e., those that are in contact with the environment through matter, energy, and/or information fluxes) far from thermo-equilibrium that are composed of numerously weakly interacting microscopic elements [4]. Due to the interaction among the microscopic elements, such systems may become organized in spatially as well as temporally ordered patterns that are typically macroscopic in nature and can be described by a limited number of so-called order parameters (or collective variables). Spontaneous switching of brain activity patterns (i.e., non-equilibrium phase transitions) occurs as one macroscopic pattern losses stability while another stable pattern takes hold. The stability of a system's state (or phase) implies that, if perturbed away from it, the system will tend to return to that state. When stability is lost, the system will instead tend to move away from that state and rather switch to another stable state. Close to such points of (macroscopic) instability the time to return increases tremendously. As a result, the macroscopic state evolves rather slowly in response to perturbation, whereas the underlying microscopic components maintain their individual time scale. Consequently, the time scales of their dynamics differ tremendously (time-scales separation). From the perspective of the slowly evolving macroscopic state, the microscopic components change so quickly that they can adapt instantaneously to macroscopic changes. Thus, even though the macroscopic patterns are generated by the subsystems, the former, metaphorically, enslaves the latter [4]. Ordered states can always be described by a very few variables (at least in the close vicinity of a bifurcation) and consequently, the state of the originally high-dimensional system can be summarized by a few or even a single collective variable, the order parameter(s). The order parameters then span the workspace. The circular relationship between the enslaving order parameters and the enslaved microscopic components, which generate the order parameters, is sometimes referred to as circular causality, which effectively allows for a low-dimensional description of the dynamical properties of the system [57]. The notion of circular causality and the emergence of low-dimensional dynamics in self-organizing systems are at the heart of Hermann Haken's synergetics [4].

The mathematical formalism of synergetics can be briefly sketched as follows. An N-dimensional state vector $Q(t)\in {\mathrm{R}}^{N}$ is defined as a function of time and comprises all state variables of the system. The evolution of the state vector is captured by a nonlinear ordinary differential equation:

$\begin{equation}\dot{Q}=F(Q,\left\{k\right\})+\upsilon (t),\end{equation} \tag{ 1 }$

where F is a nonlinear function capturing all possible interactions and $\left\{k\right\}$ is the set of control parameters, which control the state of the system and are time-independent. The nonlinear evolution equation includes the influence of noise υ(t), which shall be considered later, but not here. The state vector and its evolution equations can also be extended to incorporate spatial dependence of the state variables, but it suffices here to consider only spatially discrete systems. For a detailed mathematical treatment see Haken (1983).

The entry point to self-organization and emergence in synergetics is the qualitative change (state transitions) of states of the system on the macroscopic level. A given state typically corresponds to a (stationary) fixed point Q₀ or periodic (limit) cycle in state space (see figure 2) and synergetics aims at describing the flow in its local vicinity in state space (figure 3).

Zoom In Zoom Out Reset image size

The deviation ɛ(t) from the steady state Q₀ is Q(t) = Q₀ + ɛ(t) and its time evolution can be approximated by a Taylor expansion of F around the steady state Q₀

$\begin{equation*}\dot{Q}=F(Q,\left\{k\right\})=F({Q}_{0},\left\{k\right\})+L({Q}_{0},\left\{k\right\})\varepsilon +A({Q}_{0},\left\{k\right\})\varepsilon :\varepsilon +\cdots ,\end{equation*}$

where $L({Q}_{0},\left\{k\right\})$ and $A({Q}_{0},\left\{k\right\})$ are the first and second order terms of the Taylor expansion. As the control parameters $\left\{k\right\}$ are varied, the system is guided through parameter space and will encounter parameter configurations, for which bifurcations occur, leading to the change of stability of the current state and a major reorganisation of the system on the macroscopic level. The bifurcation point defines the working point k₀ and hence a second local constraint, this time in parameter space, next to the working point Q₀ in state space. Under these local conditions, at least one of the eigenvalues λ₀ of the Jacobian $L({Q}_{0},\left\{{k}_{0}\right\})$ becomes zero ² and allows the application of the local center manifold theorem. It states that a small subset of new variables $q(t)\in {\mathrm{R}}^{M},M\ll N$, the order parameters, arises, which dominate the system dynamics and enslave the remaining variables $s(t)\in {\mathrm{R}}^{N-M}$ by means of time scale separation, where the order parameters q(t) are slow and the enslaved variables s(t) are fast. In the jargon of synergetics, enslaving means that the intrinsic dynamics of the variables s(t) can be adiabatically eliminated and its long term evolution is expressed as a function of the order parameters, that is s(t) = s(q(t)). This leads then to the reduced description of the emergent dynamics as

$\begin{equation*}\dot{q}={\lambda }_{0}q+P\left(q,s\left(q\right),\left\{{k}_{0}\right\},{Q}_{0}\right),\quad q(t)\in {\mathrm{R}}^{M},M\ll N\end{equation*}$

$\begin{equation*}s(t)=s(q(t)),s(t)\in {\mathrm{R}}^{N-M},\end{equation*}$

where P is a nonlinear function parametrized by $\left\{{k}_{0}\right\},{Q}_{0}$ and can be analytically computed from F by use of the decomposition into order parameters and enslaved variables.

Although synergetics is conceptually not limited to local working points in state space, practically this has always been the case. In biology this has placed an enormous limitation on the utility of the synergetic framework, as the exploration of the flow in state space appears to be a fundamental activity of the organism. To this end, symmetry offers itself as another guiding principle to define the working point in parameter space [9, 16], with the objective to overcome the limitation of confinement to local regions in state space. Dynamical systems with symmetries are called equivariant dynamical systems [59–62], for which group theoretical methods provide a natural language [59, 60]. To demonstrate this, let us consider equation (1) with the state variables formally expressed in $q(t)\in {\mathrm{R}}^{M},M\ll N$ rather than Q(t) $\in {\mathrm{R}}^{N}$, for reasons. We return to this point at the end of section 4.2. Further, for simplicity, let us consider only one control parameter k. Let then Γ be a group acting on its solutions. The equation is said to be Γ-equivariant if F commutes with the group action of Γ, that is $F(\gamma \cdot q)=\hat{\gamma }\cdot F(q)$ for all γ ∈ Γ. This symmetry shall exist for a critical control parameter value k₀. An important consequence of Γ-equivariance is that if a solution q(t) solves the ordinary differential equation, then so does γ ⋅ q(t) for all γ ∈ Γ. If the symmetry is continuous, that is γ ⋅ q = q + δq, then the corresponding group is a Lie group, whose elements have the topology of an M-dimensional smooth manifold M in state space and whose group operation is a smooth function of the elements. Then the stationary solutions of equation (1) span a smooth manifold M, which is defined by

$\begin{equation*}\dot{q}=F(\gamma \cdot {q}_{0})=F({q}_{0}+\delta q)=\hat{\gamma }\cdot F({q}_{0})=0\end{equation*}$

allowing a continuous displacement δq along the manifold. In the case of symmetry breaking, k = k₀ + μ, where μ is small, the solutions of the system can be approximated by perturbations of the fully symmetric solution

$\begin{equation*}\dot{q}=F(q,\left\{k\right\})=\underset{=0}{\underbrace{F({q}_{0},{k}_{0})}}+\underset{M}{\underbrace{\frac{\partial F}{\partial q}(q,k){\left\vert \right.}_{{q}_{0},{k}_{0}}}}\cdot (q-{q}_{0})+\underset{N}{\underbrace{\frac{\partial F}{\partial k}(q,k){\left\vert \right.}_{{q}_{0},{k}_{0}}}}\cdot \underset{\mu }{\underbrace{(k-{k}_{0})}}+\cdots ,\end{equation*}$

where M is the smooth invariant manifold of stationary solutions. A zero eigenvalue is associated with the tangent space to M which establishes a time scale hierarchy to the other degrees of freedom as is typically known from synergetics. For full symmetry, k = k₀, all points on the manifold are stable fixed points if M is stable. For small symmetry breaking, μ ≪ 1, a slow flow emerges along the manifold, which is slow as compared to the fast dynamics orthogonal to the manifold. These two cases are illustrated in figure 4 for a circular manifold M : 0 = 1 − x² − y².

Zoom In Zoom Out Reset image size

The structure of the flow on the manifold is entirely determined by the symmetry breaking through k in N. In the neuroscience context, e.g. in large-scale brain network models, such a symmetry breaking result from variations in the connectivity between brain regions or the local properties of the individual regions. There is no associated immense dimension reduction as we know it from the traditional framework in synergetics, causing a compression from N to M dimensions for the order parameters $q(t)\in {\mathrm{R}}^{M},M\ll N$. Rather the symmetry is obliged to be defined within a subspace spanned $q(t)\in {\mathrm{R}}^{M},M\ll N$, which then is completed by the remaining N − M variables $s(t)\in {\mathrm{R}}^{N-M}$ as in the traditional synergetic framework. Assuming that these variables are not passing through an instability at the critical point, k = k₀, then they can be adiabatically eliminated via enslaving in the usual way and their dynamics is expressed as a function of the order parameters q(t). The difference between the traditional synergetic framework and SFMs is that in the former, the full N-dimensional system is considered allowing the full dimension reduction naturally, whereas in the latter the added constraint needs to be assumed or imposed. These considerations lead us to the following form of system equations:

$\begin{equation*}\dot{q}=M\left(q,s,\left\{{k}_{0}\right\},{Q}_{0}\right)\cdot (q-{q}_{0})+\mu N\left(q,s,\left\{{k}_{0}\right\},{Q}_{0}\right)\end{equation*}$

$\begin{equation*}s(t)=s(q(t)),\end{equation*}$

where $q(t)\in {\mathrm{R}}^{M},M\ll N$, and $s(t)\in {\mathrm{R}}^{N-M}$. This set of equations established the basic mathematical framework expressed by SFMs, in which we will develop the basic network equations next.

Neural mass models are reduced mathematical representations of the collective activity of neurons. They are typically derived from a population of neurons, which are represented by coupled point neuron models. Under assumptions about the statistics of the distribution of action potentials and/or coupling between neurons, mean field theory is applied to derive equations of collective variables, capturing the evolution of mean, variance and higher statistical momenta of the population. Notable examples include the Brunel Wang model assuming Poisson distributed spikes [63], Zerlaut et al model using master equation and transfer function formalisms [64], and the Stefanescu–Jirsa model [65, 66] exploiting the heterogeneity of neuronal parameters leading to synchronized clusters of neurons. The mean field derivation of Montbrio et al is particularly attractive from the theoretical perspective [67], because it is exact under the assumption of Lorentzian ansatz. It derives two collective state variables, the mean firing rate r and the mean neuronal membrane potential V. The corresponding phase flows are plotted in figure 5(A), the equations in figure 5(B). These and all other neural mass models have in common that they reduce the mean field dynamics to a low-dimensional representation, often in 2 dimensions. The neural mass dynamics typically comprises a down state corresponding to a low firing rate, an upstate corresponding to a high firing rate, and the capacity to show oscillatory dynamics prevalently in the upstate. Ignoring the oscillatory dynamics for the time being, we can conceptually reduce the co-existing up and down states and their stability changes via bifurcations to the phaseflow of a single variable, x, as illustrated in figure 5(C). Here the phase flow is reduced to the separation of the two stable states, which may lose stability via saddle-node bifurcations under the change of an external control parameter k. This model is mathematically identical to Wong–Wang model, which was derived under adiabatic approximation from the Brunel Wang model. We wish to emphasize that this representation is not specific to the Brunel Wang neural mass model, but captures the essential dynamic features of all neural mass models, that is co-existence of a low- and high-firing state for a range of intermediate control parameter values, loss of the low/high-firing state via bifurcations at high/low values of the control parameter. For these reasons we here use this reduced neural mass as the basic building block of the equivariant brain network model in the next section.

Zoom In Zoom Out Reset image size

Symmetry breaking in a network composed of bistable neural masses naturally leads to the creation of SFMs [16]. This can be understood as follows. Let us first consider an intuitive toy network of two nodes with state variables x and y. The underlying equations read

$\begin{equation*}\dot{x}=f\left(x,k\right)+\mathrm{v}=x(1-{x}^{2})+k+\mathrm{v}\end{equation*}$

$\begin{equation*}\dot{y}=f\left(y,k\right)+\mathrm{v}=y(1-{y}^{2})+k+\mathrm{v}\end{equation*}$

and form a system of coupled neural mass models, where k is the local excitability and v the noise, following the notation of equation (1). Figure 6 shows the phase flow in state space for this situation, illustrating four stable fixed points, one unstable fixed point and four saddle points. The red and the green lines identify the nullclines. The above equations can now be formally rewritten as

$\begin{equation*}\dot{x}={f}_{x}\left(x,y,G,k\right)+\mathrm{v}=x(1-{x}^{2}-{y}^{2})+G{y}^{2}x+k+\mathrm{v}\end{equation*}$

$\begin{equation*}\dot{y}={f}_{y}\left(x,y,G,k\right)+\mathrm{v}=y(1-{x}^{2}-{y}^{2})+G{x}^{2}y+k+\mathrm{v},\end{equation*}$

where the first term on the rhs represents the same circular invariant manifold as seen in the previous section 3.3. For G = 1, these equations are identical to the uncoupled system. As G varies from 1 to 0, the nullclines change their form continuously from infinity (straight lines) through elliptic to a closed circular shape (see figure 6). G acts as a second control parameter quantifying the degree of mutual connectivity, ranging from fully disconnected (G = 1) to all-to-all coupling (G = 0) topology. For a highly interconnected network, that is G ≈ 0, these intermediate values constrain the phase flow on a closed manifold, centered around the origin, creating a symmetric flow and time scale separation. If other forms of symmetry breaking are introduced, such as regional variance of k or asymmetric connections, then this provides a means of systematically controlling the flow on the manifold.

Zoom In Zoom Out Reset image size

Extending this network to three network nodes, the equations read

$\begin{equation*}\dot{x}=f\left(x,y,z,G,k\right)+\mathrm{v}=x(1-{x}^{2}-{y}^{2}-{z}^{2})+G({y}^{2}+{z}^{2})x+k+\mathrm{v}\end{equation*}$

$\begin{equation*}\dot{y}=f\left(x,y,z,G,k\right)+\mathrm{v}=y(1-{x}^{2}-{y}^{2}-{z}^{2})+G({x}^{2}+{z}^{2})y+k+\mathrm{v}\end{equation*}$

$\begin{equation*}\dot{z}=f\left(x,y,z,G,k\right)+\mathrm{v}=y(1-{x}^{2}-{y}^{2}-{z}^{2})+G({x}^{2}+{y}^{2})z+k+\mathrm{v},\end{equation*}$

where the corresponding phase flow is plotted for G = 0 in state space shown in figure 7. This situation corresponds to the unconnected network nodes as in the previous example. Again, the origin is an unstable fixed point, located symmetrically within a cube comprising 1 unstable fixed point, 8 stable fixed points, separated by 12 saddle points, that are aligned along the nullclines defining the edges of the cube (see figure 7). As connectivity is being established between the three nodes, G reduces to smaller values toward the fully connected network with G = 0. For the latter, the invariant manifold is a 2-sphere, centered at the origin with radius 1 and zero flow.

Zoom In Zoom Out Reset image size

This situation can be formally generalized to a network composed of N basic units

$\begin{equation*}{\dot{x}}_{i}=f\left({x}_{i},{x}_{j},G,{k}_{i}\right)+{\mathrm{v}}_{i}={x}_{i}\left(1-\sum\limits _{j=1}^{N}{{x}_{j}}^{2}\right)+G\sum\limits _{j\ne i}^{N}{c}_{ij}{{x}_{j}}^{2}{x}_{i}+{k}_{i}+{\mathrm{v}}_{i}.\end{equation*}$

The network establishes an equivariant matrix in N-dimensional state space, in which for the absence of any coupling, 2^N stable fixed points exist, separated by the same number of unstable fixed points, all centered around the origin. As connections are established, the system's phase flow is more and more constrained to the proximity of the (N − 1)-sphere ${S}^{N-1}=\left\{x\in {\mathrm{R}}^{N}:{\Vert}x{\Vert}=1\right\}$, which is the manifold M as discussed in section 2.2 and composed entirely of stable fixed points for G = 0, c_ij = 1 ∀ i, j, in full analogy to the 2 and 3 dimensional cases. In case of more elaborate symmetry breaking, for instance through the introduction of a connectome, here c_ij ≠ 1, or variation in regional excitability k_i and local noise v_i , a large range of structured flows can be established on this manifold.

Connectome based network modeling has been performed extensively for resting state activity [68–75], but no rigorous theoretical perspective has yet been offered beyond descriptive statistics (such as functional connectivity, functional connectivity dynamics, multiscale entropy) and hypothesized mechanisms (subcriticality, stochastic resonance). The symmetry breaking of the equivariant matrix through the connectome adds an attractive alternative explanation. However, it still lacks an important explanatory argument, as there is no substantial dimension reduction provided by the symmetry breaking of the equivariant matrix. The invariant manifold remains N − 1 dimensional. As discussed in section 3.3 in general, and now here specifically for neuroscience, it poses questions in terms of how such a dimension reduction of the variables Q(t) inN-dimensional space to variables q(t) in M-dimensional space can be accomplished, as M should be much smaller than N.

The following is a possible realisation. In more realistic neuroscience models, the fixed points will not be identical. Generally the down state is more stable than the upstate. Although such differences will inevitably introduce gradients into the flow of the network, it will not suffice to reduce the low-dimensional manifolds that are required for task-specific processes and are known to exist in the human brain activity. Discussions in the literature evoke the possibility of decorrelation as an important mechanism for information processing [76–78] and may here be realized by oscillations in the upstate. Decoupling via frequency separation or averaging is a well known mechanism in dynamical system theory (e.g. rotating wave approximation) and exploited in large-scale brain networks for the organization of synchrony [79] and cross-frequency coupling [80]. Such capacity for oscillatory behavior is an excellent candidate for the purpose of disconnecting an equivariant M-dimensional manifold from its N − M-dimensional complementary subspace with M ≪ N, offering a clear path forward for future scientific investigation.

The free energy can be seen to mirror the evolution of structured flows on the manifold created by the deterministic features present in the system. In the context of the large scale brain network models, we demonstrated that the symmetry breaking via the connectome qualifies as a candidate mechanism underlying the emergence of these features, thereby establishing the link to biophysical processes such as Hebbian learning and other plasticity mechanisms. The system is then entropically driven by stochastic fluctuations exploring the manifold. Each point on the manifold is associated with a probability distribution. As a limit case in absence of stochastic forces, these distributions are sharply peaked approximating delta-functions and the structured flow is fully deterministic [16]. The link between these deterministic and stochastic influences is the Fokker–Planck equation, which prescribes the temporal evolution of the probability distribution functions [4]. By restricting our attention to the prediction of its stationary properties, the interpretation of the statistical properties is rendered time independent, and we can refer to states, otherwise the interpretation of the state of the system at time t is to be carried out on the basis of measurements performed at time t only.

We wish to illustrate this along a two examples. The first example pertains to the situation previously discussed in section 4.1, in which two nodes were coupled. Let us place the initial focus on only one of the two stable fixed points. We linearize its flow and obtain the following expression

$\begin{equation*}\dot{x}={f}_{1}(x,y)=-x+2\beta y+\mathrm{v}=-\frac{\partial V}{\partial x}+\mathrm{v}\end{equation*}$

$\begin{equation*}\dot{y}={f}_{2}(x,y)=-(y-{y}_{0})+2\beta x+\mathrm{v}=-\frac{\partial V}{\partial y}+\mathrm{v},\end{equation*}$

where V(x, y) is the potential, v the noise, y₀ a deplacement and β the coupling strength. The deterministic and stochastic influences express themselves in the experimentally accessible correlations, establishing the constraints that need to be satisfied when constructing the probability distribution function p(x, y). The expectation value of a function g(x, y) is ⟨g(x, y)⟩ which can be widely represented by the statistical momenta ⟨x⟩, ⟨y⟩, ⟨x²⟩, ⟨y²⟩, ⟨xy⟩, .... The stationary probability distribution function

$\begin{equation*}p(x,y)=N{\mathrm{e}}^{-2V(x,y)/Q}\end{equation*}$

is the time independent solution to the Fokker–Planck equation ³

$\begin{equation*}\dot{p}+\sum\limits _{i}\frac{\partial }{\partial {x}_{i}}\left({f}_{i}(x,y)-\frac{1}{2}Q\sum\limits _{i}\frac{\partial p}{\partial {x}_{i}}\right)=0,\end{equation*}$

where x₁ = x and x₂ = y, and N is the normalization constant. It can be expressed by the ansatz

$\begin{equation*}p(x,y)=\mathrm{exp}({\lambda }_{0}+{\lambda }_{1}{x}^{2}+{\lambda }_{2}x+{\lambda }_{3}xy+{\lambda }_{4}y+{\lambda }_{5}{y}^{2}),\end{equation*}$

where all the Lagrange multipliers λ_i can be estimated explicitly in the usual way from the experimentally accessible correlations ⟨g(x, y)⟩ and the normalization condition ${\int }_{-\infty }^{\infty }\enspace p(x,y)\mathrm{d}x\enspace \mathrm{d}y=1$. Without limiting the generality, we can write the above for simplicity assuming λ₁ = −1, λ₂ = 0 and λ₃ = −2β, λ₄ = 2y₀, λ₅ = −1 and obtain

$\begin{equation*}p(x,y)=C\enspace \mathrm{e}\mathrm{x}\mathrm{p}(-{x}^{2}+2\beta xy+2{y}_{0}y-{y}^{2}).\end{equation*}$

The normalization factor is obtained from the term with λ₀. This joint probability distribution function can be formally rewritten to reflect the Bayes theorem and link back to our initial discussion, that is

$\begin{equation*}p(x,y)=p(x\vert y)p(y)={N}_{x}\enspace \mathrm{exp}(-{(x-\beta y)}^{2}/Q){N}_{y}\enspace \mathrm{exp}(-((1-{\beta }^{2}){y}^{2}+2{y}_{0}y-{y}_{0}^{2})/Q),\end{equation*}$

where N_x , N_y are the normalization constants of p(x|y), p(y). If the two nodes are independent, then naturally the momenta factorize ⟨xy⟩ = ⟨x⟩⟨y⟩ and the deterministic coupling β = 0. The conditional probability p(x|y) becomes independent of y, that is p(x|y) = p(x), and both p(x) and p(y) are fully Gaussian. The expressionF = (x − βy)²/Q is the free energy and intuitively captures the interaction between the nodes as a deformation of the stationary probability density as illustrated in figure 8.

Zoom In Zoom Out Reset image size

The second illustrative example comprises the equivariant matrix in section 4.2. For simplicity, we restrain the discussion again to two nodes. Following the same mathematical steps as in the previous example, the potential function for the coupled nodes reads

$\begin{equation*}V(x,y)=-\frac{1}{2}{y}^{2}+\frac{1}{4}{y}^{4}-2ky-\frac{1}{2}{x}^{2}+\frac{1}{4}{x}^{4}-2kx+\frac{1}{2}\beta {x}^{2}{y}^{2},\end{equation*}$

where β = G − 1 is the coupling strength, with β = 0 for no coupling. With the same arguments as before, we can now write the stationary probability function as p(x, y) = p(x|y)p(y) with

$\begin{equation*}p(y)={N}_{y}\enspace \mathrm{exp}\left(-\left(-\frac{1}{2}{y}^{4}+{y}^{2}+4ky\right)/Q\right)\end{equation*}$

$\begin{equation*}p(x\vert y)={N}_{x}\enspace \mathrm{exp}\left(-\left(-\frac{1}{2}{x}^{4}+{x}^{2}+4kx-\beta {x}^{2}{y}^{2}\right)/Q\right).\end{equation*}$

In full analogy, the free energy is explicitly given by F = $(-\frac{1}{2}{x}^{4}+{x}^{2}+4kx-\beta {x}^{2}{y}^{2})/Q$ and provides an illustrative example of Friston's free energy principle, where the essential deterministic and stochastic features of the correlations between x and y are captured via the interaction term with β. For β = 0 the two nodes are statistically independent and display each a bimodal distribution illustrated in figure 9. Each individual fixed point can be locally approximated by a Gaussian, which can be analytically derived by a Taylor expansion of the full probability density function around the fixed point. As the coupling strength β increases toward 1, the stationary probability distribution changes its shape approaching the circular manifold and structuring its flow on the manifold. Distributions for various β values are plotted in figure 9. Here the structured flow comprises the four stable and four unstable fixed points.

Zoom In Zoom Out Reset image size