Global dynamics of neural mass models

Abstract

Neural mass models are used to simulate cortical dynamics and to explain the electrical and magnetic fields measured using electro- and magnetoencephalography. Simulations evince a complex phase-space structure for these kinds of models; including stationary points and limit cycles and the possibility for bifurcations and transitions among different modes of activity. This complexity allows neural mass models to describe the itinerant features of brain dynamics. However, expressive, nonlinear neural mass models are often difficult to fit to empirical data without additional simplifying assumptions: e.g., that the system can be modelled as linear perturbations around a fixed point. In this study we offer a mathematical analysis of neural mass models, specifically the canonical microcircuit model, providing analytical solutions describing slow changes in the type of cortical activity, i.e. dynamical itinerancy. We derive a perturbation analysis up to second order of the phase flow, together with adiabatic approximations. This allows us to describe amplitude modulations in a relatively simple mathematical format providing analytic proof-of-principle for the existence of semi-stable states of cortical dynamics at the scale of a cortical column. This work allows for model inversion of neural mass models, not only around fixed points, but over regions of phase space that encompass transitions among semi or multi-stable states of oscillatory activity. Crucially, these theoretical results speak to model inversion in the context of multiple semi-stable brain states, such as the transition between interictal, pre-ictal and ictal activity in epilepsy.

Introduction

The surface of the human brain is covered by a cortical layer of grey matter. This cortex contains histologically distinguishable layers with laminar-specific types of neural cells that show patterned interlaminar connectivity [1]. The cortex is also composed of cortical columns, that show denser interlaminar synaptic connectivity—within a column—than horizontally among adjacent columns [2]. Ensemble neuronal activity within each cortical column generates extracellular currents, which can be measured using microelectrodes with relatively high spatial precision; or more coarsely as net average currents through local field potential recordings [3]. Please see glossary list, if required, for technical terms in the paper (S2 Text).

There are numerous approaches to modelling the dynamics of neuronal activity at the scale of cortical columns. One-dimensional models describing simplified single-neuron dynamics with threshold firing; namely, integrate-and-fire neurons [4] can model the impact of interneuron connectivity on larger scale dynamics, such as those observed within and between cortical columns. However, these equations often result in integral-differential equations or partial differential equations (please see S2 Text for definition), which are usually difficult to solve, at least in closed form, and especially when considering larger neuronal systems [5,6].

Using approximation methods from statistical physics, the collective activity (or ensemble activity) of cortical microcircuits comprising populations of distinct cell types have been modelled as neural mass point processes (please see S2 Text for definitions) [7–12]. There are broadly two sets of neural mass models: (1) Conduction-based models based on the Hodgkin-Huxley model, with specific modelling of the intrinsic dynamics of ion channels [13], which can allow for multi-compartment extensions [14,15]; and (2) convolution-based models where simple synaptic kernels are used to estimate the postsynaptic depolarization from the presynaptic input. These convolution-based models are computationally efficient whilst allowing a range of complex dynamics to be simulated. The canonical microcircuit model (CMC) is a typical example of these convolution-based neural mass models. It comprises 4 populations and has been used extensively to model EEG and MEG data [16]. In this setting, an idealized cortical column is modelled using a set of 2nd-order differential equations with nonlinear coupling, using a sigmoid function to map the effect of the neuronal potential of one population on the depolarization of another (please see S2 Text for definitions).

In this study, we present an analysis of the dynamics of convolution-based neural mass models, specifically, the canonical microcircuit model [16]. When considering the spatially coarse-grained dynamics—typically recorded in human neurophysiological recordings—these kinds of models afford a balance between accuracy and complexity. In other words, they are sufficiently expressive to capture neuronal dynamics, at the scale of a cortical column, while being sufficiently simple to preclude overfitting, when used to explain empirical timeseries.

Electrical activity of the cortex can be readily recorded using electro- or magnetoencephalography (EEG/MEG). Visual analysis of the ensuing timeseries reveals a mixture of apparently stochastic features, such as paroxysmal discharges, intermixed with recurrent rhythms. The spectrum of rhythms and paroxysmal discharges vary in frequency, location and progression and are seen in a variety of states of the normal and diseased brain [3]. Quantitative analysis of brain states suggest that these recurrent, recognizable patterns can be modelled as semi-stable states, i.e., states with local but not global stability [10,17–20]. Importantly, these dynamical features are sensitive markers for whole-brain (dys)function and can indicate pharmacological and pathophysiological changes of neuronal connectivity at the synaptic level [21–23]. Thus, developing biophysical models of cortical dynamics may help identify the synaptic and connectivity processes that shape paroxysmal and time-varying dynamic features of brain activity.

The dynamical structure of these systems—such as neural mass models—can be characterized using a phase space representation (please see S2 Text for definition). Simulation of the dynamics of neural mass models of the cortical column have disclosed complex structures in this phase space [24–26]. These include stationary points, limit cycles and chaotic attractors (please see S2 Text for definitions). Yet many popular current approaches of fitting neural mass models to empirical data, such as dynamic causal modelling, assume that the dynamics exist in a quasi-steady state, without transitions between distinct dynamical modes [27]. Thus, there is a need to develop the tools necessary to link biophysical models of cortical function to more complex, time-varying, and paroxysmal brain dynamics (e.g., seizure activity in epilepsy).

The theoretical and empirical applications of neural mass models—with separation of time scales—to analyze electrophysiological data, especially seizure activity, are diverse [28–33]. A seminal paper, by Jirsa et al., elegantly modelled epileptic seizure activity, in terms of bifurcations, using a neural mass model with two neural masses and a “permittivity” variable which introduced a separation of time scales. The slow permittivity variable reflects chemical changes, such as variations in electrolyte or metabolic homeostasis. This work has further been elaborated in computational pipelines that can be applied to specific patient data recorded in the evaluation for epilepsy surgery [34]. The Epileptor model differs from the canonical micro-circuit model, as the latter does not have explicit multiple time scales. However, we show in this paper that neural mass models do have an implicit time scale variation, which could explain the phenomena described by the Epileptor model. The method we use to study the implicit time scales of neural mass models are standard techniques in the study of dynamical systems [35,36] and was alluded to by the authors of the Epileptor model [37].

In what follows, we will estimate the semi-stable (i.e., multi-stable) states of the canonical microcircuit (CMC) model. We use an adiabatic approximation to integrate over fast changing variables—to derive the implicit slow time scale dynamics. To this slowly evolving state we add Brownian noise (to model unknown variations of the model, see S2 Text for definition) allowing for transitions between semi-stable states. More specifically, the dynamics of mean activity become gradient flows on a potential function (please see S2 Text for definitions). Crucially, this formulation can be used to generate key measures such as the frequency and duration of transient oscillatory modes, ie. features of itinerant activity. Using this formulation, we briefly consider variational methods that could be used to estimate synaptic connectivity from empirical measures of occupancy and mean exit times from distinct modes of activity. The theoretical results speak to a potential characterization of dynamical itinerancy seen in disorders of the brain—including epilepsy—in computationally tractable ways that can be incorporated in established (variational) model inversion schemes, such as dynamic causal modelling.

Methods

This section describes the formal results upon which the proposed analysis rests. This analysis can be summarised, in narrative form, along the following lines:

* First, identify a sufficiently expressive neural mass model, whose equations of motion are parameterised in terms of synaptic connectivity.

* With a suitable transformation of variables, separate the fast (oscillatory) dynamics from slow (amplitude) fluctuations, in the spirit of an adiabatic approximation.

* Formulate the dynamics of the slow variables in a relatively simple mathematical form. More exactly, define the dynamics as a gradient flow on a potential function. This potential function can now be taken as the negative logarithm of a steady-state density; namely, the solution of the Fokker Planck equation describing density dynamics under random fluctuations (please see S2 Text for definitions).

* From the steady-state density, evaluate the probability of occupying various fixed points—in the state space of the slow variables—to generate the statistics of dwell times and transitions.

* Use the above steps as a generative model—that generates the statistics of dwell times from synaptic connectivity—to infer connectivity from empirically observed statistics of (slow) transitions among (fast) oscillatory modes of activity.

In what follows, we describe the mathematical basis of the steps above and present a brief, illustrative application, using simulated data to recover the connectivity giving rise to transient bursts of oscillatory activity.

CMC-model dynamics

The electrical activity of the cortical column can be modelled using the CMC-model as described in [16,27]. To keep the derivations simple, and to ensure readability, we will rename some parameters here, without changing the dynamics of the system. The model considers four distinct neuronal populations modelled using a set of four 2nd order differential equations. These ordinary differential equations can be written equivalently as eight first order equations.

Where the potential of the i-th neural population is given by qi and the rate of change of this is given by pi (i.e. a current term). We have introduced 2 variables μ and ε to control the effect of perturbation of the non-connected model where the connectivity matrix is given by gij. A sigmoid function (S) is used to parametrize the interaction between the potential and current terms of the neuronal populations. The P-function parametrizes spike rate variability.

The original variables (qi, pi) can be written as a complex variable (∈ℂ) which will keep the derivations relatively concise and easy to follow.

The closed loops of the unperturbed motion will be given by,

The equation of unperturbed motion will be given by,

The solution will be given by,

The motion of the trajectories can then be given by the modulus (R) and argument (θi) of the complex number (zi), Fig 1.

[Figure omitted. See PDF.]

Fig 1. Schematic figure showing the relation between the variables (q,p) and the amplitude and angle variables.

The blue curve shows the trajectory of an unperturbed trajectory (no damping and no coupling to other neuronal populations). The red curve depicts the trajectory of a perturbed path. The amplitude (R) and angle variable (θ) for a point of the new trajectory (red) is shown.

https://doi.org/10.1371/journal.pcbi.1010915.g001

We introduce a new phase variable φi describing the deviation of the phase flow around its value when μ and ε = 0, i.e. no perturbation.

The resulting equations of motion can then be given by,(1)

Adiabatic approximation

The above equations of motion can be simplified under the assumption that the time scale for phase dynamics (which is determined by ωi) is small in relation to the time scale for amplitude dynamics. This assumption has empirical validity when modelling cortical activity, as the amplitude modulation of cortical rhythms or spike amplitude evolves over a slower time scale (at least 10–100 times slower) than that of the phase dynamics [28,38].

In the above formulation this adiabatic assumption can be satisfied by limiting the range of ε and μ. The resulting adiabatic approximation is a simplified model of the amplitude variables (i.e., a 4-dimensional model), which permits further analysis of the dynamics in terms of closed form solutions. The usefulness of the ensuing adiabatic approximation can be assessed by the ability of the model to explain empirical data (see final section).

An adiabatic approximation can be specified with a time duration (T), where each neuronal population will complete—to 0th order (see Perturbation analysis)—a series of cycles. The angle variables move on a 4-dimensional torus (please see S2 Text for definition), and T is chosen such that the motion over the torus completes (or approximately completes) one cycle and the angle variables can be eliminated from the dynamical flow [35].

Results

Further analysis of the dynamics

A perturbation analysis of Eq 1 can now be pursued using the following expansion for φi and Ri (in modulus-argument form, please see S2 Text for definition).

Where the two terms will be given by,

From the flow of states, we get the following,

Which can be written as,

The resulting equation can now be integrated over the closed (or almost closed) path over the torus, such that the angle variables are eliminated. To perform this integration we expand the perturbative functions S and P as polynomial series. We are now in a position to derive the change in Ri to 1st order, by combining the results derived in the Appendix (S.1.1–3):

The dynamics of Ri is determined by 2 coupling terms: the potential-to-current coupling (S) and the current-to-current coupling (P). The nonlinearity of the coupling functions (S and P) result in a series expansion in multiples of Ri. The complexity of the dynamics for the amplitude parameters is determined by the connectivity matrix and the coupling functions.

These parameters are defined in the cortical column by reference to the synaptic connections between the neuronal populations, where the gain of each connection is given by gij and their physiological effect on the target neuron parametrized in S and P. Note that the above derivation is valid for any set equal to or larger than two neuronal populations.

Dynamics of amplitude modulations

The electrical or current recordings of cortical activity is often modelled as a linear combination of the activity of the potential of each neuronal population [39]. This can be used to derive a 1-dimensional equation of motion for the amplitude dynamics of the recorded data. As can be seen above the amplitude dynamics will be governed by the self-connection terms, gii. The self-interactions within the cortical column will result in a set of non-interacting limit cycles (taking up-to 1st order contribution from the perturbation).

Amplitude envelope dynamics

The measured potential (e.g., EEG) is given by (assuming a linear lead field).

(2)

The amplitude power can be estimated easily, provided the ωi are distinct.

We can now normalize qi s.t

R will then be the radius of a sphere if all ai are constant in the vector space given by Ri. The rate of change of the amplitude function R is given as below.

As we are interested in the average rate of change of R we will take a spatial average over ellipsoids with constant R2.

R will be constant on the surface,

Hyperspheroidal coordinates (please see S2 Text for definition) are used for the integration (see S1 Text, S1.4). Let Rs be the average of R over the 3-sphere,(3)

As was discussed in the appendix, cross connectivity terms (or off-diagonal terms) will affect the dynamics of Rs when higher order terms are included in the perturbation expansion. The functional form of this equation means that the amplitude dynamics can be expressed as a gradient flow on a potential, U that depends upon the synaptic connections gij and the coupling functions:

In other words, different connectivity matrixes and functions correspond to different potentials (U): see Fig 2. This spectrum then defines different possibilities of dynamics supported by the CMC model. Panel A shows minima at 0 and 2. Panel B has three stable points, at 0, 1 and 2 (see S1 Text, S1.1–3). This suggests that the CMC model—with certain synaptic connections—features continuous high amplitude oscillations and a point of stability at the origin. The oscillations would be a mixture of sinusoidal activity, with frequencies of the four neural populations, where this mixture is determined by the lead field (Eq 2).

[Figure omitted. See PDF.]

Fig 2. Potential distributions for the dynamics of the total amplitude (R) for different connectivity between the neuronal subpopulations.

Here different connectivity matrices and couplings are selected to illustrate the spectrum of potential functions. Minima of U in these functions are stable point attractors (A) The first example has two stable points at the origin and at R = 2. (B) The second figure has three stable states at the origin and at R = 1 and 2.

https://doi.org/10.1371/journal.pcbi.1010915.g002

The potential distributions U in Fig 2 are generated by one self-connectivity term in the superficial pyramidal cell layer. Up until this point, we have been considering deterministic dynamics in the absence of noise. In what follows, we show that the same dynamical structures persist in the presence of random fluctuations, via an analysis of the density dynamics.

A schematic of the steps taken to derive the dynamical equations of the amplitude modulation of the EEG is shown in Fig 3.

[Figure omitted. See PDF.]

Fig 3. A schematic figure of the steps taken to derive the amplitude dynamics.

The starting point is the 4 population CMC-model. The activity of the model generates activity for each of the four populations as shown after the arrow “4-population dynamics”. The EEG lead field will then take a linear mixture of the activity of the 4 populations resulting in a “EEG” curve. The amplitude modulation of the activity is drawn by a red curve. The adiabatic and perturbation analysis discussed in this paper will define a gradient field as drawn after the arrow “perturbation expansion”. The gradient field will depend on the connectivity matrix, G. A possible trajectory of the amplitude variable is shown as a red arrow moving from one stable state to another. The corresponding change in amplitude over time is shown after the arrow “gradient field dynamics” in red above the “EEG curve”. The main contribution of this work is to find a direct relation between the intrinsic connectivity and the gradient field which in turn will define statistical measures such as the probability of being in a specific state or exit times from states.

https://doi.org/10.1371/journal.pcbi.1010915.g003

Stationary states and exit times

Fast neuronal fluctuations can be easily modelled by adding innovations to the deterministic equations of motion [40,41]. For the potential flow of the amplitude dynamics derived above we can add a white (Wiener) noise process, resulting in a stochastic differential equation (SDE) (see S2 Text, for definitions).

The trajectories are then described by the following SDE,

In the above equation, I and Bt are random variables while U and σ are real valued functions. The associated Fokker Planck equation, describing the density dynamics over the amplitude, is given by:

The stationary state of this equation is given by the following, where Z is the partition function (please see S2 Text for definition) or equivalently, a normalizing constant for p).

With a connectivity matrix that features several stable states, we can now estimate the probability of finding the dynamical system in that state by simply integrating the stationary probability density function (i.e., solution to the Fokker Planck equation above) over the interval defining the state: see Fig 4. This provides an analytic estimate of noise induced multistability (c.f., topological supersymmetry breaking), in terms of the probability of finding the amplitude dynamics near one of the fixed points or stable (oscillatory) states.

[Figure omitted. See PDF.]

Fig 4. Regions attracted to different minima of the potential (shaded differently).

The probability of finding the system in one of these stable states is given by integrating the probability density function over the attracting interval (state 1 = light grey), 0 to the first maxima (Rmax) of the potential curve or Rmax to ∞ (state 2 = dark grey). The red arrow depicts how the trajectories pass between state 1 and 2 due to the underlying noise.

https://doi.org/10.1371/journal.pcbi.1010915.g004

For the first state the probability is given by,

The mean duration of time (t) the system is in a given state can be approximated using Kramer’s escape rate,

Finally, because we have closed form expressions for both the state probability and mean exit times, we are able to estimate the connectivity matrices from empirical measures of occupancy and mean exit times from discernible oscillatory states. In other words, the closed form expressions above provide the basis for a generative model that can be inverted given empirical measures of transitions among oscillatory modes evident in any empirical data. See Fig 5 for inversion of simulated data of exit times, for different noise amplitudes, using the variational procedures used in dynamic causal modelling. A simulated EEG trace is shown in Fig 6 illustrating how the system is moving between two semi-stable states resulting in a bursting appearance of the EEG.

[Figure omitted. See PDF.]

Fig 5. Inversion of the generative model based on simulated data.

The left panel shows the mean exit time (on the vertical axis) from state 2 (Fig 4) for different noise amplitudes (plotted along the horizontal axis). The dotted black curve indicates the synthetic data that was inverted. The blue curve shows the expected data under the prior values of the connectivity matrix. The red curve indicates the expected exit times with the posterior values of the connectivity matrix. The right panel shows the change in the connectivity matrix between the prior and posterior estimates of the connectivity matrix, suggesting that the connectivity parameters used to generate the data can be recovered from mean exit time data features.

https://doi.org/10.1371/journal.pcbi.1010915.g005

[Figure omitted. See PDF.]

Fig 6. Simulation of EEG trace with bi-stable activity.

The trace shows a simulated EEG trace with a connectivity matrix giving a bistable system as described in Fig 4. The system is bursting with high amplitude oscillations (state 2) with shorter suppressions of attenuation (state 1) in between.

https://doi.org/10.1371/journal.pcbi.1010915.g006

In this proof of principle example, synthetic data were generated from using prior values of the connectivity in the CMC model provided by the SPM software. Inversion of the synthetic data was performed using standard (variational Laplace) inversion routines available in the SPM software package. The code used in this paper can be obtained from https://github.com/gercoo/exittimes2synapticconnectivity. This example shows that the connectivity parameters used to generate transitions among oscillatory states or modes can be recovered from mean exit times alone. The implication here is that it is possible to both explain and leverage multi-stability in terms of underlying synaptic connections and the amplitude of fast neuronal fluctuations. This work could be applied when investigating spiking in epileptic and non-epileptic tissue as the inter-spike interval could be modelled as the exit time from a large amplitude transient. The inversion method described could then map spike frequency to synaptic connectivity, which could be used to uncover the underlying physiological differences between regions of the brain with different epileptogenicity.

Discussion

This paper has presented a formulation of multi-stable neuronal dynamics based on neural mass models—and an accompanying analysis framework—that enables inference about synaptic connectivity parameters that engender to observable transitions among oscillatory states. We have derived closed form expressions for the statistics of the semi-stable states, e.g. the state probability and the mean duration of a state.

Our work builds on both theoretical insights of prior neural mass modelling and empirical observations of transient brain states. Prior simulations of neural mass models have already shown the simultaneous presence of fixed point and limit cycle dynamics in the same cortical model [25,42]. Furthermore, empirically there are several instances—in both the healthy and diseased brain—where the mode of cortical and subcortical dynamics shifts repeatedly among states; including physiological rhythms like sleep spindles, mu-rhythms, and pathological states, such as epileptic seizure activity [3,43–47]. Indeed, the implicit itinerancy may be the hallmark of all self-organising dynamical systems; ranging from the stochastic chaos that characterises nonequilibrium steady-state [48–55] through to heteroclinic cycles [56–58].

The neural mass models we have considered here show complex dynamics that can be divided into explicit fast and implicit slow variables. This separation in time scales allows for an adiabatic approximation, which eliminates the fast modes producing a new set of dynamical equations of slow variables. This new set of slow variables now describes slowly evolving amplitude dynamics, rather than the fast oscillatory activity itself. This reduces the dimensionality of the system making it more tractable for simulation and inference. In this reduced system, several types of fixed points—representing stationary states or oscillatory limit cycles—can be characterised using measurable data features such as dwell and mean exit times.

In this work, we have presented closed form equations that link patterns of cortical microcircuit synaptic connectivity to statistical data features. This generative model allows for Bayesian model inversion and comparison, allowing empirical comparison of different mechanistic explanations of intermittent cortical dynamics. This might be particularly interesting when modelling relatively fast changing dynamics like on and off rhythmic spiking or seizure onset. As seizures progress there are changes in several variables including electrolytic or metabolic changes which would have a large but secondary effect on cortical dynamics and might change some key connectivity parameters, which would alter both fast and slow dynamics in the system described here [59]. These slower changes can also be modelled effectively using Adiabatic DCM where the effect of slow changes in synaptic parameters is included in the model inversion [60].

The amplitude dynamics of a single-source CMC model (i.e., cortical column) has been shown to be sufficiently approximated by a second order perturbation expansion. In contrast, extrinsic connections between cortical columns may require higher order perturbation terms and thus be of weaker strength compared to the intrinsic dynamics (at least within the parameter range where the perturbation expansion is valid). This is due to the number of neurons required to create an interconnected loop between two cortical columns [2,61]. The intrinsic (i.e., within-source or column) dynamics is thus instantiated in cortical columns through the intrinsic connectivity matrix, while the extrinsic connectivity results in a weak but richer set of dynamics, as the extrinsic (i.e., between-source or column) connectivity can be diverse. This is in line with the architectural setup of cortical processing, where cortical columns represent the smallest units that process data as conditionally independent units and interact with other columns, using a weaker set of interactions [62]. We will further develop the perturbation theory presented here to include higher order terms to see whether the interactions between cortical columns can be characterised.

Our perturbation expansion is sufficient to explain semi-stable states of the cortical modes of activity (i.e. stationary states and limit cycles of the full dynamics) within cortical columns in the absence of changes of connectivity. Using a perturbation expansion, we have managed to derive the characteristics of semi-stable states of the cortical columns from the connectivity matrix describing connectivity within the cortical column. This is an important step, as we can now map statistical characterisations of complex cortical dynamics onto synaptic connectivity. Understanding larger cortical networks from a modelling perspective can be further implemented using a discrete simplification of the above model. It is possible—through a projection mapping of the dynamics of the cortical columns onto the semi-stable state it occupies at any given time—to simplify things considerably. This should preserve the important link between connectivity and the semi-stable states that are generated. Connecting different cortical columns could then result in a lattice model (similar to an Ising or Pott’s model) of the cortex, where it would be relatively straightforward to link the intrinsic and extrinsic connectivity to measurable dynamics [63,64]. Such an expansion of our work would prove particularly relevant to the characterisation of epileptogenic networks and their role in seizure initiation and spread. Furthermore, the effect of medication could be addressed in terms of the changes in the intrinsic connectivity and in ceiling changes in network dynamics. Finally, the effect of epilepsy surgery or disconnection surgery could be modelled by changes in extrinsic connectivity; c.f. [65–67].

It is important to compare this work with other published methods using neural mass models to investigate itinerant brain dynamics. The Epileptor model is frequently used to characterize epileptic seizure onset, propagation and offset [28,34,66]. As described previously, the timescales are explicit in the functional form of the equations of the Epileptor model; while they are only implicit in the physiologically grounded (canonical microcircuit) model described in this study. The explicit parameterization in terms of distinct rate (i.e., time) constants simplifies the mathematical modelling of seizures and provides a compelling taxonomy of seizures in terms of bifurcations for the field of computational epilepsy. The CMC-model, in contrast, requires extensive mathematical derivations to uncover the underlying slow processes; however, the biophysics of slow modulation can be shown to depend on the synaptic connectivity within a cortical column. This may be important in terms of understanding the mechanisms (e.g., synaptopathies) that underwrite epilepsy.

We have thus shown the possibility of a cortical column evincing slow dynamics in the absence of slowly fluctuating variables. This dynamical mechanism fits well with the principle of self-organization, where a system shows patterns in time or spatial scale larger than the local dynamics governing the system. The analysis provided in this paper further allows for the statistics of transient electrophysiological activity (e.g., spike frequency and seizure duration) to be mapped onto the synaptic connectivity, in contrast—or as a complement—to the permittivity variable in the Epileptor model.

Finally, there is a key difference between accounts exemplified by the Epileptor model, and the model described in this paper. The dynamics of slow variables (e.g. amplitude modulation) derived in this paper are dependent on fast variables (e.g. high frequency activity within neural populations) only via the adiabatic approximation. In other words, the fast activity pushes the system configuration slowly giving rise to the implicit slow dynamics. We describe a causal flow from fast to slow timescales, in contrast to a feed-back system of the kind found in the Epileptor model, where the dynamics of the slow variables directly affect faster variables. The dependence of slow variables on fast variables in this paper allows us to define an “arrow of causality” that is similar to that found in physics, where macroscopic variables can be estimated from the underlying microscopic variables but not vice versa. [65,68–71]. However, we note that the feed-back from slow to fast timescales may be required in certain cases (e.g., periodic seizures or seizures of longer duration), perhaps using a similar set-up as the Epileptor model [28]. This speaks to some interesting issues about the nature of circular causality, and bottom-up versus top-down causation [66, 72]. These issues could, in principle, be addressed using Bayesian model comparison of appropriate dynamic (causal) models.

The approach illustrated here may have relevance for our understanding of paroxysmal disorders of cortical function such as epilepsy and epileptic seizures. The epilepsies are characterised by recurrent, sudden onset of rhythmic epileptic activity, on a background of spontaneous ‘normal’ interictal ongoing activity. In certain cases, cortical activity in epilepsy can alternate rapidly between more than two states, 1) rhythmic spiking, 2) seizure activity and 3) spontaneous background activity [17,73]. The transition between these states has been described in terms of phase transitions and bifurcations, where the synaptic parameters that shape microcircuit dynamics exceed a threshold, resulting in a sudden change in phase space dynamics [39]. Our model, similar to previous modelling endeavours, offers an alternative perspective in that it does not require a critical transition (as in transcritical bifurcations) to explain the change in cortical activity [74]. Rather, the mechanism resulting in intermittently observable dynamic phenomena is predicated on noise induced itinerancy, of the kind found in stochastic chaos and topological supersymmetry breaking. In other words, the stochastic nature of the model will move the dynamics of the system between semi-stable states. This form of transition between states has been described in terms of transcritical activity [39]. What is new in this study is that we can now map the observed changes to synaptic connectivity within the cortex.

Observations of intermittent dynamic features could—in principle—give support for either bifurcations or transitions between semi-stable states. Note that we do not provide any evidence for the absence of bifurcations in these intermittent systems, as these still exist within the framework presented (i.e. changing the connectivity matrix can still cause transcritical bifurcations): rather that these are not always required to explain observable state transitions. To distinguish between the two proposals; bifurcations or semi-stable state transitions, the statistics of transitions among states may need to be further explored. It is through tools such as the invertible, reduced models presented here that future work will be able to specifically address these questions empirically; i.e., through the (Bayesian) model comparison of models with and without changes in connectivity.

In summary, we have introduced the theoretical backdrop for a model of cortical multi-stability and itinerancy that could be used to infer the underlying synaptic connectivity from statistical descriptions of transient dynamic features. We hope to elaborate the theoretical work presented above for the characterisation of global dynamics based on measurable electromagnetic responses in health and disease.

Supporting information

S1 Text. Appendix with in-depth mathematical analysis.

https://doi.org/10.1371/journal.pcbi.1010915.s001

(DOCX)

S2 Text. Glossary list with technical terms used in the paper.

https://doi.org/10.1371/journal.pcbi.1010915.s002

(DOCX)

S1 Fig. 2D-ellipoids representing the surfaces with constant EEG amplitude.

The spatial average of the amplitude flow will be done over the positive sector (Ri>0) of each of these surfaces.

https://doi.org/10.1371/journal.pcbi.1010915.s003

(TIF)

S2 Fig. 2D-spherical and 2D-ellipsoidal area unit.

The spherical (yellow) and the ellipsoidal (blue) area unit are separated for better visualization. Using the spherical area unit instead of the elliptic will simplify the integrals considerably allowing for analytical solutions.

https://doi.org/10.1371/journal.pcbi.1010915.s004

(TIF)

Citation: Cooray GK, Rosch RE, Friston KJ (2023) Global dynamics of neural mass models. PLoS Comput Biol 19(2): e1010915. https://doi.org/10.1371/journal.pcbi.1010915

References

1. Jiang X, Shen S, Cadwell CR, Berens P, Sinz F, Ecker AS, et al. Principles of connectivity among morphologically defined cell types in adult neocortex. Science. 2015 Nov 27;350(6264):aac9462. pmid:26612957

2. Douglas RJ, Martin KA, Whitteridge D. A canonical microcircuit for neocortex. Neural computation. 1989 Dec;1(4):480–8.

3. Niedermeyer E, da Silva FL, editors. Electroencephalography: basic principles, clinical applications, and related fields. Lippincott Williams & Wilkins; 2005.

4. Abbott LF, Kepler TB. Model neurons: from hodgkin-huxley to hopfield. In Statistical Mechanics of Neural Networks: Proceedings of the Xlth Sitges Conference Sitges, Barcelona, Spain, 3–7 June 1990 2005 Jul 10 (pp. 5–18). Berlin, Heidelberg: Springer Berlin Heidelberg.

5. Coombes S., Waves bumps, and patterns in neural field theories. Biological cybernetics. 2005 Aug;93(2):91–108. pmid:16059785

6. Brunel N, Hakim V. Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural computation. 1999 Oct 1;11(7):1621–71. pmid:10490941

7. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological cybernetics. 1995 Sep;73(4):357–66. pmid:7578475

8. David O, Friston KJ. A neural mass model for MEG/EEG: coupling and neuronal dynamics. NeuroImage. 2003 Nov 1;20(3):1743–55. pmid:14642484

9. Wilson HR, Cowan JD. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik. 1973 Sep;13(2):55–80. pmid:4767470

10. Lopes da Silva FH, Hoeks A, Smits H, and Zetterberg LH. Model of brain rhythmic activity: The alpha-rhythm of the thalamus. Kybernetik, 15:27–37, 1974.

11. Wendling F, Benquet P, Bartolomei F, Jirsa V. Computational models of epileptiform activity. Journal of neuroscience methods. 2016 Feb 15;260:233–51. pmid:25843066

12. Dafilis MP, Frascoli F, Cadusch PJ, Liley DT. Chaos and generalised multistability in a mesoscopic model of the electroencephalogram. Physica D: Nonlinear Phenomena. 2009 Jun 15;238(13):1056–60.

13. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of physiology. 1952 Aug 28;117(4):500. pmid:12991237

14. Rall W. Branching dendritic trees and motoneuron membrane resistivity. Experimental neurology. 1959 Nov 1;1(5):491–527. pmid:14435979

15. Traub RD. Motorneurons of different geometry and the size principle. Biological Cybernetics. 1977 Sep;25(3):163–76. pmid:836920

16. Bastos AM, Usrey WM, Adams RA, Mangun GR, Fries P, Friston KJ. Canonical microcircuits for predictive coding. Neuron. 2012 Nov 21;76(4):695–711. pmid:23177956

17. Cooray CN, Carvalho A, Cooray GK. Noise induced quiescence of epileptic spike generation in patients with epilepsy. Journal of Computational Neuroscience. 2021 Feb;49(1):57–67. pmid:33420615

18. Wendling F, Bellanger JJ, Bartolomei F, Chauvel P. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biological cybernetics. 2000 Sep;83(4):367–78. pmid:11039701

19. Freeman WJ. Models of the dynamics of neural populations. Electroencephalography and clinical neurophysiology. Supplement. 1978 Jan 1(34):9–18. pmid:285858

20. Liu F, Wang J, Liu C, Li H, Deng B, Fietkiewicz C, et al. A neural mass model of basal ganglia nuclei simulates pathological beta rhythm in Parkinson’s disease. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2016 Dec 15;26(12):123113. pmid:28039987

21. Boly M, Moran R, Murphy M, Boveroux P, Bruno MA, Noirhomme Q, et al. Connectivity changes underlying spectral EEG changes during propofol-induced loss of consciousness. Journal of Neuroscience. 2012 May 16;32(20):7082–90. pmid:22593076

22. Smailovic U, Koenig T, Kåreholt I, Andersson T, Kramberger MG, Winblad B, et al. Quantitative EEG power and synchronization correlate with Alzheimer’s disease CSF biomarkers. Neurobiology of aging. 2018 Mar 1;63:88–95. pmid:29245058

23. Willoughby JO, Fitzgibbon SP, Pope KJ, Mackenzie L, Medvedev AV, Clark CR, et al. Persistent abnormality detected in the non-ictal electroencephalogram in primary generalised epilepsy. Journal of Neurology, Neurosurgery & Psychiatry. 2003 Jan 1;74(1):51–5. pmid:12486266

24. Ermentrout GB, Kopell N. Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM journal on applied mathematics. 1986 Apr;46(2):233–53.

25. Grimbert F, Faugeras O. Bifurcation analysis of Jansen’s neural mass model. Neural computation. 2006 Dec;18(12):3052–68. pmid:17052158

26. Touboul J, Wendling F, Chauvel P, Faugeras O. Neural mass activity, bifurcations, and epilepsy. Neural computation. 2011 Dec 1;23(12):3232–86. pmid:21919787

27. Moran RJ, Pinotsis DA, Friston KJ. Neural masses and fields in dynamic causal modeling. Frontiers in computational neuroscience. 2013 May 28;7:57. pmid:23755005

28. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014 Aug 1;137(8):2210–30.

29. Ermentrout B, Terman DH. Mathematical foundations of neuroscience. New York: springer; 2010 Jul 1.

30. Hindmarsh JL, Rose RM. A model of neuronal bursting using three coupled first order differential equations. Proceedings of the Royal society of London. Series B. Biological sciences. 1984 Mar 22;221(1222):87–102. pmid:6144106

31. Cressman JR, Ullah G, Ziburkus J, Schiff SJ, Barreto E. The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics. Journal of computational neuroscience. 2009 Apr;26(2):159–70.

32. Ullah G, Cressman JR Jr, Barreto E, Schiff SJ. The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: II. Network and glial dynamics. Journal of computational neuroscience. 2009 Apr;26(2):171–83.

33. Köksal Ersöz E, Modolo J, Bartolomei F, Wendling F. Neural mass modeling of slow-fast dynamics of seizure initiation and abortion. PLoS computational biology. 2020 Nov 9;16(11):e1008430. pmid:33166277

34. Jirsa VK, Proix T, Perdikis D, Woodman MM, Wang H, Gonzalez-Martinez J, et al. The virtual epileptic patient: individualized whole-brain models of epilepsy spread. Neuroimage. 2017 Jan 15;145:377–88. pmid:27477535

35. Arnold VI. Mathematical methods of classical mechanics. Springer Science & Business Media; 2013 Apr 9.

36. Shilnikov A, Kolomiets M. Methods of the qualitative theory for the Hindmarsh–Rose model: A case study–a tutorial. International Journal of Bifurcation and chaos. 2008 Aug;18(08):2141–68.

37. El Houssaini K, Bernard C, Jirsa VK. The epileptor model: a systematic mathematical analysis linked to the dynamics of seizures, refractory status epilepticus, and depolarization block. Eneuro. 2020 Mar;7(2).

38. Novak P, Lepicovska V, Dostalek C. Periodic amplitude modulation of EEG. Neuroscience letters. 1992 Mar 2;136(2):213–5. pmid:1641193

39. Mosher JC, Leahy RM, Lewis PS. EEG and MEG: forward solutions for inverse methods. IEEE Transactions on biomedical engineering. 1999 Mar;46(3):245–59. pmid:10097460

40. Friston KJ, Li B, Daunizeau J, Stephan KE. Network discovery with DCM. Neuroimage. 2011 Jun 1;56(3):1202–21. pmid:21182971

41. Oksendal B. Stochastic differential equations: an introduction with applications. Springer Science & Business Media; 2013 Mar 9.

42. Spiegler A, Kiebel SJ, Atay FM, Knösche TR. Bifurcation analysis of neural mass models: Impact of extrinsic inputs and dendritic time constants. NeuroImage. 2010 Sep 1;52(3):1041–58. pmid:20045068

43. Alvarado-Rojas C, Valderrama M, Fouad-Ahmed A, Feldwisch-Drentrup H, Ihle M, Teixeira CA, et al. Slow modulations of high-frequency activity (40–140 Hz) discriminate preictal changes in human focal epilepsy. Scientific reports. 2014 Apr 1;4(1):1–9. pmid:24686330

44. Brown P. Bad oscillations in Parkinson’s disease. Parkinson’s Disease and Related Disorders. 2006:27–30. pmid:17017505

45. Canolty RT, Edwards E, Dalal SS, Soltani M, Nagarajan SS, Kirsch HE, et al.High gamma power is phase-locked to theta oscillations in human neocortex. science. 2006 Sep 15;313(5793):1626–8. pmid:16973878

46. Cooray GK, Sengupta B, Douglas P, Englund M, Wickstrom R, Friston K. Characterising seizures in anti-NMDA-receptor encephalitis with dynamic causal modelling. Neuroimage. 2015 Sep 1;118:508–19. pmid:26032883

47. Sterio D, Berić A, Dogali M, Fazzini E, Alfaro G, Devinsky O. Neurophysiological properties of pallidal neurons in Parkinson’s disease. Annals of Neurology: Official Journal of the American Neurological Association and the Child Neurology Society. 1994 May;35(5):586–91. pmid:8179304

48. Ao P. Potential in stochastic differential equations: novel construction. Journal of physics A: mathematical and general. 2004 Jan 7;37(3):L25.

49. Agarwal S, Wettlaufer JS. Maximal stochastic transport in the Lorenz equations. Physics Letters A. 2016 Jan 8;380(1–2):142–6.

50. Friston K, Heins C, Ueltzhöffer K, Da Costa L, Parr T. Stochastic chaos and markov blankets. Entropy. 2021 Sep 17;23(9):1220. pmid:34573845

51. Jiang Da-Quan, and Jiang Donghua. Mathematical theory of nonequilibrium steady states: on the frontier of probability and dynamical systems. Springer Science & Business Media, 2004.

52. Kim EJ. Investigating information geometry in classical and quantum systems through information length. Entropy. 2018 Aug 3;20(8):574. pmid:33265663

53. Lasota A, Mackey MC. Chaos, fractals, and noise: stochastic aspects of dynamics. Springer Science & Business Media; 1998 Apr 1.

54. Nieolis G, Prigogine I. Self-Organization in Nonequilibrium Systems. New York: Wiley. 1977.

55. Yan H, Zhao L, Hu L, Wang X, Wang E, Wang J. Nonequilibrium landscape theory of neural networks. Proceedings of the National Academy of Sciences. 2013 Nov 5;110(45):E4185–94. pmid:24145451

56. Friston KJ. The labile brain. II. Transients, complexity and selection. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences. 2000 Feb 29;355(1394):237–52. pmid:10724458

57. Rabinovich M, Huerta R, Laurent G. Transient dynamics for neural processing. Science. 2008 Jul 4;321(5885):48–50.

58. Rabinovich MI, Afraimovich VS, Bick C, Varona P. Information flow dynamics in the brain. Physics of life reviews. 2012 Mar 1;9(1):51–73. pmid:22119154

59. Meyer JS, Gotoh F, Favale E. Cerebral metabolism during epileptic seizures in man. Electroencephalography and Clinical Neurophysiology. 1966 Jul 1;21(1):10–22. pmid:4165360

60. Jafarian A, Zeidman P, Wykes RC, Walker M, Friston KJ. Adiabatic dynamic causal modelling. NeuroImage. 2021 Sep 1;238:118243. pmid:34116151

61. Kiebel SJ, Garrido MI, Moran RJ, Friston KJ. Dynamic causal modelling for EEG and MEG. Cognitive neurodynamics. 2008 Jun;2(2):121–36. pmid:19003479

62. Creutzfeldt OD. Generality of the functional structure of the neocortex. Naturwissenschaften. 1977 Oct;64(10):507–17. pmid:337161

63. Ising E. Beitrag zur theorie des ferro-und paramagnetismus [Doctoral dissertation]. [Hamburg]: Grefe & Tiedemann; 1924.

64. Potts RB. Some generalized order-disorder transformations. In Mathematical proceedings of the cambridge philosophical society 1952 Jan (Vol. 48, No. 1, pp. 106–109). Cambridge University Press.

65. Jirsa VK, Friedrich R, Haken H, Kelso JS. A theoretical model of phase transitions in the human brain. Biological cybernetics. 1994 May;71(1):27–35. pmid:8054384

66. Proix T, Jirsa VK, Bartolomei F, Guye M, Truccolo W. Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy. Nature communications. 2018 Mar 14;9(1):1–5.

67. Roy D, Sigala R, Breakspear M, McIntosh AR, Jirsa VK, Deco G, et al. Using the virtual brain to reveal the role of oscillations and plasticity in shaping brain’s dynamical landscape. Brain connectivity. 2014 Dec 1;4(10):791–811. pmid:25131838

68. Haken H. Synergetics: An introduction. Non-equilibrium phase transition and self-selforganisation in physics, chemistry and biology. Berlin: Springer Verlag;1983

69. Frank T. D. Nonlinear Fokker-Planck Equations: Fundamentals and Applications. Springer Series in Synergetics. Berlin: Springer; 2004.

70. Haken H, Portugali J. Information and selforganization: A unifying approach and applications. Entropy. 2016 Jun;18(6):197.

71. Kelso JS. Unifying large-and small-scale theories of coordination. Entropy. 2021 Apr 27;23(5):537. pmid:33925736

72. Ellis GF, Noble D, O’Connor T. Top-down causation: an integrating theme within and across the sciences?. Interface Focus. 2012 Feb 6;2(1):1–3.

73. Lopes da Silva FH, Blanes W, Kalitzin SN, Parra J, Suffczynski P, Velis DN. Dynamical diseases of brain systems: different routes to epileptic seizures. IEEE transactions on biomedical engineering. 2003 May 13;50(5):540–8. pmid:12769430

74. Goodfellow M, Schindler K, Baier G. Intermittent spike–wave dynamics in a heterogeneous, spatially extended neural mass model. NeuroImage. 2011 Apr 1;55(3):920–32. pmid:21195779