The future
Here we describe potential lines for future developments related to metastability in neuroscience. First, we discuss avenues for theoretical improvements, and we then focus on possible empirical work.
Beyond metastability — Turbulence
Turbulence, when a smooth flow breaks up into whorls and eddies, has been found to provide a fundamental principle governing optimal mixing and efficient transfer of energy/information over space and time [
116].
|
“Big whorls have little whorls Which feed on their velocity, And little whorls have lesser whorls And so on to viscosity.”
[117]
|
Recent research has shown how the brain is turbulent and how the scale-free nature of turbulence provides a dynamical regime where hierarchical information cascades allow the brain to function optimally despite its relative slowness [
118,
119,
120]. Turbulence has been demonstrated in fast local field potentials in local brain regions [
121] as well as in whole-brain dynamics measured with magnetoencephalography [
122].
A key insight relating metastability and turbulence comes from Kuramoto's pioneering research, which showed that models of coupled oscillators can be used to capture turbulence in many other systems [
32]. He showed that the Kuramoto order parameter can be extended to include information about space to provide a measure of spatiotemporal metastability, as the variability across spacetime of the
local Kuramoto order parameter R [
123]. Specifically,
Rn(t), is defined as the modulus of the local order parameter for a given brain node as a function of time:
Where
are the phases of the timeseries and
Cnp is the anatomical distance rule connectivity matrix
and
r(n,p) is the Euclidean distance between regions
n and
p, and
the decay.
In other words, Rn defines local levels of synchronisation at a given scale, , as function of space and time. This measure captures what could be called brain vortex space, Rn, over time, inspired by the rotational vortices found in fluid dynamics, but of course not identical.
As such, the level of amplitude turbulence
D can be defined as the standard deviation of the modulus of Kuramoto local order parameter and can be applied to the empirical data of any physical system. The amplitude turbulence [
124],
D, corresponds to the standard deviation across time and space of
R:
where the brackets < > denotes averages across space and time. This measure has been used to show turbulence in relatively slow neuroimaging fMRI data using a parcellation of 1000 regions (Cruzat et al., 2022; Deco, Kemp, et al., 2021; Deco, Sanz Perl, et al., 2021; Escrichs et al., 2022).
In other words, local metastability is measuring the local vorticity and is in essence an extension of the concept of global metastability (Kawamura et al., 2007), introduced in neuroscience to measure the variability across time of the global level of synchronisation of the whole system described here (Cabral et al., 2014; Kitzbichler, Smith, Christensen, & Bullmore, 2009; Kuramoto, 1984; Shanahan, 2010; Tognoli & Kelso, 2014; Wildie & Shanahan, 2012).
However, in order to show turbulence in fast neuronal dynamics of magnetoencephalography which has less good spatial resolution, a new edge metastability measure was introduced (Deco et al., 2023). This metric captures phase difference rather than mean phase and is defined as the standard deviation of the edge centred matrix
E. In this matrix, each column corresponds to a time point
t, where the column is defined as a vector combining all pairwise combinations of the differences of the phases’ state at time
t. With
N parcels this results in
N(N-1)/2 pairs. The difference of the phase state at a given time point
t for the different pairs is given by:
Where i,j are the parcels in a given parcellation. The edge metastability is thus the standard deviation across space and time of the edge-centric measure.
As such, the concepts of local and edge metastability in brain dynamics provide a way of generalising vorticity in fluid dynamics of the brain.
Beyond first order coupling in the Kuramoto model
The Kuramoto model [
32] consists of a population of
N coupled phase oscillators
i (
t) with natural frequencies
i , whose dynamics are governed by the following equation:
As discussed in the previous section, the variance of the instantaneous phase synchrony (i.e., the variance of the modulus of the complex-valued Kuramoto order parameter) is one of the most popular signatures of metastability [
33]. However, by construction, the corresponding order parameter is a mean-field value for in-phase synchrony — in other words, this measure of metastability only takes into account mean phase synchronisation, and other types of patterns (e.g., anti-phase) are neglected [
53,
125]. Therefore, the switching within the metastable regime considered by this signature only includes changes between disorder and partial order where in-phase cluster synchronisation and desynchronisation coexist, and ignores partial order where in-phase and antiphase clusters coexist with low desynchronised activity [
68,
82]. We will see in the next section that it is necessary to expand the Kuramoto model used in [
33] to allow for antiphase cluster synchronisation.
Kuramoto model with 2 Fourier modes
In the Hansel-Mato-Meunier model [
63], the canonical Kuramoto model was extended to include the 2nd Fourier mode in a model of weakly excitatory, synaptically coupled Hodgkin–Huxley neurons with phase lag. A similar modification was applied in a Hindmarsh-Rose model of neuronal bursting with-time delayed coupling [
64]. The general form of the extended model is:
A coupling function of the form:
was used by Hansel et al. (1993), where
is relative phase, α de-phases the relative position of the two modes which is similar to introducing a time delay, r modifies the contribution of the second order antiphase locking term, and the positive sign serves to destabilise the antiphase attractor. Three types of dynamics were found in the model: a fully synchronised state with one cluster, a totally incoherent state, and pairs of two-cluster saddle states linked by heteroclinic loops. Switching between in-phase and anti-phase synchrony in a 2-cluster state was observed in the transient dynamics before the system stabilised into 2-phase-locked clusters with a fixed phase-offset. The switching behaviour was dependent on the control parameter
which controlled the phase lag, or equivalently the time delay of the oscillators. Adding a small noise term led to a slow periodic switching between the two-cluster states, where the frequency of oscillations was proportional to the logarithm of the noise intensity (see Supplementary Video 3).
Building on the work of Hansel [
63], Kori & Kuramoto [
64] developed a second phase reduced model based on the Hindmarsh-Rose oscillator model of neuronal bursting [
65,
126] which included a time-delay as a phase shift [
64]. The coupling function of the phase-reduced model was determined numerically. In this model, when the oscillators clustered non-symmetrically, a saddle-node bifurcation transferred stability to the unstable dimension of the saddle-node. A heteroclinic loop then linked these two solutions (
Figure 6). They also found that slow switching could be induced by randomly distributing the time delay and/or breaking the uniformity in the coupling.
Generalised Haken-Kelso-Bunz model
Finally, we return to the roots of metastability in human movement coordination. The extended HKB model [
23,
127] was initially developed to model metastable coordination between two oscillatory components (for a review see [
128]). More recently, driven by experimental observations of in-phase and antiphase coordination and metastability in multiagent interaction [
49], a N-dimensional generalisation of the extended HKB model was developed [
50,
125], here referred to as the generalised HKB model. The equation of motion for the generalised HKB model has a coupling function that includes the first two Fourier modes:
where
is the phase of the i-th oscillator,
is the natural frequency of the i-th oscillator,
is relative phase, and a and b are the first and second coupling parameter, respectively. The model was initially developed as a phase-oscillator model in the form above [
50]. Later theoretical work confirmed that this phase-oscillator model can be derived from the N-dimensional van der Pol-Rayleigh oscillators model via phase reduction, thus consistent with the initial derivation of the dyadic HKB models [
129]. Note that the second order coupling is negative (c.f. [
63]), which stabilises the antiphase attractor. Adding the second order coupling was key to capture all experimental observations, which the canonical Kuramoto model could not (compare Supplementary videos 1 (Shanahan model) and 2 (generalised HKB model)). See the Technical Appendix for additional discussion of this model.
There is an assumption that the phenomenon of ghosts in dyadic interactions holds for higher dimensional interactions. The metastability in the Hansel et al., and Kori and Kuramoto models, was dynamically realised with heteroclinic cycles where the system slows down progressively as it approaches a saddle for longer periods between cycling. Interestingly, the intermittent converging behaviour of the General HKB model is not dissimilar from the Hansel et al. model, suggesting a similarity in the dynamics (see
Figure 7).
Avenues for future empirical work
Metastability has been proposed as a principle of brain function that reflects the simultaneous tendencies for global integration and functional segregation. It would be interesting to see how signatures of metastability correlate with these tendencies in clinical groups, pharmacological challenge studies, and studies in disorders of consciousness.
Although many questions are still open on the nature of metastability, an interesting dilemma relates to the saddle-type feature that facilitates metastable dynamics in the brain - is it a ghost left after a saddle-node bifurcation, a destabilised Milnor attractor, or a conventional saddle, or even all three, when viewed from different levels of observation? [
9]. Perhaps extending the Kuramoto model to include a second term to account for antiphase synchronisation provides a modelling framework where this dilemma could potentially be investigated.
Markers of metastability in empirical data tend to be time-averaged measures, and the underlying system has been shown to contain mono-, multi- and metastable regimes and transitions between them. Perhaps it would make sense to calculate these markers over windows of brain activity and compare the values with changes in similar windowed spatiotemporal patterns of activity, such as in time-varying functional connectivity. Investigating a time-varying marker of metastability at local, i.e., community level, and comparing across all communities, could identify if any particular community leads the transitions between mono-, multi- and metastability, if such regimes exist in the data.
The catchphrase terms of metastability and multistability are often used interchangeably and indiscriminately, as noted by [
4,
68] and [
31] respectively. Before describing some observation as metastable, it would be helpful to provide the rationale for this description. For example, in [
97], the brain waves found in a computational model were classified as metastable as they appeared and changed without any perturbation or added noise to the system, and transitions among metastable waves were identified through a interhemispheric cross-correlation function. Whilst a visual confirmation of metastability is possible with a limited number of oscillatory components as in [
49], this solution is not tractable for whole-brain models with a large number of parcellated brain regions. However, it would be interesting to investigate the dynamics of the generalised HKB model [
50,
125] to look for evidence of ghosts or heteroclinic cycles, and subsequently develop a large-scale brain model using a network of coupled neural masses as in [
97] to investigate if the dynamics are similar at larger scales (see [
130]).
There are also many open questions about the relationship between metastability and biophysical properties of the brain. With the release of NeuroMaps [
131], differences in local metastability could be mapped to myelination, metabolism, receptor and transmitter density, developmental expansion, or genomics. This could potentially provide insights into the biophysical underpinnings of metastability as it differs across psychiatric conditions, throughout healthy ageing, across disorders of consciousness and within development.
Specifically, we would like to establish the molecular, metabolic, physiological and cellular determinants of metastability, keeping in mind that one would expect these neurobiological factors to display not only reciprocal interaction but also bilateral relations with higher functions and environmental stimuli [
132]. For example, pronounced propagation delays have been posited to hinder information integration [
133]. Reduced integration has been found to strongly correlate with increased metastability [
53]. Propagation delays are dependent on myelination, which in turn is linked to metabolism, and in particular oxygen tension [
134]. Higher oxygen tension appears to favour the growth of neurons and oligodendrocytes, and a reduction in partial pressure results in thinner myelin sheaths [
135]. Reduced myelination in left temporal white matter was found in patients with recent onset psychosis [
136], and numerical density of oligodendrocytes was found to be significantly reduced in post-mortem studies, together with signs of compromised energy and redox metabolisms [
137]. However, myelin plasticity [
138] is modulated by astrocytes and microglia, and so one possible research question is to what extent glia contribute to metastability.
Another cellular question relates to the role of parvalbumin-expressing GABAergic interneurons in the generation of gamma-band oscillations which are altered in schizophrenia [
139,
140,
141,
142] and which drive neurovascular coupling [
143]. The E/I ratio reflects a homeostatic plasticity that keeps the relative levels of excitatory and inhibitory drive within a narrow range for healthy functioning [
144]. Alterations in the E/I ratio have been consistently found in schizophrenia [
145,
146,
147,
148,
149], and the E/I ratio is posited to be dependent on continual adaptation of the synaptic strength of PV GABA interneurons [
150,
151,
152]. The high metabolic demand required to maintain gamma oscillations puts strains on PV GABA interneurons making them susceptible to neuroinflammation and oxidative stress [
134,
151]. Immune gene expression analysis identified an immune fingerprint that could discriminate early psychosis from chronic schizophrenia and health controls [
153], and both reduced antioxidant status and a pro-inflammatory imbalance were found with first-episode psychosis [
154]. Brain metabolism in schizophrenia patients was modified with clozapine with concomitant improvements in symptoms. Metabolic activity measured with [18
F]Fluoro-deoxy-glucose PET was found to be reduced in the prefrontal, motor, and basal ganglia, and promoted in the visual cortex with improvements in disorganisation, negative and positive symptoms [
155]. Establishing whether any or all of the above cellular, molecular, and metabolic disease markers are determinants of metastability, could lead to improved understanding of the pathophysiology of schizophrenia, and potentially lead to new treatment targets.
Ultimately, the future of metastability as a theoretical construct rests upon demonstrating its causal roles in behaviour, cognition, and neuropsychiatric disorders, and the potential to leverage them in clinical settings. In the study of human movement coordination, metastability was often induced via causal manipulation of the control parameters, instrumentalized as, for example, pacing metronomes [
49,
156]. However, applications of the concept of metastability to common measures of neural activities have largely been observational, that is, without the causal manipulation of a control parameter guided by specific theoretical assumptions about the route to metastability (e.g., phase lag, delay, frequency differences). Developing methods to causally induce metastability in simultaneously measured brain activities in a theoretically predictable manner would further provide the opportunity to reveal the causal relationships between network specific metastability and cognitive functions and psychiatric symptoms. Notably, non-invasive brain stimulation, such as transcranial magnetic stimulation (TMS) and transcranial alternating current stimulation (tACS), has become an increasingly important tool for probing the causal roles of different brain oscillations in cognition and psychiatric disorders, with applications in clinical treatments (e.g., the use of alpha frequency stimulation to reduce depressive symptoms [
157,
158,
159,
160]). Currently, rhythmic brain stimulations are primarily designed to engage or entrain endogenous brain oscillations at their natural frequency rather than to induce metastability. However, existing techniques in non-invasive brain stimulation are fully equipped to control delay in coupling and frequency differences between the stimulator and endogenous oscillations, which provides promising grounds for inducing metastability. How to control the metastable interaction between existing endogenous oscillations using one or more external stimulators remains an open theoretical question.
In summary, the aim of future research in this area should strive to predict the behaviour of metastability with computational models, and then, supported with biological understanding, look for ways to return metastability to a healthy regime through appropriate pharmacological, psychotherapeutic, or non-invasive brain stimulation control strategies.
Glossary
Dynamical systems theory
A branch of mathematics that studies how the state of
systems evolve over time based on either an analytical (pencil and paper), or a
geometric (shapes), or a numeric (approximations using a computer) study of
deterministic evolution equations.
The state of a dynamical system
A vector representation that fully determines the value of
its variables (e.g., the position and velocity of a given particle). The
temporal evolution of the state of a dynamical system can be mathematically
described by a difference (i.e., discrete) or differential (i.e., continuous)
equation.
Phase space
The set of all possible states, and hence contains all the
allowed combinations of values of the system’s variables (also known as state
variables).
Trajectory of the system
A sequence of states within the phase space that satisfy
the dynamics of the system as defined by its differential equation.
A phase portrait
A geometric representation of the trajectories of a
dynamical system.
Potential landscape
The state
space of some dynamical systems can be illustrated by a potential landscape,
where valleys/dwells represent stable fixed points which correspond to
attractors, peaks represent unstable fixed points which correspond to
repellers, and balls represent the states of the system.
Manifold
A geometric model that can represent a broad range of
shapes. While at each point of a manifold it looks like a flat, n-dimensional
(euclidian) space, their overall structure can be highly non-trivial. The
dimensionality of the manifold, n, is the same for all its points.
Equilibria or fixed points
Points in the state space where the rate of change of the
system with respect to time is equal to zero, corresponding to states at which
the system remains unchanged unless perturbed.
Attractor
When many trajectories converge to a set of states, that
set is called an attractor.
Chaotic attractor
An attractor that holds dynamics that are highly sensitive
to their initial conditions.
Basin of attraction
All the points in phase space that flow onto the attractor.
Repeller
A set of states from which many trajectories migrate.
Saddle Point
A fixed point that is stable in one direction but unstable
in another, that is it behaves as an attractor for some trajectories and as a
repellor for others.
Heteroclinic orbit or cycle
A path in phase space that links saddles
Homoclinic cycle (or loop)
A path that links a saddle back onto itself.
Control parameters
A parameter that modifies the system of differential or
difference equations, hence deforming the corresponding flows through phase
space.
Critical point
The value of a control parameter at which a bifurcation
occurs.
Order parameter
A single variable that captures the collective or macro
behaviour of a system composed of microscopic elements.
Saddle-node bifurcation
When a stable and an unstable fixed point collide and
annihilate each other.
Bifurcation
A qualitative change in dynamics produced by varying a
control parameter in a dynamical system.
Bifurcation diagram
A diagram follows how the landscape changes with the
control parameter showing the location of attractors and repellers.
Phase transition
A concept originated from statistical physics. In general,
a system is said to undergo a phase transition when a small change in a control
parameter (e.g., temperature) causes a large collective change (e.g., the
system going from liquid state to solid).
Non-equilibrium phase transition
A phase transition that occurs in a physical system—like
the brain—that is far from equilibrium.
Synchronisation manifold
A smooth surface onto which the orbits of synchronised
chaotic attractors converge.
Metastability
A specific type of dynamics that may take place in a system
characterised by patterns that recur either in repeatable sequences
(pattern) or flexible alternation (no pattern).
Global cohesion
Mean functional connectivity computed using Pearson’s
correlation.
Spectral radius
The maximum of the absolute values of its eigenvalues for a
square matrix.
Relative coordination
Where occasional slippages between coordinating components
are balanced by an intrinsic attraction to certain preferred phase relations.
Critical fluctuations
As a control parameter approaches a critical point,
fluctuations increase, and this enhancement of fluctuations signals an upcoming
phase transition.
Bifurcation memory
Fixed points disappear after a saddle-node bifurcation, but
the “memory” of the fixed point remains attractive for the system and the
dynamics become very slow due to a bifurcation delay.
Ghost
Fixed points disappear after a saddle-node bifurcation, but
the “memory” of the fixed point remains attractive for the system and the
dynamics become very slow due to a bifurcation delay.
Chaos
A form of dynamical behaviour that can arise from time-invariant
nonlinear system. Chaos is characterized by sustained aperiodic (nonrepeating)
oscillations, leading to extreme sensitivity of future states to small changes
in present values of the system.
Chaotic itinerancy
The behaviour of complicated systems with weakly attracting
sets, where destabilised attractors allow the system to leave its basin of
attraction for another through heteroclinic orbits.
Milnor attractor
An attractor that has lost asymptotic stability, that is,
there exists a number of repelling sets.
Stable heteroclinic channel
A sequence of successive metastable (saddle) states in the
phase space
Chimera state
A state that represents the co-existence of coherent and
incoherent dynamics
Edge centric matrix
A matrix that contains the pairwise phase difference
between brain regions stored in upper triangle format.