Animal experiments in this study were conducted in accordance with the U.K. Animal (Scientific Procedures) Act 1986, with approval from Home Office (license PPL70-13691) and the local ethics committee at the Institute of Neurology, University College London. The data were recorded in a chemoconvulsant animal model of seizures (using 4-aminopyridine drug) from awake, unanaesthetised mice. Recently developed state-of-the-art graphene-based transistor arrays were implemented for wide-bandwidth electroencephalography recordings from the surface of cortex (epicortical grids of 4 × 4 which are placed over the somatosensory cortex), as well as laminar recordings through a cortical column, in the visual cortex (intracortical depth electrode with 14 recording sites) (see Bonaccini Calia et al. (2021) for further details).

Animals were group-housed (to acclimatize for at least 1 week before surgery) on a 12-h/12-h dark/light cycle, where food and water were given ad libitum. A surgery was performed to implant a head bar on the mouse skull to allow stable attachment of the animal to a Neurotar chamber for reliable (and effectively movement artifact free) data recordings both from epidural and intracortical columns with a sampling frequency of 9.6 kHz (Bonaccini Calia et al., 2021). Focal injection of 4-aminopyridine (4-AP) (50 mM; 350 nL into the somatosensory cortex region) was performed at a depth of ~500 μm into the cortex using a 33-gauge needle (see Figure 1).

The epi-/intracortical electrodes capture the neuronal activity during the entire experiment. The spread of drug is locally restricted to the site of injection (Rossi et al., 2017). For the aim of this study, we use data from the superficial layers of the cortex located close to the injection site. The trace of the wide-bandwidth and high-pass-filtered data is shown in Figure 1. An interesting feature in these data is that some seizures are accompanied with DC shifts and some without DC shifts. Therefore, a key hypothesis for the data in this study is that there could be differences between biological generators of seizures that are accompanied with or without DC shifts.

A neural mass model (NMM) describes the electrical activity of a cortical column (which is captured using electroencephalography techniques) through a low-dimensional biological informed dynamical systems (Wilson and Cowan, 1972; Jansen and Rit, 1995). Assumptions that contribute to driving NMM formulation are as follows: (i) synchronized firing of neurons in a cortical column induce observable electrical activity (Mountcastle, 1957; Hubel and Wiesel, 1963; Felleman and Van Essen, 1991), (ii) neurons within cortical columns can be clustered into few populations (because at each layer of the cortex, effectively one type of neurons resides), mean activity of each of which could be modeled independently, and (iii) the mean activity of different populations interacts and shape the mesoscale electroencephalogram recordings (this assumption is supported by the statistical mean field theory) (Cowan, 1969; Deco et al., 2008; Faugeras et al., 2009).

NMMs are minimally biologically informed dynamical systems that approximate interactions between (collective) activities of neurons in a cortical column with few neuronal populations. NMMs are low-dimensional, with few biological parameters (compared to networks of interconnected neurons) and can well recapitulate electrophysiological data (Jansen and Rit, 1995). Crucially, there are formal mathematical links between an interconnected network of single neurons and an NMM (Deco et al., 2008; Faugeras et al., 2009; Veltz and Faugeras, 2010; Faugeras and Inglis, 2015). Collective dynamics of neurons are well approximated by NMMs, while they retain biological realism. The low-dimensional representation of NMMs is well suited for either designing biologically motivated control systems to suppress seizures (Schiff, 2011) or understanding interactions between different brain regions in cognitive tasks (Friston et al., 2011, 2019; Schwartenbeck and Friston, 2016; Shaw et al., 2017; Jafarian et al., 2019c). Therefore, NMMs are suitable for a variety of translational neuroscience applications and can be implemented to study biologically motivated hypothesis regarding underlying generators of brain activity.

One could show (both theoretically and experimentally) that dynamics of a population of many neurons can be governed (or summarized) by static conversion of input mean synaptic activity to firing rates (average of action potentials) and firing rates to mean synaptic activity (scaled by anatomical connection strengths between populations or cortical layers), which is input to other populations (Wilson and Cowan, 1972; Freeman, 1975; Jansen and Rit, 1995).

Mathematically, synaptic potential (V) into firing rate conversion is modeled by a sigmoid transformation (σ(.)) as previously described (Wilson and Cowan, 1972).

In Equation (1), e₀ is the maximum firing rate, ρ is the slope of the transformation, and V_th is the firing threshold (when input potentials reach half maximum firing rates).

The conversion of the firing rate to the mean postsynaptic membrane potential is modeled by a second-order low-pass filter with an impulse response h(t)=AtTe-t T where t≥0, A is a maximum postsynaptic potential (also known as synaptic gain), and T is a synaptic time constant (Freeman, 1975). Postsynaptic potential V that is generated by firing rate σ(.) can be calculated simply by convolving the input firing with synaptic kernels, V(t) = h ⊗ σ (the symbol ⊗ represent the convolution), which is equivalent to the following second-order differential Equation (note the second-order differential equation can be equivalently written in terms of two first-order systems)

The input–output relations that govern the dynamics of a neuronal population are largely similar among different NMMs, although the number of neural populations and their connections (scale factor to synaptic potentials) to each other can be varied between NMMs. In this study, we use an NMM that was developed by Jansen and Rit (1995) and is shown in Figure 2.

This NMM explains the activity of a cortical column through interactions of three populations, namely, (i) excitatory spiny stellate cells (denoted by ex in model equations) situated in layer 4 of a typical cortical column, and they receive endogenous input from other regions; (ii) inhibitory cells (denoted by i in model equations) that are distributed across the columnar organization; and (iii) pyramidal cells (denoted by p in model equations) whose activity is predominantly captured by recording electrodes. The equation of the NMM associated with the structure in Figure 2 can be written as follows:

In Equation (3), the vector [xp, xp*, xex, xex*,xi, xi*]∈ R6 represents the physiological states, either the synaptic activity of the populations (x_p, x_ex, x_i) or their first-order derivative (xp*,xe*,xi*). The fixed parameters in Equation (3) are synaptic parameters for each population (A.,1T.), inter-region connections which are defined by scaling the universal constant C, and parameters of sigmoid transformations. The endogenous input to the model is ω, which can be a random, constant, or smooth function (see Table 1 for the range of variations for NMM parameters). The membrane potential of pyramidal cells is considered simulated electroencephalogram data. Hereinafter, we rewrite the NMM equation in the form of a general canonical dynamical systems as follows:

In Equation (4), x= [xp, xp*, xex, xex*,xi, xi*]∈ R6 and represents the biological states, the right hand side of equation 3 is denoted by a nonlinear function f(.)∈ R6, the vector n(t) is the endogenous fluctuations (which can be modeled as a constant number, smooth function, or noise), and θ is the set of constant parameters in the model (thereby their derivative is zero as written in the second line of the Equation (4).

The most obvious way to explore the origin of the oscillations in the brain using the NMM is forward simulation of different sets of parameters and initial conditions. For instance, in Jansen and Rit (1995), the model was simulated for different values of connections between populations (the universal connection constant C was varied from 68 to 675) and brain rhythms, such alpha and spike wave discharges, were replicated. The NMM is used to explore a path between normal and pathological activity in parameter space by changing the balance between maximum postsynaptic responses of excitatory and inhibitory connections (e.g., Wendling et al., 2002, 2005). Motivated by these studies and electrophysiological knowledge about the animal model of seizures in this study, we assume that an underlying cause for the transitions can be explained by changes of postsynaptic potentials of inhibitory populations (A_i).

Sophisticated mathematical analysis such as bifurcation can be performed to formally study the behavior of the NMM as some parameters (one or two) are changed (e.g., Grimbert and Faugeras, 2006; Spiegler et al., 2010; Touboul et al., 2011). Despite valuable information informed from bifurcation analysis, practically, it can only be applied to one of two parameters. In addition, this form of analysis can only be applied to a deterministic NMM (i.e., the input to the model is a constant). In other words, the interpretation of the bifurcation structure of random bifurcation is still in its infancy (e.g., Crauel and Flandoli, 1994, 1998; Callaway et al., 2017). However, it is possible to prove that the solution of the NMM is stationary for a constant set of parameters, with stationary noise as its input (Faugeras et al., 2009; Veltz and Faugeras, 2010, 2011; Faugeras and Inglis, 2015). In this study, we re-confirm these findings by simulating the global invariant measures of the NMM for different parameters to illustrate that only one type of activity can be viewed using an NMM with fixed parameters.

In this study, first, we estimate constant parameters of the model based on spectral features in a stationary segment of data using dynamic causal modeling (also known as variational Bayesian inversion of a nonlinear system under Laplace assumption) of spectral response (Friston et al., 2012; Jafarian et al., 2021). The assumption in dynamic causal modeling of cross-spectral response is that neural dynamics rest at a stable equilibrium and oscillations (finite deviation from baseline equilibrium) are induced due to random exogenous input (Lopes Da Silva et al., 1974; Friston et al., 2012). [See Friston et al. (2012) for details of procedures].

We used the ensuing identified model and tracked changes of postsynaptic potentials of the inhibitory population in an NMM using a continuous–discrete unscented Kalman filter (UKF) (Voss et al., 2004; Jafarian et al., 2019b) from data that exhibit transitions into and out of seizures. The UKF is a class derivative-free stochastic filter that can recursively estimate hidden states of partially observed nonlinear dynamical systems from real data. The UKF approximates posterior estimates of states using Gaussian distribution, which in turn makes this filter computationally efficient while accurate (Sitz et al., 2002; Voss et al., 2004; Sarkka, 2007).

The generative model for tracking the inhibitory gain in the NMM for a given data using the unscented Kalman filter (UKF) is as follows:

The first line of Equation (5) is similar to the general equation of the NMM. The second line expresses the dynamics of inhibitory gains, which is recovered from data using the UKF method. We consider additive uncertainty (with very small variations n * ~N(0,10-8) to the equation of motions of θ(= A_i) to allow the UKF algorithm to adjust the values of the synaptic gains from the data (Schiff, 2011; Jafarian et al., 2019a). It should be noted that in the context of state estimations from data using any form of Kalman filter, the role of noise in hidden states is interpreted as uncertainty (Sitz et al., 2002; Voss et al., 2004). The left hand side of the third line of Equation (5) is discrete real data (sampled recordings), and the right hand side expresses how the solution of the generative model is linked to observational data (here, the membrane potential of pyramidal cells which is defined by the vector H = [0 0 1 0 −1 0 0]'), and e is the random effect which has a normal distribution.

We performed forward simulation of a noise-driven NMM akin to Jansen and Rit (1995) for different values for intrinsic connections (C = 68, 128, 135, 270, 675) and replicated alpha rhythms and spike wave discharge, as shown in Figure 3. NMM is simulated using a stochastic RK method as explained in Wilkie (2004). In Figure 3, we also show the global attractor of the NMM for alpha rhythm and spike wave activity. The shape of the global attractor illustrates the region in the model phase space that is almost certainly occupied by the noise-driven model in the limit. Because the phase space that is occupied by the global attractor is similar to its finite time simulations, one could assume that the noise-driven model for different fixed parameters could generate alpha and epileptic discharge without any transitions. Finding the global attractor can be seen as a complementary analysis to recent mathematical proof that the noise-driven NMM could generate stationary solutions with fixed parameters (Faugeras et al., 2009; Veltz and Faugeras, 2010, 2011). Here, we also provide a local bifurcation plot of the NMM, similar to Grimbert and Faugeras (2006), which can be used to study the system equilibrium properties as its input is altered. As mentioned, the conventional bifurcation analysis can be applied to deterministic systems and could reveal information relevant to the underlying mathematical features that support different sorts of activities. For instance, in the bifurcation plot in Figure 3, one could argue that limit cycle activity may be associated with epileptic activity.

To simulate transition into and out of seizures, we equipped the NMM with slow dynamics of postsynaptic gain of inhibitory populations, as shown in Figure 4. In this simulation, the inhibitory gain is regulated by input firing rates and also has a period of recovery time. This effectivity resembles alteration of synapse gain and plasticity (e.g., Jafarian et al., 2021). As can be seen in Figure 4, this slow–fast NMM can show transitions to and from seizures, as well as alterations in DC shifts as the plasticity of inhibitory populations is altered.

We infer baseline parameters (noise inputs and synaptic gains), while others are fixed (see Table 2 for details) from a stationary segment of real data using dynamic causal modeling of cross-spectral density (Friston et al., 2012). We assume that all ensuing baseline parameters remain unchanged during transitions to seizures in data, except the synaptic gain of the inhibitory population. This effectively models pathological alterations of excitatory-inhibitory balance during paroxysmal transitions (Wendling et al., 2002).

We employed an unscented Kalman filter to track changes of inhibitory gain in two sets of real data (scaled by a factor of 10⁴ to make its variance in the range of the NMM output) and their high-pass-filtered version (which mimics data that are usually available from conventional electrographic devices). We start by tracking the inhibitory gains for the full-bandwidth data with a small DC shift. This effectivity implies that the high-pass-filtered version and original wide-bandwidth data are very similar. The outcome of tracking the dynamics of inhibitory gain is shown in Figure 5 for both wide-bandwidth and high-pass-filtered data. The inferred trajectories of inhibitory gain have very similar behavior. In this case, the inferred inhibitory gain is altered as the large spikes appear in the data and increase during the recurrence of seizures.

We repeat the analysis for the wide-bandwidth data with significant DC shifts and for the same-signal subsequently high-pass-filtered data. The inferred inhibitory gains for both wide-bandwidth and high-pass-filtered data are shown in Figure 6. In this simulation, there are significant differences between behavior (effectively in the opposite direction) of the recovered slow evolution of inhibitory gains from both data. In the high-pass-filtered version of the data, the synaptic gain increases before paroxysmal transitions and returns to its baseline after seizures. The recovered trajectory of the inhibitory gains from wide-bandwidth data shows a totally different behavior as it starts to decrease during the pathological activity and then increases during seizures.

These four examples may be seen as a proof of concept that all information in the recorded data needs to be considered for capturing the underlying dynamics of biological generators of data, particularly when modeling data with paroxysmal transitions into and out of seizures, which requires information about all frequency ranges from slow to very high.