Physics of brain dynamics: Fokker–Planck analysis reveals changes in EEG δ–θ interactions in anæsthesia

We use drift and diffusion coefficients to reveal interactions between different oscillatory processes underlying a complex signal and apply the method to EEG δ and θ frequencies in the brain. By analysis of data recorded from rats during anæsthesia, we consider the stability and basins of attraction of fixed points in the phase portrait of the deterministic part of the retrieved stochastic process. We show that different classes of dynamics are associated with deep and light anæsthesia, and we demonstrate that the predominant directionality of the interaction is such that θ drives δ.


Introduction
Complex signals typically contain a great deal of information about the underlying processes that generate them, but the extraction of such information can be difficult and requires methods drawn from physics and nonlinear science [1]- [4]. This is especially true of signals derived from living systems. Their complexity arises in part from nonlinear interactions between subsystems oscillating on widely differing timescales, and in part from the confusing influence of seemingly random fluctuations. Electroencephalographic (EEG) recordings of brain activity [1] exemplify the problem. It has long been appreciated [5] that they exhibit continuous widespread oscillations, but their origin remains barely understood. As in the case of cardiovascular signals [6], there are multiple rhythms enclosed within broad-band noise. In dynamical terms, these oscillations might be attributed to limit-cycle attractors, because analyses of short segments often reveal spectral peaks in particular frequency ranges. Interactions between the processes generating some of these EEG waves have been demonstrated [7] but not yet explored much.
Understanding the interdependences between complex dynamical systems is of quite general scientific importance, e.g. in physics, chemistry, biology and neuroscience; it can also be of importance for, e.g. economics and sociology. Often, the underlying equations of motion are unknown, but a detailed quantitative description of interdependences can nonetheless be achieved by time-series analysis of experimentally measured observables. In recent years, a number of methods have been developed that allow one to detect and to quantify the strengths of the possible interactions [8]. Asymmetric approaches have facilitated the detection of directional coupling from time series. At some risk of over-simplification they can conveniently be divided into three main groups: techniques based on interrelationships between parameters of oscillatory dynamics and, in particular, phase dynamics [9]- [11]; state-space-based methods [12,13] and information-theoretic approaches [14,15].
State-space-based approaches provide strong methods for detecting generalized synchronization or directionality patterns between two coupled systems. Also, the fact that there are not many free parameters involved in driving the state-space patterns makes such measures relatively easy to implement, and robust. Their main challenges are in finding welldefined quantitative measures of directionality or driver-response and in that it is often difficult to localize the reconstructed properties in time. Generally, all three groups of methods tend to make rather strict assumptions about the dynamics of the systems that generate the time series, e.g. that they should be linear systems, or in the case of self-sustained oscillators that they should be weakly-coupled. Furthermore, many of the approaches focus preferentially on the low-dimensional deterministic part of the dynamics. In reality, however, the influences of noise on nonlinear dynamical systems can be far from trivial [16]- [18] and they constitute a subject of continuing interest and importance. Real signals frequently result from dissipative dynamical systems under the influence of noise, and it is often useful to describe them in terms of Fokker-Planck formalism [19].
In this paper, we combine the multidimensional Fokker-Planck equation [16,20] with mathematical modelling of the system dynamics, enabling us to extract model parameters directly from the measured time series [2,3,20]. We concentrate on the interactions between EEG θ waves and δ waves in anaesthesia, and characterize some of their properties. Following [20,21] we apply a two-dimensional (2D) Markov analysis to EEG signals recorded [22] from rats during deep and light anaesthesia (see examples in figure 1); and we carry out a stability analysis of the deterministic part of the derived equations. We find in these experiments that the dynamics falls into different universality classes for the states of deep and light anaesthesia, enabling us to distinguish robustly between them. Although we treat brain dynamics explicitly, the method is applicable to complex systems quite generally.
In section 2, we describe succinctly how our experimental data were acquired from rats and pre-processed prior to analysis. Section 3 describes the basic method used for the analysis, and section 4 discusses how we applied it to the rat EEG signals and considers the results that emerge. Section 5 describes some simple tests undertaken to eliminate artifacts and to check the validity of the results. Finally, in section 6, we summarize and consider the main results, and draw conclusions.

Experimental methods and data acquisition
A full description of how the signals were measured has already been given elsewhere [22], so here we just summarize the salient features. The experiments were performed on 10 adults, male Wistar rats weighing 250-300 g. The animals were each anaesthetized with a single intraperitoneal injection of ketamine hydrochloride (45 mg (kg body wt) −1 ) and xylazine hydrochloride (7 mg (kg body wt) −1 ). (Measurements were also carried out on a second group of rats anaesthetized with pentobarbital [22] but these are not discussed here.) As soon as the rat could no longer hold its upright posture, 10-15 min after administration of the drug, it was placed in a darkened Faraday cage where sensors were mounted and recording started immediately. The EEG was recorded over the left and right parietal cortex with three hypodermic needles inserted under the animal's scalp to serve as electrodes. The EEG was differentially amplified by 10 3 × and low-pass filtered with a cut-off frequency of 300 Hz. The resultant signal was fed through a signal conditioning system, digitized at 1 kHz with 16-bit resolution, and stored on the hard disk of a laptop computer.
Depth of anaesthesia was assessed at 5 min intervals by a nociceptive stimulus, the skinpinch test, applied to the sole of the animal's front paw [23]. The recording started with a negative test response, i.e. when the rat stopped responding with a reflex withdrawal of the limb. The monitoring was terminated on the reappearance of a positive pinch-test response, as the animal immediately started to move thereby terminating reliable data recording. The duration of recording varied from rat to rat and was on average 87 min. The measurements were at a constant room temperature of 24 ± 1 • C.
The transition from deep to light anaesthesia was determined by application of several criteria: as the times when both cardiac and respiratory frequencies significantly increased and became significantly more variable [22]; as the times when the amplitude of δ-waves dramatically decreased; and as the times when the amplitudes of θ and γ significantly increased. These different criteria yielded consistent results. The transition is most probably associated with the effects of the ketamine wearing off faster that those of xylazine so that in the light phase the effect of xylazine is dominant.

Fokker-Planck analysis
We define a state vector q = {δ(t), θ(t)} composed of the signals of the δ and θ bands, hypothesize that it follows a stochastic process of form and approximate the fluctuations j (t) by Gaussian white noise. We suppose that the process is Markovian, an assumption that can subsequently be validated by data analysis. The drift vector D (1) describing the deterministic part of the dynamics, and the diffusion matrix D (2) determining the strength of the driving noise, are We have used the Itô interpretation of the stochastic integral and conditional expectation values · that can be determined from the experimental data; √ D (2) is to be calculated by diagonalizing the matrix D (2) , taking the square root of each of its elements, and transforming the result back into the original system of coordinates [16]. The Markov condition requires an n-point conditional probability function (CPF) P(q(t + nτ )|q(t + (n − 1)τ ), . . . , q(t + τ ), q(t)) to be equal to a two-point CPF P(q(t + nτ )|q(t + (n − 1)τ )). The timescale τ in which this condition is fulfilled is the Markov timescale (τ M ). For multidimensional stochastic variables, it is hardly possible to check the Markov condition directly by means of n-point CPFs. We therefore use three-point CPFs, for which the Markov condition corresponds [25] to being zero. Here, and x represents one of the components of q, i.e. q i . σ 2 J and σ 2 M are the variances of P(x 3 , x 2 , x 1 ) and P(x 3 |x 2 )P(x 2 , x 1 ), respectively. Figure 2 shows χ 2 calculated for both δ and θ waves from a rat EEG signal as a function of the timescale τ . In both cases, the minimum value of χ 2 occurs just after 0.1 s, after which it hardly changes, corresponding to τ M with one σ confidence level.
Using equation (1) with this value of τ M , we can estimate the drift and diffusion coefficients from δ(t) and θ (t) as shown in figure 3. In doing so, we used 1D coefficients, i.e. in deriving the drift and diffusion coefficients for each variable, the effect of the other variable was integrated. Both coefficients can be well-approximated by low-order polynomials. The D (1) i with i = δ, θ indicate an overall damping behaviour, as shown earlier [26]. Moreover, the behaviour of D (2) δδ and D (2) θθ implies that the noise is multiplicative in character: constant diffusion coefficients represent additive noise whereas, for multiplicative noise, the diffusion coefficients are functions of the dynamical variables. We tested whether the deviation from additive behaviour is attributable to the finiteness of τ M [27] by consideration [28] of a Taylor expansion of the drift coefficients. We conclude that the effect of the higher order terms in τ M cannot explain the apparently multiplicative nature of the estimated second coefficient, which implies that the noise strengths are functions of the dynamical variables.   Figures 3(a) and (c) show that there are significant differences in the δ waves between deep and light anaesthesia, consistent with earlier work [22]. On the other hand, panels (b) and (d) show that little change occurs in the θ waves. The functionality of D (2) δδ changes markedly between the two states as can be seen in figure 3(c): in the light phase, a parabola fits this coefficient very well, while for deep anaesthesia we must use a fourth-order function to describe the central peak.

Stationary solution for retrieved drift and diffusion coefficients
We can also examine qualitatively the stationary solution of the Fokker-Planck equation around the different extrema in D 2 (δ). The solution is P stationary = N 0 D 2 (δ) exp(− dδ D (1) D (2) ), where the integration constant N 0 can be derived by normalization of P stationary . For simplicity, we suppose that D (1) δ is well approximated by a line (D (1) δ = −bδ where b > 0). In light anaesthesia, D (2) δδ has a minimum at δ = 0, whereas in the deep state there is a local maximum there. Focusing on the region around this point, D (2) δδ can be approximated as D (2) δδ D (2) δδ (0) ∓ 1 2 |D (2) δδ |δ 2 , where denotes d/dδ. It is evident that the contribution of D (2) δδ to the stationary solution of the Fokker-Planck equation changes between deep and light anaesthesia. Two additional extrema (minima) appear in D (2) δδ during deep anaesthesia ( figure 3(c)), causing a double-humped PDF in deep anaesthesia (cf the single-humped PDF during light anaesthesia).
To test whether the behaviour of the system deviates from a Langevin description, we calculated the fourth-order coefficients D (4) 4 for x = δ and θ , allowing us to determine whether the driving noise process j (t) exhibits deviations from a Gaussian distribution [16]. Only if D (4) vanishes is j (t) white Gaussian. The probability density function (PDF) of the process then evolves according to a Fokker-Planck or, equivalently, Langevin equation [16]. Figure 3 shows that D (4) δ is slightly above zero, but the magnitude of this coefficient is less than one-fifth of the second coefficient D (2) δδ , for both deep and light anaesthesia. The fluctuations for the case of light anaesthesia are larger than for the deep state. This can be seen from the observation that, in the domain of [−3σ, 3σ ] in figure 3(c), the D (4) δ terms are bigger in the light case, which indicates greater fluctuation amplitude. In contrast, for the θ band signals, D (4) θ compared with D (2) θ θ was clearly above zero for both deep and light anaesthesia suggesting that a 1D Langevin description may be inadequate for modelling δ-θ activities.
In a 1D Langevin model, the effects of other components are integrated, e.g. the 2D drift and diffusion coefficients can be integrated to yield the 1D ones as D (1) i Obviously, the information about interactions and their directionality is lost. This is an important point: although we have clearly different functionality for D (2) δδ and a different pattern for D (1) δ , the differences are unclear in the case of the θ component and we should not thoughtlessly integrate out the effects of the other variable. So, next, we derive the drift and diffusion coefficients as functions of both δ and θ . Figure 4 shows all the coefficients calculated for deep anaesthesia. Figures 4(a) and (b) plot D (1) δ (δ, θ ) and D (1) θ (δ, θ). The θ -component of the drift vector shows an almost linear dependence on θ and only weak variations with δ. But interestingly, the δ-component of the drift vector has strong functionality in both δ and θ (see figure 4(a)). We used a 2D least-squares method to derive the functionality of the coefficients. In the case of θ, we have D (1) θ (δ, θ ) = a + bδ + cθ with good precision and the value of b is small (|c/b| 10 for all rats), while in the case of δ the situation changes dramatically. The best fit available is a third-order polynomial with respect to δ and θ, D (1) δ (δ, θ) = m,n a m n δ m θ n (m + n = 0, 1, 2, 3). In figures 4(c) and (d), we plot the coefficients D (2) δδ and D (2) θ θ of the diffusion matrix. It is evident that D (2) θθ , like D (1) θ , depends weakly on δ. In both cases, a second-order polynomial of δ and θ such as m,n a m n δ m θ n (m + n = 0, 1, 2) can describe the functionality of the two diffusion components. In case of light anaesthesia, we find similar functionality albeit with different parameters for the diffusion coefficients; but we find different functionality for the δ-component of the drift coefficient, which has a much weaker dependence on θ . It was found that these same patterns were repeated for all ten rats. Thus, they can be used to characterize the dynamical behaviour of the δ and θ brain activities.
The fact that D (1) δ and D (2) δδ have strong dependence on both δ and θ , whereas the D (1) θ and D (2) θθ hardly vary at all with δ, shows that the EEG θ component influences the δ-component but not vice versa. This situation persists regardless of depth of anaesthesia.  and D (2) δδ have strong dependences on both δ and θ , whereas the dependence of D (1) θ and D (2) θ θ on δ is weak.

Dynamics and fixed points (FPs)
To complete the analysis we need to derive a more quantitative description of the patterns appearing in the drift and diffusion coefficients. The generated Langevin equations give the deterministic and stochastic parts of dynamics, which can be written as ∂ t q i = ∂ D t q i + ∂ S t q i , of which ∂ D t q i = D (1) i is the deterministic part of the dynamics and ∂ S t q i = 2 j=1 D (2) i j j . Because of the finite noise we do not have clear control of the stochastic term so, for simplicity, we just compare the dynamics arising from the deterministic terms: we can derive the FPs and the stability exponents in the phase diagram from the deterministic parts of the δ and θ trends. We note that another way to determine asymmetries in coupling is from measures based on the drift and diffusion coefficients [29].
In carrying out these analyses, we used all data we had for each rat [22], i.e. for a measurement of duration 30 min we had 30 × 60 × 1000 = 1.8 × 10 6 sample data for our sampling frequency of 1000 Hz. Thus the results related to D (1) , D (2) , and D (4) are all supported by a good level of statistics. However, we have also checked our results by using different portions of data. We found that 5 min worth of data (300 k data samples) was sufficient to give us a robust estimation of D (1) , D (2) and D (4) . In relation to the Markov length scale τ M we found that, although the results were obtained for τ M = 0.1 (100 sample data), a range of 80-120 samples gave the same results. Diverging from this range, however, would drive us into a regime where modelling based on the Markov process is no longer valid.
To examine the FPs and flows in the δ-θ phase space, we have to solve together: for both deep and light anaesthesia. We find that there are two distinct regimes. Firstly, for eight out of ten rats during deep anaesthesia there are three nontrivial FPs. Two of them are stable in both the δ and θ directions, and the third one is unstable in the δ but stable in the θ direction. Figure 5(a) shows the flow diagrams related to equation (3) for a typical rat. There is a region outside the two stable FPs in δ-θ space in which the flows take any initial point to one of the stable FPs. However, when δ or θ are small enough that the initial point is close to the unstable fixed point, a small change in δ could take the flow to a different stable fixed point, hence implying that the joint PDF of the δ and θ , P deep (δ, θ), must have two distinctive peaks. Secondly, for all ten rats during light anaesthesia, we have only one attracting FP. Figure 5(b) shows that, for a typical rat, over the whole interval of δ and θ flows take any initial point to the stable fixed point. Hence, in this case, the system appears homogeneous at large timescales and we must have one peak for P light (δ, θ), unlike the previous case. To test this prediction, we plot in figures 5(c) and (d) the directly calculated PDFs P deep (δ, θ) and P light (δ, θ). They are in full agreement with the results of our stability analysis of drift coefficients. At least in the present experiments, therefore, we may conclude that the dynamics of δ-θ activity falls into different universality classes for the states of deep and light anaesthesia. However, this finding needs further investigation and will be a subject of our future work.

Tests of the validity of the results
We now describe some tests that we carried out to confirm that our results relate to the real dynamics of the system, and are not just artifacts of filtering. To do so, we processed a dummy signal consisting of Ornstein-Uhlenbeck noise [30] by filtration in the δ and θ bands and then repeated the whole analysis procedure. We used standard Ornstein-Uhlenbeck noise with an added sinusoidal interaction to simulate the external oscillatory components. The equation for where A corresponds to the coupling strength of the external oscillatory components, ξ is the noise strength, and η is the white Gaussian noise. We always obtained one attracting FP, implying the damping structure of the system. Moreover, when analysing the empirical data with 5 min moving windows for the δ and θ bands of the real signals, we observed that the change of pattern between deep (3 FPs) to light (1 FP) anaesthesia occurs very sharply: we estimated the D (1) , D (2) for one window of data, and then repeated the whole process explained in sections 4.1 and 4.2 by first deriving the functionality of D (1) and D (2) based on the dynamical variables, and then applying stability analysis to the deterministic parts of the processes (similar to the procedure explained in section 4.2). Interestingly we found that even with the short 5 min data windows, we could distinguish clearly between two different phases, having 3 and 1 FPs for deep and light anaesthesia, respectively, the transition from one to the other being related to the transition from deep to light anaesthesia, which occurs very sharply (∼1-2 min).
Finally, we investigated the question of whether or not the same results can be derived from the phase dynamics of the variables δ and θ . We used the Hilbert transform [31] to detect the corresponding phases; we calculated the time interval for each period of the phase; and we inverted these time series to derive new signals in frequency space. Note that the phases of δ and θ are defined only in a statistical sense because both waves are relatively broad in frequency. It is only δ waves in deep anaesthesia that are relatively localized in frequency. For further details see [22].
We then processed the δ and θ phases using the Fokker-Planck formalism. Reassuringly, but unsurprisingly, there were no significant differences between the drift and diffusion coefficients corresponding to the δ and θ phases in deep and light anaesthesia (not shown here).

Concluding remarks
In conclusion, our Fokker-Planck analysis of EEG signals has provided evidence of interactions between the EEG δ and θ wave activities in rats in deep and light anaesthesia. We have characterized the functionalities of δ and θ and shown that the functionality of δ alters with changes in depth of anaesthesia. It is commonly accepted that sleep or anaestheticinduced unconsciousness arises through changes in the conduction of ion channels. Ketamine inhibits [32] NMDA 4 and Muscarinic ACh-sensitive 5 channels, which may occur through two distinct mechanisms: (i) ketamine blocks open channels thereby reducing their mean open time; or (ii) it decreases the frequency of channel opening. Either mechanism would lead to neuronal hyperpolarization and reduce the activation of action potentials. The consequence of the hyperpolarization is that the function of neurons is blocked. Where exactly this occurs is currently unknown. Specific suppression of activity in the regional-thalamic and midbrain reticular formation and hyperpolarization blockage of thalamocortical neurons has been discussed [33]. On the other hand, ketamine is considered as a drug that does not reduce cortical metabolism and glutamate release or depress sensory information flow through the thalamus [34]. Our observation that θ drives δ, regardless of depth of anaesthesia is consistent with this picture. It is an important result, suggesting that the changes in δ-activity during deep anaesthesia or sleep are mainly due to reduced sensory information from the spinal cord to the thalamus, rather than on account of reduced thalamocortical interactions. However, the physiological and neurological implications of this finding remain to be further investigated.