The anesthetic propofol shifts the frequency of maximum spectral power in EEG during general anesthesia: analytical insights from a linear model

1. Introduction

General anesthesia (GA) is an important medical application in today's hospital surgery. Although GA is omnipresent in recent medicine, its underlying neural interactions have been a long-standing mystery. In the last decades, the anesthetic phenomena have attracted theoreticiens, e.g., (Steyn-Ross et al., 1999; Bojak and Liley, 2005; Hutt, 2011), who aim to describe mathematically some major experimental phenomena by population models (Steyn-Ross et al., 2004, 2012; Bojak and Liley, 2005; Hutt and Longtin, 2009; Hindriks and van Putten, 2012), or spiking-neuron models (McCarthy et al., 2008; Ching et al., 2010). Most theoretical studies aim to explain signal features of electroencephalographic (EEG) data observed during anesthesia. Such features comprise the diminution of α-activity accompanied by a subsequent enhancement of δ-activity while increasing anesthetic concentration (Gugino et al., 2001; Cimenser et al., 2011; Murphy et al., 2011) and the power enhancement of activity induced by some anesthetics (McCarthy et al., 2008; Ching et al., 2010). Another example is the increase of the frequency of maximum EEG-power to higher values as observed experimentally in several studies (Gugino et al., 2001; Ching et al., 2010; Murphy et al., 2011; Boly et al., 2012; Hindriks and van Putten, 2012). The current work focusses on the power enhancement and the frequency shift of maximum power while increasing the anesthetic concentration and gives insights into its origin by the analytical treatment of a linear neural field model.

One of the major objectives of this work is to answer the question whether it is possible to explain spectral EEG-features observed during GA by a low-dimensional linear model. The advantage of such a reduced model is the analytical tractability and an identification of underlying neural interactions or even the origin of the spectral feature. Here the difficulty is to find a simple model, that, however, still involves important, i.e., realistic and neural interactions. We are convinced that such a model has been found in a previous work (Hutt and Longtin, 2009). The present work will simplify further this spatio-temporal model while taking into account the biophysical effects of the anesthetic propofol on synaptic receptors and hence retaining the neurobiological plausibility.

The simplicity of the model will allow to reveal the effect of different actions of the anesthetic propofol on synaptic receptors on the frequency of maximum spectral power. Moreover, the work gives criteria under which conditions the frequency of maximum power increases with increasing propofol concentration and when it may decreases. In the analytical treatment, we will see that the power enhancement for larger propofol concentrations may be explained by an oscillatory instability and we predict that it always occurs while the frequency of maximum power increases.

2. Methods

2.1. The Model

The neural field model under study (Hutt and Longtin, 2009) describes the evolution of the mean membrane potential of a neural population in a small spatial patch at spatial location x and at time t. Similar models have been derived and studied before (Wilson and Cowan, 1972; Amari, 1977; Ermentrout, 1998) and applied successfully to explain spatio-temporal neural activity observed experimentally (Ermentrout and Cowan, 1979; Huang et al., 2004; Angelucci and Bressloff, 2006; Schwabe et al., 2006). The population includes both excitatory and inhibitory neurons and takes into account excitatory and inhibitory synapses. Assuming that excitatory and inhibitory neurons exhibit identical effective membrane potentials, i.e., an identical difference between excitatory and inhibitory post-synaptic potentials, the mean excitatory and inhibitory post-synaptic potentials V_e(x, t) and V_i(x, t), respectively, obey

$\begin{array}{l}{\widehat{L}}_e{V}_e(x,t)=a_e{\displaystyle {\int }_{\mathcal{D}}{K}_e}(x-y)S_e[V_e(y,t)-V_i(y,t)]dy+I(x,t)\\ \begin{array}{cc}{\widehat{L}}_i{V}_i(x,t)=a_i{\displaystyle {\int }_{\mathcal{D}}{K}_i}(x-y)S_i[V_e(y,t)-V_i(y,t)]dy& (1)\end{array}\end{array}$

with the circular spatial population domain $\mathcal{D}$ of length L, i.e., assuming periodic boundary conditions. The model under study differs from some other previous models, e.g., by Liley and Bojak (2005), by the implementation of synaptic action, generation of action potentials or axonal connectivity [see also the work of Coombes et al. (2007) for a comparison of the current model and other models]. The functionals S_e[·] and S_i[·] are continuously increasing and represent the population firing rate of excitatory and inhibitory neurons, respectively. In the population the single neurons are connected by a complex system of axons from neuron somata to synapses. The kernels K_e(x) and K_i(x) are the probability density of such connections in the population. Here axonal transmission delay is neglected for simplicity although it is straightforward to include it in this type of model (Hutt and Longtin, 2009). In addition, the model considers excitatory and inhibitory synapses, ${\widehat{L}}_e={\widehat{L}}_e(d/dt)$ and ${\widehat{L}}_i={\widehat{L}}_i(d/dt)$ denote functional operators describing the corresponding temporal synaptic response phase and the factors a_e, a_i represent the corresponding synaptic efficacies.

Mathematically, the differential operators are the inverse of the integral operators in

$\begin{array}{cc}V(t)={\displaystyle {\int }_{-\infty }^{t}h}(t-\tau )P(\tau ) d\tau & (2)\end{array}$

where h(t) is the synaptic response function, or more precisely the electric current response in the synaptic receptor to an impact of binding neurotransmitters (Koch, 1999). The function P(τ) > 0 is the mean pulse activity arriving at the synapses. In a reasonable approximation, the response function reads

$h(t)=\frac{a}{\tau }{e}^{-t/\tau }$

with the decay time constant τ and the synaptic efficacy a > 0. Then the response amplitude is h(0) = a/τ and the charge transferred in the receptor ρ = a is the time integral over the current flow. The corresponding differential operator stipulates $\widehat{L}V(t)=aP(t)$ leading to

$\widehat{L}\left(\frac{\partial }{\partial t}\right)=\tau \frac{\partial }{\partial t}+1$

and Equation (2) re-casts to

$\tau \frac{\partial V(t)}{\partial t}+V(t)=aP(t).$

These expressions hold for excitatory and inhibitory synapses.

The synaptic receptors are major targets of anesthetic agents. The present work considers the action of propofol on inhibitory synaptic and extra-synaptic GABA_A-receptors. The former receptor is supposed to be a major anesthetic target (Franks and Lieb, 1994) and there is growing evidence that extra-synaptic inhibitory receptors may play an important role in anesthesia as well (Orser, 2006; Hutt, 2012). The subsequent sections consider effects on synaptic receptors due to the well-established experimental evidence. Hence, the synaptic parameters of inhibitory synaptic receptors depend on the anesthetic concentration and are parameterized by the factor p ≥ 1 (Steyn-Ross et al., 2001), i.e., the decay time of inhibitory synapses τ₂ = τ₂(p) and the corresponding synaptic efficacy a_i = a_i(p) depend on p.

The input in Equation (1) fluctuates randomly in space and time with ξ(x, t) about a constant value I₀ = const, i.e., I(x, t) = I₀ + ξ(x, t). The random fluctuations are independent in space and time and thus obey 〈ξ(x, t)〉 = 0, 〈ξ(x, t)ξ(y, T)〉 = 2Dδ(t − T)δ(x − y), where 〈·〉 denotes the ensemble average.

Considering the latter definitions of synaptic properties, anesthetic action and external input, the final model equations read

$\begin{array}{cc}\begin{array}{l}{\tau }_1\frac{\partial V_e(x,t)}{\partial t}=-V_e(x,t)+a_e{\displaystyle {\int }_{\mathcal{D}}{K}_e}(x-y)S_e[V_e(y,t)\\ -V_i(y,t)]dy+I_0+\xi (x,t)\\ {\tau }_2(p)\frac{\partial V_i(x,t)}{\partial t}=-V_i(x,t)+a_i(p){\displaystyle {\int }_{\mathcal{D}}{K}_i}(x-y)S_i[V_e(y,t)\\ -V_i(y,t)]dy\end{array}& (3)\end{array}$

with the decay time of excitatory synapses τ₁.

Assuming that the random fluctuations are small and do not affect the stationary state [in contrast to recent results gained from non-linear systems (Hutt et al., 2007; Hutt, 2008)], the stationary state V_e(x, t) = V⁰_e = const, V_i(x, t) = V⁰_i = const obeys V⁰_e = a_eS_e[V₋] + I₀, V⁰_i = a_i(p)S_i[V₋] with V₋ = V⁰_e − V⁰_i = a_eS_e[V₋] − a_i(p)S_i[V₋] + I₀ (Hutt and Longtin, 2009).

2.2. Theoretical Power Spectrum

To compute the power spectrum, we employ the method of Greens function. Let us assume the activity variable vector x(t) ∈ ℛ^N, the matrix A, the external input vector ξ(t) ∈ ℛ^N, the Greens function matrix G(t) ∈ ℛ^{N × N} and

$\dot{x}(t)=Ax+\xi (t).$

Then, for t → ∞, the solution of the system is

$\begin{array}{cc}x(t)={\displaystyle {\int }_{-\infty }^{\infty }G}(t-\tau )\xi (\tau )d\tau .& (4)\end{array}$

and the Greens function obeys

$\dot{G}-AG(t)=1\delta (t)$

with te unitary matrix 1 ∈ ℛ^{N × N}. Applying the Fourier transform

$\begin{array}{cc}G(t)=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\tilde{G}}(\omega )e^{i\omega t}d\omega .& (5)\end{array}$

yields

$\begin{array}{cc}\tilde{G}(\omega )=\frac{1}{\sqrt{2\pi }}{(i\omega 1-A)}^{-1}& (6)\\ =\frac{1}{\sqrt{2\pi }}\frac{F(i\omega )}{P(i\omega )}& (7)\end{array}$

with the matrix F(iω) and the characteristic polynom P(iω). The matrix F in Equation (7) includes the matrix elements of the inverse of iω1 − A and the characteristic polynom P represents the corresponding matrix determinant.

Inserting Equation (6) into (5) allows to compute G(t) by the residue theorem in functional analysis

$\begin{array}{cc}G(t)=2\pi i{\displaystyle \sum _{n=1}^{r}Res}(z_n,t)\Theta (t)& (8)\end{array}$

with the Heaviside function Θ(t) and the residues matrix Res(z_n, t) of F(z)/P(z) at the roots z_n of the characteristic equation P(z) = 0. The condition t > 0 considered by the Heaviside function is the mathematical condition for the validity of Equation (8) while it also guarantees the causality of the system response. Equation (8) together with Equation (4) determines the time dependence of the solution x and is computed explicitly in section 3.4.

Finally, the power spectral density matrix S(ω) of x is the Fourier transform of the auto-correlation function matrix 〈x^t(t)x(t − T)〉 (Wiener-Khinchine Theorem) leading to

$S(\omega )=2D\sqrt{2\pi }\tilde{G}(\omega ){\tilde{G}}^t(-\omega ),$

where the high index t denotes the transposed vector or matrix.

3. Results

3.1. Effect of Propofol

The effect of the anesthetic propofol on neural properties is manifold (Alkire et al., 2008). It affects properties of membrane ion channels, synaptic receptors and extra-synaptic receptors, see Franks and Lieb (1994) for a review. Kitamura et al. (2002) have revealed in an experimental study how propofol affects post-synaptic phasic responses of inhibitory synapses to spontaneous neurotransmitter release. They have found that the decay time constant decreases with increasing anesthetic blood concentration, the charge transfer increases while the amplitude of the responses remains constant. Some previous studies (Hutt and Longtin, 2009; Hindriks and van Putten, 2012) have implemented these effects for a bi-exponential synaptic response function. The present work considers an exponential decay phase due to its mathematical simplicity, which nevertheless reflects the major anesthetic impact. To this end, the phasic response at inhibitory synapses is determined by the decay time constant and the response amplitude. Introducing the parameter p = τ₂(p)/τ₂(1), p = 1 reflects the absence of anesthetic agents and increases of the anesthetic concentration yields an increase of p. In addition Kitamura et al. (2002) have shown that the amplitude in cortical neurons remains constant or changes slightly only, i.e., h(0) = a/τ₂(p)≈ const and the charge transfer ρ(p) = a increase with increasing anesthetic concentration. In the case of constant the model

$\begin{array}{cc}{\tau }_2(p)={\tau }_2(1)p,a_i(p)=H_{0}p.& (9)\end{array}$

with the constant H₀ > 0. These choices of anesthetic actions reflect different synaptic mechanism. The first relation reflects the increase of the time constant with increasing anesthetic concentration, and the second one both the constant amplitude and the resulting increasing charge transfer. However, synapses may have different properties in different brain areas, e.g., synapses in cortico-thalamic connections increase their amplitude with increasing propofol concentrations (Ying and Goldstein, 2005). Hence the relations (Equation 9) are specific properties of inhibitory synapses on cortical neurons only.

3.2. The Linear Model

The following investigation considers the stability of stationary states and spectral properties of small deviations about them. These small deviations represent fluctuating currents on the dendritic trees of the neurons in the population generating an electric field on the scalp. They originate from fluctuations in the neuron membranes or from spontaneous neurotransmitter emission at synapses. The generated electric field is measured in terms of voltage differences between two spatial locations on the scalp which is the electroencephalogram (EEG) (Nunez and Srinivasan, 2006).

The small fluctuations about the stationary states u_e(x, t) = V_e(x, t) − V⁰_e, u_i(x, t) = V_i(x, t) − V⁰_i obey

$\begin{array}{l}{\tau }_1\frac{d{u}_e(x,t)}{dt}=-u_e(x,t)+a_e{S}_e^{\prime }{\displaystyle {\int }_{\mathcal{D}}{K}_e}(x-y)(u_e(y,t)\\ -u_i(y,t))dy+\xi (x,t)\\ {\tau }_2(p)\frac{d{u}_i(x,t)}{dt}=-u_i(x,t)+a_i(p)S_i^{\prime }{\displaystyle {\int }_{\mathcal{D}}{K}_i}(x-y)(u_e(y,t)\\ -u_i(y,t))dy\end{array}$

with the somatic non-linear gain S′_{e, i} = dS_{e, i}(x)/dx at x = {V⁰_e, V⁰_i}. By virtue of the finite spatial domain, the small deviations about the stationary state may be expanded into a discrete infinite Fourier series. The major contribution of neuronal activity to encephalographic activity on the scalp is modeled successfully by a spatially constant Fourier mode u_e(x, t) = x(t), u_i(x, t) = y(t) (Nunez and Srinivasan, 2006). Then x(t), y(t) obey

$\begin{array}{cc}\begin{array}{l}{\tau }_1\frac{dx(t)}{dt}=(-1+N_1)x(t)-N_{1}y(t)+\gamma (t)\\ {\tau }_2(p)\frac{dy(t)}{dt}=N_{2}x(t)+(-1-N_2)y(t)\end{array}& (10)\end{array}$

with the synaptic non-linear gains $N_1=a_e{S}_e^{\prime }{\tilde{K}}_e(0)\sqrt{L}, N_2=N_2(p)=a_i(p)S_i^{\prime }{\tilde{K}}_i(0)\sqrt{L}$, the spatial Fourier transform of the kernels ${\tilde{K}}_e(k), {\tilde{K}}_i(k)$ and the spatial Fourier transform of the external noise at zero wavenumber γ(t). We point out that 〈γ(t)〉 = 0, 〈γ(t)γ(T)〉 = 2Dδ(t − T).

3.3. Stability Analysis

To study the dynamics about the stationary state, at first let us neglect the external input ξ(t) since in a first approximation the stability of the linear system does not depend on the external input. Then the stationary state is asymptotically stable if the characteristic equation of Equation (10)

$\begin{array}{cc}{\lambda }^2-\lambda Tr+det=0& (11)\end{array}$

with λ ∈ ℂ and

$\begin{array}{cc}Tr=\frac{N_1-1}{{\tau }_1}-\frac{N_2+1}{{\tau }_2}& (12)\end{array}$

$\begin{array}{cc}det=\frac{N_1{N}_2}{{\tau }_1{\tau }_2}-\frac{N_1-1}{{\tau }_1}\frac{N_2+1}{{\tau }_2}.& (13)\end{array}$

has solutions Re(λ) < 0. Here Tr and det are the trace and determinant of the linear matrix in Equation (10), respectively. In addition, the stationary state is a stable focus if Im(λ) = Ω ≠ 0, Ω ∈ ℝ and

${\Omega }^2=\frac{N_1{N}_2}{{\tau }_1{\tau }_2}-\frac{1}{4}{\left(\frac{N_1-1}{{\tau }_1}+\frac{N_2+1}{{\tau }_2}\right)}^2.$

We observe immediately from Equation (12), that the stable focus is asymptotically stable if Tr < 0 or N₁ < 1 for all other parameters and the system may lose stability only if N₁ > 1. Since the present work aims to give conditions for certain oscillation frequencies in the population, the subsequent part of the work considers stable foci only. In addition the inhibitory synaptic time scale τ₂ = τ₂(p) depends on the anesthetic concentration, but not τ₁. Thus τ₁ is treated as a constant.

The new variables a = (N₁ − 1)², b = N₁N₂ − N₁ + N₂ + 1 and c = (N₂ + 1)² depend solely on N₁, N₂, and simplify the notation in the following analysis. Then for N₁ > 1, stable foci stipulate

$\begin{array}{cc}Tr=\frac{\sqrt{a}}{{\tau }_1}-\frac{\sqrt{c}}{{\tau }_2}<0\text{ or }\frac{{\tau }_2}{{\tau }_1}<\sqrt{\frac{c}{a}},& (14)\end{array}$

$\begin{array}{cc}\Omega =\frac{1}{2}\sqrt{\frac{1}{{\tau }_1}\left(\frac{2b}{{\tau }_2}-\frac{a}{{\tau }_1}-\frac{c{\tau }_1}{{\tau }_2^2}\right)}\in ℝ.& (15)\end{array}$

The last equation implies that the determinant is positive definite, i.e.,

$d({\tau }_2)=a{\tau }_2^2-2b{\tau }_1{\tau }_2+c{\tau }_1^2<0.$

A comparison to Equation (15) reveals that d(τ₂) = 0 leads to Ω = 0 which permits to compute the range of τ₂ for which the system exhibits stable foci:

$\begin{array}{cc}\begin{array}{l}\frac{{\tau }_2^{-}}{{\tau }_1}\le \frac{{\tau }_2}{{\tau }_1}\le \frac{{\tau }_2^{+}}{{\tau }_1}\\ \frac{{\tau }_2^{-}}{{\tau }_1}=\frac{b}{a}\left(1+\sqrt{1-\frac{ac}{b^2}}\right), \frac{{\tau }_2^{+}}{{\tau }_1}=\frac{b}{a}\left(1-\sqrt{1-\frac{ac}{b^2}}\right)\end{array}& (16)\end{array}$

implying (b/a)² > c/a. Condition (16) constrains the relation of both synaptic time scales τ₂/τ₁ by N₁ and N₂, see Figure 1 for the corresponding parameter space.

FIGURE 1

Figure 1. Illustration of the parameter areas which exhibit stable, unstable, and no oscillations. The upper and lower solid lines denote τ⁻₂/τ₁ and τ⁺₂/τ₁, respectively. The dotted line denotes the maximum frequency f_max = 2πω_m with τ₂/τ₁ = b/c, cf. Equation (20), and the stability threshold is given by ${\tau }_2/{\tau }_1=\sqrt{c/a}$, see Equation (14).

To learn more about the dynamics of the model, we consider parameters yielding strong oscillations with a predefined frequency, such as f = 4 Hz (δ-band), f = 10 Hz (α-band), or f = 15 Hz (β-band) as observed in experiments (Cimenser et al., 2011). To this end, we fix the frequency Ω = 2πf in Equation (15). Then inserting the threshold condition ${\tau }_2/{\tau }_1=\sqrt{c/a}$ given in Equation (14) into (15) for a fixed τ₁ yields a relation between N₁ and N₂, see Figure 2A. Similarly the condition Ω = 0 given by Equation (16) determines the values of N₁ and N₂ at the threshold of oscillations. Figures 2B–E plot the parameter space τ₂/τ₁ − N₂ where the system exhibits stable oscillations. The figure also shows how the values change when increasing p (arrows in panels), i.e., increasing the propofol concentration, and it turns out that the system always approaches the stability threshold.

FIGURE 2

Figure 2. Conditions for Hopf-instabilities and the effect of propofol action. (A) threshold of Hopf-instability at 4, 10, and 15 Hz dependent on N₁, N₂. Panels (B,C) show parameters for Hopf-instabilities at 4 Hz dependent on N₂ for which the corresponding values of N₁ obey N₁ < N_1,min = 1.12 and (C) N₁ ≥ 1.12, respectively. Panels (D,E) show parameters for Hopf-instabilities at 10 Hz for N₁ < N_1,min = 1.45 and N₁ ≥ 1.45, respectively. In each panel (B–E) the top line reflects no oscillations, i.e., f = 0 Hz, and the bottom line denotes the corresponding stability threshold. Hence the shaded areas include parameters for which the system exhibits stable oscillations. The black solid and dotted lines reflect the situation of no anesthetic action, i.e., p = 1.0. The red-dashed line denotes the values of τ₁/τ₂ × p with p = 1.3. The arrows illustrate how the parameters change when anesthetic action is increased, i.e., when p is increased.

3.4. The Spectral Power

To compute the power spectrum, we employ the method of Greens function. Equations (10) show that the solutions of the system obey

$x(t)={\displaystyle {\int }_{-\infty }^{\infty }\text{}}{G}_{11}(t-\tau )\gamma (\tau )d\tau ,y(t)={\displaystyle {\int }_{-\infty }^{\infty }{G}_{21}}(t-\tau )\gamma (\tau )d\tau$

with the Greens functions

$\begin{array}{cc}{G}_{11}(t)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }\frac{i\omega +(N_2+1)/{\tau }_2}{P(i\omega )}}{e}^{i\omega t}d\omega & (17)\end{array}$

$\begin{array}{cc}{G}_{21}(t)=-\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }\frac{N_2/{\tau }_2}{P(i\omega )}}{e}^{i\omega t}d\omega & (18)\end{array}$

and P(iω) = −ω² − iωTr + det using the definitions of Tr and det in Equations (12) and (13). In fact, P(iω) in Equations (17) and (18) is the characteristic polynom in Equation (11) for λ = iω.

In the following in a good approximation (Nunez and Srinivasan, 2006) we assume that the experimental encephalographic data that we want to model originates from the excitatory synapses, i.e., x(t). This assumption is reasonable since the encephalographic activity is observed due to aligned apical dendritic branches and more excitatory synapses than inhibitory synapses are located on the apical branches of dendrites on cortical neurons. However, it is also possible to derive the power spectral density for the difference of excitatory and inhibitory potentials x(t) − y(t) (Hutt and Longtin, 2009). Now the application of the residue theorem allows to compute the integrals in Equation (17) by

$G_{11}(t)=\frac{1}{2\pi i}{\displaystyle \sum _{n=1}^2\text{Res}}(z_n)\Theta (t)$

with the residues

$\begin{array}{l}\text{Res}(z_1)=\frac{1}{2\pi i}\frac{{\lambda }_1+(N_2+1)/{\tau }_2}{{\lambda }_1-{\lambda }_2}{e}^{{\lambda }_{1}t},\\ \text{Res}(z_2)=\frac{1}{2\pi i}\frac{{\lambda }_2+(N_2+1)/{\tau }_2}{{\lambda }_2-{\lambda }_1}{e}^{{\lambda }_{2}t}\end{array}$

and the roots λ₁ = λ*₂ = R + iΩ of the characteristic polynom P(λ) = 0 in Equation (11). This solution is valid if and only if Res(λ_n) = R < 0. These yields

$G_{11}(t)=e^{Rt}\left(\frac{R+Z}{\Omega }\sin (\Omega t)+\cos (\Omega t)\right)\Theta (t)$

with Z = −(N₂ + 1)/τ₂, i.e., the Greens function of x, and hence x(t) itself oscillates with frequency Ω and is damped with the factor |R|. In other words, the solution x oscillates with the imaginary part of the root of the characteristic equation, and this frequency is already determined in the stability analysis.

Finally, the power spectral density S(ω) of x is the Fourier transform of the auto-correlation function 〈x(t)x(t − T)〉:

$\begin{array}{cc}\begin{array}{l}S(\omega )=2D\sqrt{2\pi }|{\tilde{G}}_{11}(\omega ){|}^2\\ =2D\sqrt{2\pi }\frac{Z^2+{\omega }^2}{{(R^2+{\Omega }^2-{\omega }^2)}^2+4R^2{\omega }^2}\end{array}& (19)\end{array}$

where ${\tilde{G}}_{11}(\omega )$ is the Fourier transform of the Greens function G₁₁(t).

3.5. The Frequency of Maximum Power Spectral Density

According to the reversed-engineering approach motivated in the previous section, this section aims to derive further conditions on model constants for certain oscillations close to instabilities. This vicinity to the stability threshold guarantees a small damping factor R and thus the power peak is located close to Ω.

There is a maximum frequency ω_m of the oscillation frequency with respect to τ₂ given in Equation (15) and reached at

$\begin{array}{cc}\frac{{\tau }_2}{{\tau }_1}=\frac{b}{c} \to {\omega }_m=\frac{1}{2}\sqrt{c/{\tau }_2^2-a/{\tau }_1^2}\text{ if }\frac{{\tau }_2}{{\tau }_1}<\sqrt{c/a}.& (20)\end{array}$

The last condition is identical to the stability condition (14), i.e., stable systems always have a non-vanishing maximum frequency ω_m > 0. Together with Equation (16):

• If $b/c<\sqrt{c/a}$, then the frequency Ω may increase or decrease while increasing τ₂ with

$\begin{array}{cc}\frac{d\Omega }{d{\tau }_2}>0\text{ for }\frac{{\tau }_2^{-}}{{\tau }_1}\le \frac{{\tau }_2}{{\tau }_2}<\frac{b}{c} ,& (21)\end{array}$

$\begin{array}{cc}\frac{d\Omega }{d{\tau }_2}\le 0\text{ for }\frac{b}{c}\le \frac{{\tau }_2}{{\tau }_2}\le \sqrt{\frac{c}{a}} .& (22)\end{array}$

• If $b/c\ge \sqrt{c/a}$, then the frequency Ω increases only while increasing τ₂ with

$\begin{array}{cc}\frac{d\Omega }{d{\tau }_2}>0\text{ for }\frac{{\tau }_2^{-}}{{\tau }_1}\le \frac{{\tau }_2}{{\tau }_2}\le \sqrt{\frac{c}{a}} .& (23)\end{array}$

Figure 1 shows the case $b/c<\sqrt{c/a}$. There increasing τ₂ by increasing p from small frequencies on the left border for constant τ₁ increases the oscillation frequency until reaching the dotted line, i.e., dΩ/dτ₂ > 0. Then a further increase of p decreases the frequency of the oscillations again, dΩ/dτ₂ < 0. Although this reasoning assumes that N₂(p) does not change with p, it gives a first insight into the dependence of the systems oscillation frequencies on p.

Now considering the power spectral density (Equation 19), for R⁴ + R²(2Ω² − 3Z²) + Ω²(Ω² + Z²) > 0 S(ω) has a global maximum at

${\Omega }_{\text{peak}}=\sqrt{{\Omega }^2+R^2\left(1-4\frac{Z^2}{R^2+{\Omega }^2+Z^2}\right)}.$

We learn that a vanishing real part of the characteristic root |R| yields a global maximum of the power spectral density at the imaginary part of the characteristic root Ω. From a reversed-engineer point of view, one could say that a strong peak in the experimental power spectral density reflects a characteristic root in the underlying linear system with a small real part. This way to interpret the spectral results allows to find analytical conditions for physiological parameters, as will be seen below.

If R⁴ + R²(2Ω² −3Z²) + Ω²(Ω² + Z²) < 0, then the system exhibits a global maximum at Ω_peak = 0. Since R and Ω depend on the anesthetic concentration, i.e., the parameter p, it is interesting to examine how Ω_peak depends on R and Ω:

$\begin{array}{cc}\frac{d{\Omega }_{\text{peak}}}{dR}>0\text{ for }{R}^2>Z^2(p)-{\Omega }^2 \text{or} R^2<Z^2(p)-{\Omega }^2& (24)\end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{d{\Omega }_{\text{peak}}}{d\Omega }>0\text{ for }{R}^2>\sqrt{2}|Z(p)|\Omega -{\Omega }^2+Z^2(p) \\ \text{ or }{R}^2<-\sqrt{2}|Z(p)|\Omega -{\Omega }^2+Z^2(p)\end{array}& (25)\end{array}$

Figure 3 illustrates these conditions and shows that Ω_peak is increased or decreased subjected to values of R and Ω. Importantly, we observe in Figure 3 that Ω_peak increases with R and Ω for large frequencies Ω.

FIGURE 3

Figure 3. Illustration of the conditions (24) and (25). The regions are marked as follows: (+): dΩ_peak/dΩ > 0, dΩ_peak/dR > 0; (∓): dΩ_peak/dΩ < 0, dΩ_peak/dR > 0; (−): dΩ_peak/dΩ < 0, dΩ_peak/dR < 0; (±): dΩ_peak/dΩ > 0, dΩ_peak/dR < 0.

It remains to examine how R² and Ω² depend on p to finally gain the full description how Ω_peak changes with p. To this end, we consider the specific assumption (Equation 9) on the anesthetic action yielding the specific dependence of Ω² and R² to p

$\begin{array}{l}\frac{d^{2}R}{dp}=-2\frac{R}{{\tau }_2(1)}\frac{1}{p^2}>0\\ \frac{d{\Omega }^2}{dp}=-\frac{2}{{\tau }_2(1)p^2}\left(A-\frac{1}{{\tau }_2(1)p}\right)\end{array}$

with A = (1 − N₁)/τ₁ − N₂(1)/τ₂(1), N₂(p) = N₂(1)p, and τ₂ = τ₂(1)p. Defining p₀ = (N₂(1) + τ₂(1)(N₁ − 1)/τ₁)⁻¹ and re-calling p ≥ 1, then there exist two distinct cases:

• If p₀ < 1, then dΩ²/dp > 0 for all p. Specifically, this stipulates

$N_2(1)>1+\frac{{\tau }_2(1)}{{\tau }_1}-\frac{{\tau }_2(1)}{{\tau }_1}{N}_1.$

• If p₀ ≥ 1, then there is a small interval 1 ≤ p ≤ p₀ for which dΩ²/dp ≤ 0. For larger p, dΩ/dp > 0. For instance, for p₀ = 1.3, i.e.,

$N_2(1)=1/1.3+\frac{{\tau }_2(1)}{{\tau }_1}-\frac{{\tau }_2(1)}{{\tau }_1}{N}_1,$

the frequency decreases with increasing p for clinically reasonable concentrations, i.e., 1 ≤ p ≤ 1.3 (Hutt and Longtin, 2009).

Figure 4 illustrates the parameter space where dΩ²/dp has different signs and we observe that the system may exhibit either increasing or decreasing frequencies while increasing p.

FIGURE 4

Figure 4. Parameter space for different signs of dΩ²/dp. The red line denotes p₀ = 1.3 and the parameter space of clinical doses (red) implies dΩ²/dp < 0 for clinically reasonable drug doses, i.e., p ≤ 1.3. An additional parameter is r = τ₂(1)/τ₁ = 1.11.

To elucidate the systems behavior at the stability threshold, we set τ₂(1)/τ₁ = (N₂(1) + 1)/(N₁ − 1) according to Equation (14) and find p₀ = 1/(2N₂(1) + 1) < 1. Consequently, oscillations close to the stability threshold always increase their frequencies with increasing p. In contrast, for 1 < τ₂(1)/τ₁ « (N₂(1) + 1)/(N₁ − 1), i.e., systems far from the stability threshold, may exhibit values p₀ > 1.

Finally, Figure 5 illustrates the temporal dynamics in two different frequency bands close to the corresponding Hopf-instabilities and shows the effect of increased propofol concentration. The power increases and the peak of maximum power moves to larger frequencies as predicted by the theory and as observed in experiments (Hindriks and van Putten, 2012).

FIGURE 5

Figure 5. Simulated time series of x(t) + V⁰_e and the corresponding power spectrum of x(t). (A) N₁ = 1.1, N₂ = p × 0.25128, τ₂ = τ₁ × p/0.10 generating a maximum power in the δ-band. (B) N₁ = 1.1, N₂ = 0.2236 × p, τ₂ = τ₁ × p/0.10 generating a maximum in the α-band. Other parameters are τ₁ = 2ms, the noise strength κ = 0.01 mV and V⁰_e = −60 mV. The Equation (10) is simulated with a Euler–Maryuama method for 200 s, time step was Δt = 0.05 ms.

4. Discussion

The introduced model considers first-order synaptic responses and take into account experimental findings on the propofol effect in synaptic GABA_A-receptors in cortical neurons. It describes the evolution of neural populations on a mesoscopic level involving major properties of underlying neurons and synapses on the microscopic description level.

The analytical study reveals that the frequency of maximum power may increase or decrease with increasing anesthetic concentration subjected to the physiological constants, cf. Figure 3. In detail, close to the oscillatory instability the frequency of maximum power always increases with increasing p as observed in EEG (Gugino et al., 2001; Feshchenko et al., 2004; Hindriks and van Putten, 2012), whereas far from the stability threshold the maximum power frequency may also decrease as observed recently in EEG (Ching et al., 2010; Cimenser et al., 2011), cf. Figures 3, 4. The analytical treatment shows clearly that these two findings depend strongly on the physiological parameters, which are derived analytically in section 3.5, i.e., the phenomena depend on the brain area in which they are generated.

Moreover, the analytically predicted increase of the power at higher frequencies explains the power enhancement in the α- and β-band in anesthesia (Gugino et al., 2001; McCarthy et al., 2008) by a dynamic oscillatory instability. In fact, the analytical treatment in the present work suggests that power enhancement always starts from oscillatory activity at lower frequencies and are generated at slightly higher frequencies with increased concentration (Figure 5). This is in accordance to previous EEG-studies (Gugino et al., 2001; Hindriks and van Putten, 2012) showing power enhancement induced in the α- and β-band.

The analytical discussion in section 3.3 also predicts that decreasing the inhibitory time constant always moves the system toward an oscillatory instability, cf. Figure 1, and hence increases the spectral power, whereas increasing the charge transfer yields a stabilization of the system due to dR/dP < 0 and consequently a decrease of power. Hence the balance between decay prolongation and increased charge transfer decides on the change of the spectral power and the shift of the frequency peaks. Since many anesthetics share this balance in the major target GABA_A-receptor (Alkire et al., 2008) and exhibit similar EEG-change (Gugino et al., 2001; Kuizenga et al., 2001), the presented work suggests that this balance reflects one of the major underlying mechanisms during the sedation phase in GA.

Previous studies (Bojak and Liley, 2005; Hindriks and van Putten, 2012) already have explained the power enhancement in anesthesia by an oscillatory instability in high-dimensional neural models. As one of the first, the present work gives analytical conditions on physiological parameters for this effect, while Bojak and Liley (2005) and Hindriks and van Putten (2012) mainly performed numerical studies. The recent work of Hindriks and van Putten (2012) resembles in some aspects the analytical approach of the present work by discussing in some detail the dynamics of superimposed oscillation modes subjected to the propofol concentration. However, no analytical conditions are given due to the higher model complexity.

It is important to point out that the current model is low-dimensional, physiologically reasonable and analytically treatable but still able to explain the neural phenomenon of the frequency shift to larger values. Bojak and Liley (2005) and Hindriks and van Putten (2012) have not performed a detailed analytical study of this phenomenon and have not derived analytical conditions under which it may occur. The current work shows that already a rather simple coupling of excitation and inhibition in cortical neural networks is sufficient to explain this phenomenon. This is concluded partially by Hindriks and van Putten (2012) based on a small numerical study, whereas the present work shows this explicitly. However, Hindriks and van Putten (2012) also argue that the cortico-thalamic feedback should be negative to gain this effect. The presented model does not need the thalamic feedback loop for the explanation.

In principle, the present work extends the work of Bojak and Liley (2005) studying just numerically a rather complicated model with tens of unknown parameters, while the current model allows to achieve some insights into the effect of few parameters. For instance, the maximum of spectral power and the corresponding frequency are highly sensitive to modification of the relation between excitatory and inhibitory synaptic time scales. This confirms the general observation that many different anesthetics share effects on the activity spectrum.

The work both supports the hypothesis of cortical generation of α-activity and predicts the presence of oscillating neural circuits where each circuit generates a certain oscillation observed experimentally. This can be observed in Figure 5 showing emerging δ- and α-activity for different parameters. Hence two neural circuits with different properties may explain the occurrence of both δ- and α-activity observed in experimental data. This interpretation of the results complements the findings of Hindriks and van Putten (2012) showing implicitly a linear decomposition into eigenmodes with corresponding eigenvalues and manifests the notion of interacting oscillation modes generated by interacting networks as observed experimentally (Fries, 2009; Spaak et al., 2012).

Of course the present model is limited since it cannot explain the increase of activity in the δ-band as observed in experiments and modeled by the previous studies. In the present work, we observe clearly that this synchronous modeling of two rhythms is not possible since the model is too low-dimensional and just can describe a single rhythm such as the α- or the δ-rhythm. Consequently further neural elements should be considered to gain this additional rhythm such as the thalamic loop. This may be possible due to the linear superposition of oscillatory activity originating from different networks. This linear superposition proposes the interaction of different sub-networks each oscillating in a certain frequency band. Future work will consider such entangled neural networks on the basis of the presented neural field to explain the spontaneous emergence or diminution of spectral peaks in experimentally observed data such as the δ-rhythms or transient phenomena such as paradoxical excitation.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The author thanks Meysam Hashemi for fruitful discussions and acknowledges funding from the European Research Council for support under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement No. 257253. The author thanks the referees for valuable comments.

References

Alkire, M., Hudetz, A., and Tononi, G. (2008). Consciousness and anesthesia. Science 322, 876–880.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Amari, S. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biol. Cybern. 27, 77–87.

Pubmed Abstract | Pubmed Full Text

Angelucci, A., and Bressloff, P. (2006). The contribution of feedforward, lateral and feedback connections to the classical receptive field center and extra-classical receptive field surround of primate v1 neurons. Prog. Brain Res. 154, 93–121.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bojak, I., and Liley, D. (2005). Modeling the effects of anesthesia on the electroencephalogram. Phys. Rev. E 71:041902. doi: 10.1103/PhysRevE.71.041902

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Boly, M., Moran, R., Murphy, M., Boveroux, P., Bruno, M. A., Noirhomme, Q., et al. (2012). Connectivity changes underlying spectral eeg changes during propofol-induced loss of consciousness. J. Neurosci. 32, 7082–7090.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ching, S., Cimenser, A., Purdon, P. L., Brown, E. N., and Kopell, N. J. (2010). Thalamocortical model for a propofol-induced-rhythm associated with loss of consciousness. Proc. Natl. Acad. Sci. U.S.A. 107, 22665–22670.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cimenser, A., Purdon, P. L., Pierce, E. T., Walsh, J. L., Salazar-Gomez, A. F., Harrell, P. G., et al. (2011). Tracking brain states under general anesthesia by using global coherence analysis. Proc. Natl. Acad. Sci. U.S.A. 108, 8832–8837.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Coombes, S., Venkov, N., Shiau, L., Bojak, I., Liley, D., and Laing, C. (2007). Modeling electrocortical activity through improved local approximations of integral neural field equations. Phys. Rev. E 76, 051901–051908.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ermentrout, G. (1998). Neural networks as spatio-temporal pattern-forming systems. Rep. Progr. Phys. 61, 353–430.

Ermentrout, G., and Cowan, J. (1979). A mathematical theory of visual hallucination patterns. Biol. Cybern. 34, 137–150.

Pubmed Abstract | Pubmed Full Text

Feshchenko, V. A., Veselis, R. A., and Reinsel, R. A. (2004). Propofol-induced alpha rhythm. Neuropsychobiology 50, 257–266.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Franks, N., and Lieb, W. (1994). Molecular and cellular mechanisms of general anesthesia. Nature 367, 607–614.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fries, P. (2009). Neuronal gamma-band synchronization as a fundamental process in cortical computation. Annu. Rev. Neurosci. 32, 209–224.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gugino, L. D., Chabot, R. J., Prichep, L. S., John, E. R., Formanek, V., and Aglio, L. S. (2001). Quantitative eeg changes associated with loss and return of conscious- ness in healthy adult volunteers anaesthetized with propofol or sevoflurane. Br. J. Anaesth. 87, 421–428.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hindriks, R., and van Putten, M. J. A. M. (2012). Meanfield modeling of propofol-induced changes in spontaneous eeg rhythms. Neuroimage 60, 2323–2344.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Huang, X., Troy, W., Schiff, S., Yang, Q., Ma, H., Laing, C., et al. (2004). Spiral waves in disinhibited mammalian neocortex. J. Neurosci. 24, 9897–9902.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hutt, A. (2008). Additive noise may change the stability of nonlinear systems. Europhys. Lett. 84:34003. doi: 10.1209/0295-5075/84/34003

CrossRef Full Text

Hutt, A. (ed.). (2011). Sleep and Anesthesia: Neural Correlates in Theory and Experiment. Springer Series in Computational Neuroscience, Vol. 15. New York, NY: Springer.

Hutt, A. (2012). The population firing rate in the presence of gabaergic tonic inhibition in single neurons and application to general anaesthesia. Cogn. Neurodyn. 6, 227–237.

Hutt, A., and Longtin, A. (2009). Effects of the anesthetic agent propofol on neural populations. Cogn. Neurodyn. 4, 37–59.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hutt, A., Longtin, A., and Schimansky-Geier, L. (2007). Additive global noise delays turing bifurcations. Phys. Rev. Lett. 98:230601. doi: 10.1103/PhysRevLett.98.230601

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kitamura, A., Marszalec, W., Yeh, J., and Narahashi, T. (2002). Effects of halothane and propofol on excitatory and inhibitory synaptic transmission in rat cortical neurons. J. Pharmacol. 304, 162–171.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Koch, C. (1999). Biophysics of Computation. Oxford: Oxford University Press.

Kuizenga, K., Wierda, J., and Kalkman, C. (2001). Biphasic eeg changes in relation to loss of consciousness during induction with thiopental, propofol, etomidate, midazolam or sevoflurane. Br. J. Anaesth. 86, 354–360.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liley, D., and Bojak, I. (2005). Understanding the transition to seizure by modeling the epileptiform activity of general anaesthetic agents. J. Clin. Neurophysiol. 22, 300–313.

Pubmed Abstract | Pubmed Full Text

McCarthy, M. M., Brown, E. N., and Kopell, N. (2008). Potential network mechanisms mediating electroencephalographic beta rhythm changes during propofol-induced paradoxical excitation. J. Neurosci. 28, 13488–13504.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Murphy, M., Bruno, M.-A., Riedner, B. A., Boveroux, P., Noirhomme, Q., Landsness, E. C., et al. (2011). Propofol anesthesia and sleep: a high-density eeg study. Sleep 34, 283–291.

Pubmed Abstract | Pubmed Full Text

Nunez, P., and Srinivasan, R. (2006). Electric Fields of the Brain: The Neurophysics of EEG. Oxford, NY: Oxford University Press.

Orser, B. (2006). Extrasynaptic gabaa receptors are critical targets for sedative-hypnotic drugs. J. Clin. Sleep Med. 2, S12–S18.

Pubmed Abstract | Pubmed Full Text

Schwabe, L., Obermayer, K., Angelucci, A., and Bressloff, P. (2006). The role of feedback in shaping the extra-classical receptive field of cortical neurons: a recurrent network model. J. Neurosci. 26, 9117–9126.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Spaak, E., Bonnefond, M., Maier, A., Leopold, D. A., and Jensen, O. (2012). Layer-specific entrainment of gamma-band neural activity by the alpha rhythm in monkey visual cortex. Curr. Biol. 22, 2313–2318.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Steyn-Ross, M., Steyn-Ross, D., and Sleigh, J. (2004). Modelling general anaesthesia as a first-order phase transition in the cortex. Prog. Biophys. Molec. Biol. 85, 369–385.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Steyn-Ross, M., Steyn-Ross, D., Sleigh, J., and Wilcocks, L. (2001). Toward a theory of the general-anesthetic-induced phase transition of the cerebral cortex: I. a thermodynamic analogy. Phys. Rev. E 64:011917. doi: 10.1103/PhysRevE.64.011917

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Steyn-Ross, M., Steyn-Ross, D., and Sleigh, J. W. (2012). Gap junctions modulate seizures in a mean-field model of general anesthesia for the cortex. Cogn. Neurodyn. 6, 215–225.

Steyn-Ross, M., Steyn-Ross, D., Sleigh, J. W., and Liley, D. T. J. (1999). Theoretical electroencephalogram stationary spectrum for a white-noise-driven cortex: evidence for a general anesthetic-induced phase transition. Phys. Rev. E 60, 7299–7311.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wilson, H., and Cowan, J. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ying, S., and Goldstein, P. (2005). Propofol suppresses synaptic responsiveness of somatosensory relay neurons to excitatory input by potentiating gabaa receptor chloride channels. Mol. Pain 1:2. doi: 10.1186/1744-8069-1-2

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text