Impact Factor 1.821

The Frontiers in Neuroscience journal series is the 1st most cited in Neurosciences

Original Research ARTICLE

Front. Comput. Neurosci., 05 February 2013 | https://doi.org/10.3389/fncom.2013.00002

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

  • INRIA CR Nancy - Grand Est, Team CORTEX, Villers-les-Nancy, France

The work introduces a linear neural population model that allows to derive analytically the power spectrum subjected to the concentration of the anesthetic propofol. The analytical study of the power spectrum of the systems activity gives conditions on how the frequency of maximum power in experimental electroencephalographic (EEG) changes dependent on the propofol concentration. In this context, we explain the anesthetic-induced power increase in neural activity by an oscillatory instability and derive conditions under which the power peak shifts to larger frequencies as observed experimentally in EEG. Moreover the work predicts that the power increase only occurs while the frequency of maximum power increases. Numerically simulations of the systems activity complement the analytical results.

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 Ve(x, t) and Vi(x, t), respectively, obey

L^eVe(x,t)=aeDKe(xy)Se[Ve(y,t)Vi(y,t)]dy+I(x,t)L^iVi(x,t)=aiDKi(xy)Si[Ve(y,t)Vi(y,t)]dy    (1)

with the circular spatial population domain 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 Se[·] and Si[·] 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 Ke(x) and Ki(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, L^e=L^e(d/dt) and L^i=L^i(d/dt) denote functional operators describing the corresponding temporal synaptic response phase and the factors ae, ai represent the corresponding synaptic efficacies.

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

V(t)=th(tτ)P(τ) dτ(2)

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)=aτet/τ

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 L^V(t)=aP(t) leading to

L^(t)=τt+1

and Equation (2) re-casts to

τV(t)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 GABAA-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 τ2 = τ2(p) and the corresponding synaptic efficacy ai = ai(p) depend on p.

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

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

     τ1Ve(x,t)t=Ve(x,t)+aeDKe(xy)Se[Ve(y,t)                          Vi(y,t)]dy+I0+ξ(x,t)τ2(p)Vi(x,t)t=Vi(x,t)+ai(p)DKi(xy)Si[Ve(y,t)                          Vi(y,t)]dy(3)

with the decay time of excitatory synapses τ1.

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 Ve(x, t) = V0e = const, Vi(x, t) = V0i = const obeys V0e = aeSe[V] + I0, V0i = ai(p)Si[V] with V = V0eV0i = aeSe[V] − ai(p)Si[V] + I0 (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

x˙(t)=Ax+ξ(t).

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

x(t)=G(tτ)ξ(τ)dτ.(4)

and the Greens function obeys

G˙AG(t)=1δ(t)

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

G(t)=12πG˜(ω)eiωtdω.(5)

yields

G˜(ω)=12π(iω1A)1(6) =12πF(iω)P(iω)(7)

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ω1A 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

G(t)=2πin=1rRes(zn,t)Θ(t)(8)

with the Heaviside function Θ(t) and the residues matrix Res(zn, t) of F(z)/P(z) at the roots zn 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 〈xt(t)x(tT)〉 (Wiener-Khinchine Theorem) leading to

S(ω)=2D2πG˜(ω)G˜t(ω),

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 = τ2(p)/τ2(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) = a2(p)≈ const and the charge transfer ρ(p) = a increase with increasing anesthetic concentration. In the case of constant the model

τ2(p)=τ2(1)p, ai(p)=H0p.(9)

with the constant H0 > 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 ue(x, t) = Ve(x, t) − V0e, ui(x, t) = Vi(x, t) − V0i obey

   τ1due(x,t)dt=ue(x,t)+aeSeDKe(xy)(ue(y,t)                         ui(y,t))dy+ξ(x,t)τ2(p)dui(x,t)dt=ui(x,t)+ai(p)SiDKi(xy)(ue(y,t)                          ui(y,t))dy

with the somatic non-linear gain Se, i = dSe, i(x)/dx at x = {V0e, V0i}. 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 ue(x, t) = x(t), ui(x, t) = y(t) (Nunez and Srinivasan, 2006). Then x(t), y(t) obey

     τ1dx(t)dt=(1+N1)x(t)N1y(t)+γ(t)τ2(p)dy(t)dt=N2x(t)+(1N2)y(t)(10)

with the synaptic non-linear gains N1=aeSeK˜e(0)L, N2=N2(p)=ai(p)SiK˜i(0)L, the spatial Fourier transform of the kernels K˜e(k), 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δ(tT).

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)

λ2λTr+det=0(11)

with λ ∈ ℂ and

Tr=N11τ1N2+1τ2(12)
det=N1N2τ1τ2N11τ1N2+1τ2.(13)

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

Ω2=N1N2τ1τ214(N11τ1+N2+1τ2)2.

We observe immediately from Equation (12), that the stable focus is asymptotically stable if Tr < 0 or N1 < 1 for all other parameters and the system may lose stability only if N1 > 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 τ2 = τ2(p) depends on the anesthetic concentration, but not τ1. Thus τ1 is treated as a constant.

The new variables a = (N1 − 1)2, b = N1N2N1 + N2 + 1 and c = (N2 + 1)2 depend solely on N1, N2, and simplify the notation in the following analysis. Then for N1 > 1, stable foci stipulate

Tr=aτ1cτ2<0 or τ2τ1<ca,(14)
Ω=121τ1(2bτ2aτ1cτ1τ22).(15)

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

d(τ2)=aτ222bτ1τ2+cτ12<0.

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

τ2τ1τ2τ1τ2+τ1τ2τ1=ba(1+1acb2), τ2+τ1=ba(11acb2)(16)

implying (b/a)2 > c/a. Condition (16) constrains the relation of both synaptic time scales τ21 by N1 and N2, see Figure 1 for the corresponding parameter space.

FIGURE 1
www.frontiersin.org

Figure 1. Illustration of the parameter areas which exhibit stable, unstable, and no oscillations. The upper and lower solid lines denote τ21 and τ+21, respectively. The dotted line denotes the maximum frequency fmax = 2πωm with τ21 = b/c, cf. Equation (20), and the stability threshold is given by τ2/τ1=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 τ2/τ1=c/a given in Equation (14) into (15) for a fixed τ1 yields a relation between N1 and N2, see Figure 2A. Similarly the condition Ω = 0 given by Equation (16) determines the values of N1 and N2 at the threshold of oscillations. Figures 2B–E plot the parameter space τ21N2 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
www.frontiersin.org

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 N1, N2. Panels (B,C) show parameters for Hopf-instabilities at 4 Hz dependent on N2 for which the corresponding values of N1 obey N1 < N1,min = 1.12 and (C) N1 ≥ 1.12, respectively. Panels (D,E) show parameters for Hopf-instabilities at 10 Hz for N1 < N1,min = 1.45 and N1 ≥ 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 τ12 × 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)=G11(tτ)γ(τ)dτ, y(t)=G21(tτ)γ(τ)dτ

with the Greens functions

G11(t)=12πiω+(N2+1)/τ2P(iω)eiωtdω(17)
G21(t)=12πN2/τ2P(iω)eiωtdω(18)

and P(iω) = −ω2iω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

G11(t)=12πin=12Res(zn)Θ(t)

with the residues

Res(z1)=12πiλ1+(N2+1)/τ2λ1λ2eλ1t,Res(z2)=12πiλ2+(N2+1)/τ2λ2λ1eλ2t

and the roots λ1 = λ*2 = 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

G11(t)=eRt(R+ZΩsin(Ωt)+cos(Ωt))Θ(t)

with Z = −(N2 + 1)/τ2, 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(tT)〉:

S(ω)=2D2π|G˜11(ω)|2        =2D2πZ2+ω2(R2+Ω2ω2)2+4R2ω2(19)

where G˜11(ω) is the Fourier transform of the Greens function G11(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 τ2 given in Equation (15) and reached at

τ2τ1=bc  ωm=12c/τ22a/τ12 if τ2τ1<c/a.(20)

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<c/a, then the frequency Ω may increase or decrease while increasing τ2 with

dΩdτ2>0 for τ2τ1τ2τ2<bc ,(21)
dΩdτ20 for bcτ2τ2ca .(22)

• If b/cc/a, then the frequency Ω increases only while increasing τ2 with

dΩdτ2>0 for τ2τ1τ2τ2ca .(23)

Figure 1 shows the case b/c<c/a. There increasing τ2 by increasing p from small frequencies on the left border for constant τ1 increases the oscillation frequency until reaching the dotted line, i.e., dΩ/dτ2 > 0. Then a further increase of p decreases the frequency of the oscillations again, dΩ/dτ2 < 0. Although this reasoning assumes that N2(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 R4 + R2(2Ω2 − 3Z2) + Ω22 + Z2) > 0 S(ω) has a global maximum at

Ωpeak=Ω2+R2(14Z2R2+Ω2+Z2).

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 R4 + R2(2Ω2 −3Z2) + Ω22 + Z2) < 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 Ω:

dΩpeakdR>0 for R2>Z2(p)Ω2 or R2<Z2(p)Ω2(24)
dΩpeakdΩ>0 for R2>2|Z(p)|ΩΩ2+Z2(p)                    or R2<2|Z(p)|ΩΩ2+Z2(p)(25)

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
www.frontiersin.org

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 R2 and Ω2 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 Ω2 and R2 to p

d2Rdp=2Rτ2(1)1p2>0dΩ2dp=2τ2(1)p2(A1τ2(1)p)

with A = (1 − N1)/τ1N2(1)/τ2(1), N2(p) = N2(1)p, and τ2 = τ2(1)p. Defining p0 = (N2(1) + τ2(1)(N1 − 1)/τ1)−1 and re-calling p ≥ 1, then there exist two distinct cases:

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

N2(1)>1+τ2(1)τ1τ2(1)τ1N1.

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

N2(1)=1/1.3+τ2(1)τ1τ2(1)τ1N1,

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Ω2/dp has different signs and we observe that the system may exhibit either increasing or decreasing frequencies while increasing p.

FIGURE 4
www.frontiersin.org

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

To elucidate the systems behavior at the stability threshold, we set τ2(1)/τ1 = (N2(1) + 1)/(N1 − 1) according to Equation (14) and find p0 = 1/(2N2(1) + 1) < 1. Consequently, oscillations close to the stability threshold always increase their frequencies with increasing p. In contrast, for 1 < τ2(1)/τ1 « (N2(1) + 1)/(N1 − 1), i.e., systems far from the stability threshold, may exhibit values p0 > 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
www.frontiersin.org

Figure 5. Simulated time series of x(t) + V0e and the corresponding power spectrum of x(t). (A) N1 = 1.1, N2 = p × 0.25128, τ2 = τ1 × p/0.10 generating a maximum power in the δ-band. (B) N1 = 1.1, N2 = 0.2236 × p, τ2 = τ1 × p/0.10 generating a maximum in the α-band. Other parameters are τ1 = 2ms, the noise strength κ = 0.01 mV and V0e = −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 GABAA-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 GABAA-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

Keywords: general anesthesia, propofol, neural fields, power spectrum, EEG

Citation: Hutt A (2013) The anesthetic propofol shifts the frequency of maximum spectral power in EEG during general anesthesia: analytical insights from a linear model. Front. Comput. Neurosci. 7:2. doi: 10.3389/fncom.2013.00002

Received: 12 December 2012; Accepted: 19 January 2013;
Published online: 05 February 2013.

Edited by:

Peter Beim Graben, Humboldt-Universität zu Berlin, Germany

Reviewed by:

Moira Steyn-Ross, University of Waikato, New Zealand
Rikkert Hindriks, Universitat Pompeu Fabra, Spain

Copyright © 2013 Hutt. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.

*Correspondence: Axel Hutt, INRIA CR Nancy - Grand Est, Team CORTEX, 615 rue du Jardin Botanique, 54602 Villers-les-Nancy, France. e-mail: axel.hutt@inria.fr