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

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.


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 longstanding 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, 2012Bojak and Liley, 2005;Hutt and Longtin, 2009;Hindriks and van Putten, 2012), or spikingneuron 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.

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 spatiotemporal 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, obeŷ 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 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, L e =L e (d/dt) andL i =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 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 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 stipulatesLV(t) = aP(t) leading toL and Equation (2) re-casts to 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 τ 2 = τ 2 (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 0 = const, i.e., I(x, t) = I 0 + ξ(x, t). The random fluctuations are independent in space and time and thus obey ξ( Considering the latter definitions of synaptic properties, anesthetic action and external input, the final model equations read 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 sta- (Hutt and Longtin, 2009).

THEORETICAL POWER SPECTRUM
To compute the power spectrum, we employ the method of Greens function. Let us assume the activity variable vector x(t) ∈ R N , the matrix A, the external input vector ξ(t) ∈ R N , the Greens function matrix G(t) ∈ R N×N anḋ Then, for t → ∞, the solution of the system is (4) and the Greens function obeyṡ with te unitary matrix 1 ∈ R N × N . Applying the Fourier transform yieldsG 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 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 where the high index t denotes the transposed vector or matrix.

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

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 ( 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 (Nunez and Srinivasan, 2006).
with the synaptic non-linear gains N 1 = a e S eK e (0) √ L, the spatial Fourier transform of the ker-nelsK e (k),K i (k) and the spatial Fourier transform of the external noise at zero wavenumber γ (t). We point out that

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) with λ ∈ C and 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, ∈ R and We observe immediately from Equation (12), that the stable focus is asymptotically stable if Tr < 0 or N 1 < 1 for all other parameters and the system may lose stability only if N 1 > 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.
A comparison to Equation (15) reveals that d(τ 2 ) = 0 leads to = 0 which permits to compute the range of τ 2 for which the system exhibits stable foci:  (20), and the stability threshold is given by τ 2 /τ 1 = √ c/a, see Equation (14).
implying (b/a) 2 > c/a. Condition (16) constrains the relation of both synaptic time scales τ 2 /τ 1 by N 1 and N 2 , see Figure 1 for the corresponding parameter space.
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 N 1 and N 2 , see Figure 2A. Similarly the condition = 0 given by Equation (16) determines the values of N 1 and N 2 at the threshold of oscillations. Figures 2B-E plot the parameter space τ 2 /τ 1 − N 2 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.

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 and P(iω) = −ω 2 − 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 with the residues Res(z 1 ) = 1 2πi Res(z 2 ) = 1 2πi 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 with Z = −(N 2 + 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(t − T) : whereG 11 (ω) is the Fourier transform of the Greens function G 11 (t).

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 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 • If b/c ≥ √ c/a, then the frequency increases only while increasing τ 2 with (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 N 2 (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 4 + R 2 (2 2 − 3Z 2 ) + 2 ( 2 + Z 2 ) > 0 S(ω) has a global maximum at 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 reversedengineer 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 4 + R 2 (2 2 − 3Z 2 ) + 2 ( 2 + Z 2 ) < 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 : 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 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.
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).

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 EEGstudies (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 ; 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 lowdimensional, 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.