In this section, we use a phase network with mean field coupling to get some insights about the novel network states observed in (*K*, *I*) the parameter space of the full network model. Coupling between each phase neuron via a mean field is formulated in section 2.4. Numerically we integrate Equation (6) to compute time averaged membrane potential, mean field as in Equation (10) (see section 2), global measure of coherence *R*θ as a function of *K*, *I* a parameter combination which is used in

the detailed network model. Time evolution of the above quantiles are shown in **Figure 7** for a parameterization *K* = 0.8, *I* = 0.001. The parameter choice is the same as the full network investigations. With this combination of parameters all the neurons synchronously spikes. Mean membrane potential-like quantity *V*mean(*t*) oscillates in phase with synchronized spike activity as plotted in **Figure 7B**. Here, *n* a quantity which determines the number of spikes per burst is kept at *n* = 1. Mean field also shows up and down states (Locked in time) and act as an oscillating drive to each individual neurons. The time series of the global order parameter *R*θ(*t*) for synchronized spiking is periodic in **Figure 7E**. Next, we show in **Figure 8** temporal evolution of the mean quantities for the choice of *K* = 0.8, *I* = 0.5. For medium values of mean field coupling strength and stimulus spread network shows PO behavior, where some of the neurons are firing incoherently and others are completely silent. This means for some parameterization network has two stable branches of solutions. It is important to note PO state of the network was observed in the full network for a comparable parameterization (see **Figure 6**). Time series for 10 neurons and their order parameter evolution in time is plotted in **Figure 8**. Three neurons are completely silent while other seven neurons are bursting with variable inter-burst intervals. As there is no noise in this system and coupling magnitude is set at high values as in the case of sync, this variability must be introduced by the heterogeneity in their individual firing rate threshold via stimulus spread.

**Figure 6** shows the parameter space diagram for the full and phase models presented in Equations (5) and (6–9). Phase boundaries are calculated by computing the mean field for both full and the phase model for different combination of (*K*,*I*) values on a two-dimensional grid. In the following subsection we would lay out the details for obtaining the phase transition boundaries semi analytically. Interestingly, Over a wide range of (*K* − *I*) values the collective dynamics of the two networks primarily show three distinct regions of interest which are close to each other in the parameter space. For sufficiently large K values holding *I* fixed to a narrow range of values near zero, the system converges to a state of partial oscillations in which the some of the neurons are not firing altogether, while the others display IN oscillations. Very large *K* values result in damping of oscillation activity and all the neurons stops firing altogether. The stability state of locking is much more difficult to achieve and infact we found distinct branches in their rotation number, these states should all be regarded as variants of 1:1 locking, and therefore we lump them together in the locked region of the stability diagram. With further increase in *I*, parameter heterogeneity,

**FIGURE 6 | Phase diagram of mean synaptic action variable is shown as a function of 2D parameter space of stimulus spread** *-I* **and excitatory coupling strength** *K***.** In the partial burst regime labeled as PO in **(A)**, a subset of neurons are not firing at all as their respective drives are below their individual firing thresholds. Heterogeneous dynamics between synchronized population spiking activity and oscillation frequency death response for PO state is displayed in **(B)**. Nearly

successively more neurons peel away until eventually the entire population is IN.

## **4. PHASE DIAGRAM USING SEMI ANALYTICAL METHODS FOR MEAN FIELD PHASE MODEL**

Mean field coupled neurons in phase model is described in Equation (6). Let's rewrite the mean field equation explicitly.

$$\dot{\theta\_i} = (F(\theta\_i) - \Gamma \sin \theta\_i (\cos \theta\_i - V\_{\text{th}})) \tag{15}$$

$$\text{where } \Gamma = \frac{K}{N} \sum\_{j=1}^{N} \frac{\beta}{\left(1 + \beta + \exp\left(-\frac{\cos \theta\_j}{2}\right)\right)}.$$

In a semianalytical approach we would like to understand the phase transitions between three distinct network states discovered in two networks. For the IN states where the average firing frequency increases monotonically plotted in **Figure 10**, the θ*i* are all distributed across the closed orbit in a unit circle. This leads to the following phase evolution equation SR state may undergo instability either through parameter changes of *K* or *I* and make phase transition to either IN or PO state. Mean field approaches a stationary density as the number of neurons are increased in both PO and IN state (see **Figure 4D**). Hence,

synchronized population of spike/burst activity lumped in a regime labeled as SR [corresponding time series is displayed in **(C)**] and incoherent population spike/burst activity is lumped into a regime called IN. In the incoherent regime mean field values stay close to zero with substantial subthreshold fluctuations. In **(D)** multi stability of PO state is displayed again; now between population burst and fixed point dynamics for an entirely different parameterization.

 approaches some positive real number for these two states. When varying *K*, we consider small perturbations μ to the SR solution θ = θ*i* = 0. With θ = θ*i* = 0 + μ Equation (15) becomes θ˙ *i* = ˙μ = *F*(μ) + 2 sin (2μ) and linearization yields θ˙ *i* = (*F* (0) + )μ. Moreover, SR state may gets phase locked at θ = π (subpopulation clusters). Hence, θ = π may get destabilized as we changed the width of heterogeneity by changing *I* or the coupling strength *K*. Similarly, we consider small perturbations μ about solution θ = π. Hence, we can write θ*i* = π + μ, θ˙ *i* = ˙μ = *F*(π + μ) + 2 sin(2π + 2μ) and linearization yields θ˙ *i* = (*F* (π) + )μ. With *F*(θ) ≈ *I* − cos θ − cos θ *n* for the SR state, we find that *F* (θ) ≈ − sin(θ) − 1 *n* sin θ *n* will be generally small for θ = 0,π. SR state solution hence becomes unstable when *F* (0,π) + = 0, which suggest almost a vertical critical line between SR and PO, IN state. The bifurcation route from PO (multistable state) to IN solutions as the parameter *I* increases is less conclusive in the framework of the circular approximation, since in the previous stability analysis the only *I*-dependent term is *F* (0,π), which is very small, hence higher orders of the approximation must be considered. We use the following ansatz: If *r* is the radius of a unit circle, any smooth deformation from a unit circle can be approximated as, *r*(θ*i*) = 1 + *h*(θ*i*). Hence we can

**neurons.** All neurons are spiking in synchrony and time locked. The parameter values are *K* = 0.8, *I* = 0.001.

compute the non-linear flow contribution with the above first order correction term as *F*(θ*i*) + *H*(θ*i*). It is possible to explicitly determine *H*(θ*i*) for a certain choice of *h*(θ*i*) and moreover, *H*(θ*i*) has a periodicity of π, that is *H*(θ + π) = *H*(θ). Thus the linear stability analysis about the fixed point θ∗ *i* = 0 + μ gives

$$\dot{\mu} = (F'(0) + \epsilon H'(0) + \Gamma)\mu = (F'(\pi) + \epsilon H'(\pi) + \Gamma)\mu \quad (16)$$

From the above equation with *F* (0) ≈ 0 and the π-periodicity of *H*(θ*i*), we find that the two fixed points at 0,π lose stability at the same time for increasing *I* and as a result leads directly to the IN state. Since *H* (θ*i*) ∼ *I*, scales linearly for fixed μ, we can also estimate the critical line of transition in the parameter space in **Figure 11** which separates PO state from IN state. For the critical line: *H* (θ*i*) = *m*(*I* − *Ic*)*p* where *m* is the slope of this line and *m* > 0 allows for destabilization. Hence the critical condition is *H* (0) + = 0. By substituting the dependence of *H* (θ*i*) on (*I*, *Ic*) and in turn dependence on *I* one can write *m*(*I* − *Ic*)*p* + = 0. This implies coupling strength *K* = −(*I* + *Ic*)*p* for (*m* > 0) and *p* is some exponent representing a scaling relationship near saddle-node bifurcation. Thus the critical condition is |*K*| = *m*(*I* − *Ic*), which serves as a convenient guide to numerically compute the stability line separating PO region from IN. Next we try to obtain analytically the stability boundary between INC and PO oscillation states in the infinite-N limit. it turns out that the IN and partial oscillation states can be made steady in our system. The possibility of doing so was suggested by the numerical results. In numerics we observed that as the number of neurons *N* is increased, the order parameter *Rt* approaches a constant for both these states **Figure 8** and the oscillators tend to arrange themselves in a stationary distribution around the circle **Figure 11**. The way to approach these two states analytically, therefore, is to first write down the appropriate infinite-N analog of our model.

$$\frac{\partial}{\partial t}f + \frac{\partial}{\partial \theta} \left\{ \left[ F(\theta) - \left( \left( \int\_0^{2\pi} \int\_{I-\Delta I}^{I+\Delta I} \Gamma p(\theta', t, I') g(I') dl' d\theta' \right) \right) \right] \right. \\ \left. + \sin(2\theta) \right] \rho \right\} = 0 \tag{17}$$

$$\times \sin(2\theta) \left[ \rho \right] = 0 \tag{17}$$

The above equation is the infinite-N analog of continuity equation introduced earlier. It is a non-linear partial integrodifferential equation for the number density *f*(θ,*t*, ω). In addition we demand *f* to be non-negative, 2π periodic in θ, and we impose the normalization

$$\int\_0^{2\pi} f(\theta, t, I) d\theta = 1,\tag{18}$$

For incoherence and partial oscillation the above system tends toward a stationary distribution of phases in time. The above two states are the fixed points of the stationary density in the continuum limit. To solve for the fixed points we set ∂ ∂*t f* = 0 in Equation (10). let's assume that *f*0(θ, ω) be the stationary phase density and *v*0 = [*F*(θ) − (( 2π 0 *I* + *I I* − *I f*(θ ,*t*,*I* )*g*(*I* )*dI d*θ ))sin(2θ)] be the velocity field. Then one can write

$$\frac{\partial}{\partial \theta} (f\_0 \nu\_0) = 0 \Rightarrow f\_0 \nu\_0 = L(I) \tag{19}$$

where *L*(*I*) is a constant which is determined exactly by using normalization condition. Depending on it's applied drive *I*, neuron's steady state behavior falls in the following two categories:

Case (i) When *I* << implies

$$\nu\_0(\theta, I) = F(\theta) - \Gamma \sin(2\theta) = 0 \tag{20}$$

corresponding time evolution of *Vi*(*t*) is shown for all 10 bursting neurons. In network. In **(E)** standard deviations of *V*mean is plotted as a function of time. In **(F)** mean field vs. time for the coherent state is shown. **(G)** displays temporal dynamics of order parameter. Average firing frequency as described in Equation (14) is plotted in **(H)** for the parameter combination of *K* = 0.001, *I* = 0.15. Panel **(H)** further demonstrates phase locking behavior among all

$$f\_0(\theta, I) = \frac{L(I)}{F(\theta) - \Gamma \sin(2\theta)}\tag{22}$$

In the IN state, neurons driven by different external drives are firing at different phases, however, their collective state is close to being stationary. Every neuron belong to Case (ii) as described above. Further, it is possible to derive nearly an exact relationship between *K*, *I* that gives the transition from case (i) to

Case (ii) When *I* >> neuron fires incoherently and typically individual phases follows an uniform distribution about the unit circle. In this case the velocity field turns out to be,

**(D)** temporal response of *V*mean is shown for the coherent state of our

$$\mathcal{W}(\theta, I) = F(\theta) - \Gamma \sin(2\theta) \tag{21}$$

Fixed point solution demands that the density must be inversely proportional to the velocity:

the neurons.

case (ii) as described above. As shown before, in case of a finite size network such a relationship in the first order perturbation *m*(*I* − *Ic*)*p* + = 0 does exist. In this scenario those neurons with a minimum bound on their applied drive *I*min reach cessation of firing as we find from numerical simulations. They then fall into the Case (i) above where mean field exerts much bigger influence on the dynamics and overall effect is damping of firing activity. The first neurons to stop firing are the ones which do not cross the threshold for firing which in this case *I* > 2. Then the boundary that separates IN from PO in the phase diagram is almost a straight line given by,

$$|K| = \epsilon m(\Delta I - \Delta I\_{\mathcal{E}}) \tag{23}$$

Hence, both finite and infinite analog of our network identifies the putative transition boundary between IN and PO states. Now from numerical simulations we find Andronov– Hopf (AH) bifurcations leads to the transition from INC to SR solutions in the **Figure 6** near *K*, *I* values close to zero. It is equivalent to look at the imaginary eigensolutions that arise due to the instability of the IN state. This instability requires calculation of higher order perturbation terms of the stationary density obtained at the IN state of our network. This is out of the scope of our paper, however, we show a numerical fitting result which gives an empirical relationship between *K* and *I* to quantify the transition boundary between IN and SR states. Assuming is the perturbation to the IN

solution we can express a relationship between *K* and *I* as follows,

$$|K| = a\_0 \epsilon + a\_1 \epsilon^2 + a\_2 \epsilon^3 + O(\epsilon^4) \tag{24}$$

Equation (24) gives us an empirical relationship between parameters upto fourth order perturbations for the bifurcation of a limit cycle. Optimization of the above equation gives coefficients *a*0 = 8 π, *a*1 = 0, *a*2 = 128 π3 , respectively. Next, we substitute the amount of dispersion *I* into the perturbative term to obtain the boundary between IN and SR state. Taken together we can write,

$$|K| = \frac{8}{\pi} \Delta I + \frac{128}{\pi^3} \Delta I^3 + O(\Delta I^4) \tag{25}$$

Results are shown in **Figure 12** in the (*K*, *I*) plane using Equations (25) and 23. Critical lines obtained semi analytically qualitatively agrees well with the numerical results that captures various network states in both these models with purely excitatory coupling. In Appendix, we show a stability calculation for an inhibitory coupled mean field network in the infinite analog limit. From numerical simulations we find that the results are independent of the number of spikes *n* per burst.

# **5. DISCUSSION**

One of the most frequent assumption employed in simulations of large neural networks is that the whole network can be lumped into small aggregates of collective unit (sometimes called a "neurocomputational unit") exhibit a sufficiently similar dynamical behavior. Consequently, the network that instantiates this ensemble, consisting of thousands of excitatory and inhibitory neurons, it is considered to display a synchronized behavior with no other significant temporal features for the dynamics of the large scale network. The main reason for this assumption, is the impractical large computational time arising from too many details considered in the large network properties. In this paper, we

have analyzed the behavior of a neural network that serves as a good example of such a unit, namely a mean field coupled bursting ensemble. First, we have investigated a Hodgkin–Huxley type detailed biophysical model widely employed in theoretical and computational neuroscience with global coupling. We found that the dynamical features of the network are far more complex then the ones corresponding to synchronized or rest state behavior. The network dynamics depends critical on the balance between firing rate threshold dispersion and mean field synaptic coupling strength; in fact, the synchronized state can be found only for a specific range of parameters typically involving a large or medium values for the coupling strength and low values of dispersion. On the other hand, for large dispersion and weak coupling strength values both networks display purely IN behavior. In the IN state, individual neurons are driven by different external drives results in firing at different phases, however, their collective state is close to being stationary. This stationarity in the density distribution led us to formulate scaling relationship between coupling strength and dispersion parameter. One interesting finding is that, when mean field exerts a greater influence than parameter dispersion; it causes shutting down of the neural activity in some neurons. In this parameter range, we find interesting dynamical behavior such as partial activity. In order to address the problem of the high computational cost of such an implementation, we have further developed a selfconsistent mathematically tractable mean field coupled phase model following (Assisi et al., 2005; Ghosh et al., 2009; Jirsa and Stefanescu, 2011), but incorporating a higher degree of realism. Rather than finding the most appropriate type and number of dimensions that could minimize certain error functions or capture statistical variance in the full network, we have focused our attention on understanding a phenomenological burst generation model system which captures the most important network dynamics of bursting units at the population level. Collective activity of synapses is described by a mean field which relies on instantaneous rise and decay time (Roy et al., 2011). This mean field is then employed in the coupling to individual neurons to describe phase network. Together, we investigate this population of neurons coupled to a common mean field drive and heterogeneity in their threshold for spikes/bursts. Our detailed analysis demonstrated that the reduced representation manages to recreate correctly the topology of the mean field amplitudes of the original system for various parameter scenarios. In the full network, In the thermodynamic limit (*N* → ∞), a collective state becomes coherent if δ*V*mean(*t*) ≡ [*V*mean(*t*) − *V*mean(*t*)] is non-stationary (i.e., an oscillating global potential *V*mean appears for a coherent case) and also, the correlated mean field (*t*) activity appears oscillatory. In the phase network, global order parameter is computed by averaging the contributions of all microscopic spikes within a burst in order to obtain a similar degree of ordering of spikes/bursts as in the full model for identical parameterization. Hence, for a dynamical behavior such as coherence-incoherence transition macroscopic order parameter gives us a crude approximation of burst timing. From a more general perspective, despite its limitations because of the consideration of purely excitatory or inhibitory network, it can be readily extended to study networks with mixed coupling. Moreover, the analytical approach to estimate the scaling relationship and transition boundaries between the IN-PO-SR states is not restricted to small scale network only. With global coupling, as the size of the network grows the boundaries may shift leading to a different parameterization than the one displayed here; however, underlying bifurcations remain the same. We have demonstrated this in our work by analytically deriving a low dimensional mean field amplitude reduction for a inhibitory coupled mean field network in the continuum limit. In this case, all the relevant dynamics of an infinite dimensional network in Equations (29) and (30) is captured by a two dimensional representation of the reduced mean field population given by Equation (40). Thus, using this approach, we derive analytically a low dimensional representation of the network dynamics and we show that the main features of the neural population's collective behavior can be captured well by the dynamics of a few cortical nodes exhibiting spiking as well as bursting behavior. While it is true that strong reductionist assumptions are common (sacrificing many of the biological realism of a network node's dynamics) in large-scale network modeling, these assumptions are usually made *ad-hoc* on the network node's dynamics and limit the network dynamics to a small range. We emphasize here that because of the "near to synchrony" assumption, neural mass models cannot capture complex dynamical features such as multiclustering, oscillator death or multi-time scale synchronization. Evidently a reduced small scale network model is desirable to serve as a node in a large scale network simulation whereby displaying a sufficiently rich dynamic repertoire. Here it is of less importance to find a quantitatively precise reduced description of a neural population; rather more importantly, we seek a computationally inexpensive population model (this means typically low-dimensional) which is able to display the major qualitative dynamic behaviors (synchronization, rest state, multi-clustering, etc.) for realistic parameter ranges as observed in the total population of neurons. Our approach may offer a viable alternative to the neural mass models currently used in the literature. By comparison, our model offers the possibility to account for such features (temporal details of their spiking activity considered irrelevant for the dynamics of the large network) at a very low computational cost. Therefore, the type of reduced representation discussed in this paper qualifies as a good candidate for a "neural unit" in computational simulations of large scale neural networks.

# **REFERENCES**

# **ACKNOWLEDGMENTS**

We would like to thank Anandamohan Ghosh and Mohit Adhikari for the helpful corrections on this manuscript. The research reported herein was supported by the Brain Network Recovery Group through the James S. McDonnell Foundation and the FP7-ICT BrainScales.

When inhibition not excitation synchronizes neural firing. *J. Comput. Neurosci.* 1, 313–321.

model neurons. *Neural Comput.* 4, 84–97.

commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 28 January 2013; accepted: 05 March 2013; published online: 26 March 2013.*

*Citation: Roy D and Jirsa V (2013) Inferring network properties of cortical neurons with synaptic coupling and parameter dispersion. Front. Comput.* *Neurosci. 7:20. doi: 10.3389/fncom. 2013.00020*

*Copyright © 2013 Roy and Jirsa. 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 thirdparty graphics etc.*

# **APPENDIX**

## **MEAN FIELD REDUCTION FOR INHIBITORY SYNAPTIC COUPLING**

Here, we extend our network in the continuum limit in the presence of inhibitory coupling. *V*th is held negative. (*N* → ∞), where the state of the coupled system can be described by a density function *f*(θ,*I*,*t*), where *f* is defined such that the fraction of neurons with phases lying between θ and *d*θ and applied drive between *I* and *dI* is given by *f*(θ,*I*, *t*)*d*θ*dI* (Antonsen et al., 2008; Ott and Antonsen, 2009). The applied stimulus are drawn from a distribution *g*(*I*) such that

$$\int\_{-\infty}^{\infty} \int\_{0}^{2\pi} f(\theta, I, \mathfrak{t}) d\theta dI = 1 \tag{26}$$

$$\int\_0^{2\pi} f(\theta, I, t)d\theta = \lg(I) \tag{27}$$

For the conservation of currents *I* the continuity equation is written as

$$\frac{\partial f}{\partial t} + \frac{\partial (f \,\nu)}{\partial \theta} = 0.\tag{28}$$

In order to make the coupling amenable to analytical study we use a pulse-like function for the mean field = *a*1 + *b*1(1 + cos θ). Response to the mean field by individual neuron's *R*(θ) = *V*th sin(θ), containing only single Fourier component, a choice motivated primarily due to the tractability of the resulting model. Further, *V*th = −1 for the convenience of calculations without losing any generality of our results. The velocity *v*(θ,*I*, *t*) in Equation (28) is now written as

$$\nu(\theta, I, t) = a + \epsilon \int\_{-\infty}^{\infty} \int\_{0}^{2\pi} (1 + \cos(\hat{\theta}))(-\sin \theta) f(\hat{\theta}, \hat{I}, t) d\hat{\theta} d\hat{I}$$

$$-F \sin(\theta) - F \sin(\theta/n) \tag{29}$$

where without loss of any generality we are using sin functions instead of cos functions in Equation (6).

In the continuum limit the order parameter *z* can be defined as

$$z(t) = \int\_{-\infty}^{\infty} \int\_{0}^{2\pi} \frac{(e^{i\theta} + e^{-i\theta})}{2} f(\theta, I, t) d\theta dI \tag{30}$$

It's a linear sum of two complex order parameters and one could in principle unfold the entire dynamics of the network in any one of the manifold given above. Here,

$$z1(t) = \int\_{-\infty}^{\infty} \int\_{0}^{2\pi} e^{i\theta} f(\theta, I, t) d\theta dI \tag{31}$$

$$z2(t) = \int\_{-\infty}^{\infty} \int\_0^{2\pi} e^{-i\theta} f(\theta, I, t) d\theta dI \tag{32}$$

Using the above it is easy to see that the expression for velocity becomes

$$\nu(\theta, I, t) = I + \frac{1}{2i} \left[ \left( 1 + \frac{\epsilon}{2} (z2 + z1^\*) + F \right) e^{-i\theta} \right]$$

$$- \left( 1 + \frac{\epsilon}{2} (z1 + z2^\*) + F \right) e^{i\theta} \Big|$$

$$- F \left( \frac{e^{i\theta/n}}{2i} - \frac{e^{-i\theta/n}}{2i} \right) \tag{33}$$

∗ indicates the complex conjugate. The distribution function can be expressed as a Fourier series

$$f(\theta, I, t) = \frac{\mathbf{g}(I)}{2\pi} \left[ 1 + \sum\_{k=1}^{\infty} f\_k(I, t) e^{ik\theta} + c.c. \right] \tag{34}$$

The above infinite dimensional system is difficult to analyze. However, the "amazing" anstaz of Ott and Antonsen (2009) has been shown to be successful in obtaining the low-dimensional description of the globally coupled phase oscillators. The anstaz impose a restriction on the fourier coefficients:

$$f\_k(I, t) = \left(\psi(I, t)\right)^k\tag{35}$$

for *k* ≥ 1 and has been shown to be a reasonable guess under different scenariors (Ott and Antonsen, 2009). This restricted class of functions readily reduces our continuity equation to an θ-independent form

$$\frac{d\Psi}{dt} = \frac{1}{2}\left(1 + \frac{\epsilon}{2}z1 + F\right)^{\*} - iI\psi$$

$$-\frac{1}{2}\left(1 + \frac{\epsilon}{2}z1 + F\right)\Psi^{2}$$

$$-F\left(\frac{\Psi^{1 + 1/n}}{2} - \frac{\Psi^{1 - 1/n}}{2}\right) \tag{36}$$

with *z*1 satisfying

$$z1(t) = \int\_{-\infty}^{\infty} \Psi^\*(I, t)\varrho(I)dI. \tag{37}$$

If we assume that *g*(*I*) is a Lorentzian distribution function

$$\lg(I) = \frac{1}{\pi[(I - I\_0)^2 + 1]}.\tag{38}$$

*z*(*t*) can be evaluated by contour integration with poles at *I* = *I*0 − *i* and we obtain the exact evolution equation of order parameter *z*

$$\frac{dz1}{dt} = iI\_0 z1 - z1 + \frac{1 + \frac{\epsilon}{2}z1 + F}{2} - \frac{1 + \frac{\epsilon}{2}z1^\* + F}{2}$$

$$-F\left(\frac{z1^{1 + 1/n}}{2} - \frac{z1^{1 - 1/n}}{2}\right) \tag{39}$$

The above equation can be expressed in polar coordinates if we substitute *z*1 = ρ1exp(*i*φ1) giving evolution equations for ρ1 and φ1

$$\frac{d\rho\_1}{dt} = \frac{\epsilon}{2}\rho\_1\left(1-\rho\_1^2\right) - \rho\_1 + \frac{F}{2}\left(1-2\rho\_1^2\right)\cos\phi\_1 + \frac{1}{2}\cos\left(\phi\_1\right)$$

$$+ \frac{F\rho\_1}{2}\left(\rho\_1^{-1/n} - \rho\_1^{1/n}\right)\cos(\phi\_1/n) \tag{40}$$

$$\frac{d\phi\_1}{dt} = I\_0 - \frac{F}{2} \left(\rho\_1 + \frac{1}{\rho\_1}\right) \sin \phi\_1 - \frac{\rho\_1}{2} \sin(\phi\_1)$$

$$-\frac{F}{2} \left(\rho\_1^{1/n} + \frac{1}{\rho\_1^{1/n}}\right) \sin(\phi\_1/n). \tag{41}$$

For the Lorentzian distribution function the above equation is exact, However, we do not find any deviation of the above results for any other unimodal distributions of our firing threshold (such as uniform distribution) such as the one considered in the numerical simulations with excitatory coupling. The above two dimensional system can be solved numerically to identify the full network states and the corresponding transition boundaries.

# Ketamine, propofol, and the EEG: a neural field analysis of HCN1-mediated interactions

# **Ingo Bojak 1,2\*, Harry C. Day 1 and David T. J. Liley 3,4**

1 Centre for Computational Neuroscience and Cognitive Robotics, School of Psychology, University of Birmingham, Birmingham, UK

2 Donders Centre for Neuroscience, Donders Institute for Brain, Cognition and Behaviour, Radboud University Nijmegen Medical Centre, Nijmegen, Netherlands

3 Brain and Psychological Sciences Research Centre, Faculty of Life and Social Sciences, Swinburne University of Technology, Hawthorn, VIC, Australia

4 Cortical Dynamics Ltd., Hawthorn, VIC, Australia

### **Edited by:**

Dimitris Pinotsis, University College London, UK

### **Reviewed by:**

Julien Modolo, Lawson Health Research Institute and Western University, Canada Sacha J. Van Albada, Research Center Jülich, Germany Jamie Sleigh, University of Auckland, New Zealand

#### **\*Correspondence:**

Ingo Bojak, Centre for Computational Neuroscience and Cognitive Robotics, School of Psychology, University of Birmingham, Birmingham West Midlands B15 2TT, UK. e-mail<EMAIL_ADDRESS>
Ketamine and propofol are two well-known, powerful anesthetic agents, yet at first sight this appears to be their only commonality. Ketamine is a dissociative anesthetic agent, whose main mechanism of action is considered to be N-methyl-d-aspartate (NMDA) antagonism; whereas propofol is a general anesthetic agent, which is assumed to primarily potentiate currents gated by γ-aminobutyric acid type A (GABAA) receptors. However, several experimental observations suggest a closer relationship. First, the effect of ketamine on the electroencephalogram (EEG) is markedly changed in the presence of propofol: on its own ketamine increases θ (4–8 Hz) and decreases α (8–13 Hz) oscillations, whereas ketamine induces a significant shift to beta band frequencies (13–30 Hz) in the presence of propofol. Second, both ketamine and propofol cause inhibition of the inward pacemaker current Ih, by binding to the corresponding hyperpolarization-activated cyclic nucleotidegated potassium channel 1 (HCN1) subunit. The resulting effect is a hyperpolarization of the neuron's resting membrane potential. Third, the ability of both ketamine and propofol to induce hypnosis is reduced in HCN1-knockout mice. Here we show that one can theoretically understand the observed spectral changes of the EEG based on HCN1-mediated hyperpolarizations alone, without involving the supposed main mechanisms of action of these drugs through NMDA and GABAA, respectively. On the basis of our successful EEG model we conclude that ketamine and propofol should be antagonistic to each other in their interaction at HCN1 subunits. Such a prediction is in accord with the results of clinical experiment in which it is found that ketamine and propofol interact in an infra-additive manner with respect to the endpoints of hypnosis and immobility.

**Keywords: ketamine, propofol, EEG, HCN1, neural field theory, drug interaction, anesthesia, infra-additivity**

# **INTRODUCTION**

Ketamine, a phenylcyclohexylpiperidine (PCP) derivative, is a powerful psychoactive drug that is predominantly used as a sedative and general anesthetic agent in humans and animals (Sinner and Graf, 2008). Ketamine occurs as two stereoisomers, R(−) and S(+), in which the latter is found to be some three to four times more potent (White et al., 1985), but despite such differences in potency the drug is generally made available clinically as a racemate (racemic mixture) that contains both stereoisomers in equal proportion. Ketamine is classified as a dissociative agent due to its ability to induce hallucinations and perceptual/environmental detachment (Wolff and Winstock, 2006). Because of these properties it has become popular recreationally. At sufficiently high doses it has been reported to induce a state of dissociation comparable to that of schizophrenia, and as a consequence has found use as a pharmacological model for psychosis (Bubenikova-Valesova et al., 2008; Corlett et al., 2011). More recently its therapeutic use has been re-evaluated in light of evidence suggesting that sub-anesthetic doses may aid in the treatment of bipolar affective disorder and major depression (Mathew et al., 2012; Murrough, 2012; Murrough et al., 2012).

While ketamine is widely believed to act principally through the non-competitive antagonism of *N*-methyl-d-aspartate (NMDA) receptor mediated glutamatergic activity (Irifune et al., 1992; Oye et al., 1992), two significant pieces of empirical evidence have emerged that challenge such a unitary view. Firstly, dizocilpine (also known as MK801), an even more potent non-competitive NMDA antagonist, produces no significant hypnotic effect (Kelland et al., 1993;Irifune et al., 2007). Secondly, ketamine's effect on spontaneous electroencephalogram (EEG) activity is qualitatively altered when administered in the presence of propofol, a widely used intravenous general anesthetic agent that, at clinically meaningful concentrations, has little or no effect on NMDA mediated currents. Ketamine alone has been shown to reduce spectral edge frequencies, an effect that is driven predominately by increases in absolute θ band (4–8 Hz) power at the expense of α band (8–13 Hz) power (Schuttler et al., 1987; Kochs et al., 1996), see **Figure 1A**. In contrast, ketamine administered in the presence of steady state propofol levels is associated with a definite acceleration of α band activity; increasing its peak frequency by up to 4.7 Hz (Hayashi et al., 2007; Tsuda et al., 2007), see **Figure 1C**. Propofol on its own roughly maintains the α peak frequency with an

Average spectra of EEG recorded during two sequential target concentrations of propofol in a single subject. BL = baseline, M = 1.25µg/ml propofol, H = 2.5µg/ml propofol. α band EEG recorded from parietal (P4) and occipital (O2) electrodes reveals minimal changes in peak frequency with increasing

anteriorization of power (decrease occipital, increase frontal), see **Figure 1B**; though an additional broadband "beta buzz"just above α frequencies, "biphasic" response dynamics and smooth transitions to lower frequencies can confound the picture (Schwender et al., 1996; Kuizenga et al., 1998, 2001; Feshchenko et al., 2004; Breshears et al., 2010; Cimenser et al., 2011). We assume here from previous theoretical studies (Liley et al., 2003; Hutt and Schimansky-Geier, 2008; Hutt and Longtin, 2010; Hindriks and van Putten, 2012) that these complications can be accounted for by mechanisms not considered in this work, in particular the prominent γ-aminobutyric acid type A (GABAA) agonism of propofol that affects dominantly inhibitory postsynaptic currents (Kitamura et al., 2003). Furthermore, the acceleration due to ketamine observed by Hayashi et al. (2007) and Tsuda et al. (2007) that we wish to describe occurred on top of a clear α rhythm at steady propofol concentration, see **Figure 1C**. Thus we assume in the following that the action of propofol is largely neutral concerning the α peak frequency (while unspecified concerning total α band spectral power).

Recently a number of alternative, behaviorally relevant, molecular targets for ketamine action have been identified (Schnoebel et al., 2005; Hevers et al., 2008; Chen et al., 2009). Of particular significance is the identification of hyperpolarization-activated cyclic nucleotide-gated (HCN) potassium channel subunits as a target for ketamine action (Chen et al., 2009). HCN subunits, of which there are four isoforms (HCN1–4), assemble to form a tetrameric ion channel that mediates an inward (i.e., depolarizing) hyperpolarization-activated pacemaker current *I* h implicated in neuronal rhythmogenesis (Biel, 2009; Biel et al., 2009). In

Tsuda et al. (2007). particular the HCN1 isoform has been identified as a molecular substrate for the actions of ketamine (Chen et al., 2009): ketamine causes inhibition of HCN1-mediated *I* h currents, and hence membrane hyperpolarization, in pyramidal neurons from wildtype but not HCN1-knockout mice. The potency of ketamine to provoke a loss of the righting reflex (the ability to regain footing from a back position), which is a behavioral correlate of hypnosis, is also strongly reduced in HCN1-knockout mice. Hence a causal relationship between ketamine-induced membrane hyperpolarization and its clinical effects can be made. The existence of such a causal relationship is made more likely by evidence indicating that the hypnotic potency of propofol is also reduced in HCN1-knockout mice, in approximate proportion to its ability to

inhibit HCN1-mediated membrane depolarization (Chen et al.,

band EEG activity is markedly shifted to higher frequencies. Data shown is mean PSD recorded at Fp1-A1, with an Fpz reference, from nine subjects undergoing elective abdominal surgery. Figure used with permission from

2009). It should be noted though that the hypnotic response was not abolished entirely in HCN1-knockout mice by either ketamine or propofol (Chen et al., 2009), thus other effects like the mentioned GABAA agonism will be required to fully understand the hypnotic action of these agents. However, etomidate, which has no effect on HCN1 channels, showed no loss of hypnotic effect in HCN1-knockout mice (Chen et al., 2009), suggesting the specific involvement of HCN1-mediated *I* h currents for ketamine and propofol. The identification of a shared molecular target for ketamine and propofol action could offer a new possibility to account for the qualitatively disparate electroencephalographic effects of ketamine alone and in the presence of propofol. On this basis we hypothesized that by modeling the differential effects of ketamine and propofol on neuronal membrane hyperpolarization, in the context of an established theory of resting EEG (Liley et al., 2002, 2010, 2011;Bojak and Liley, 2005),we would be able to describe the observed effects on the EEG at least qualitatively. Because current depth of anesthesia monitoring approaches are either insensitive (Faraoni et al., 2009; Nonaka et al., 2012), or respond anomalously (Hans et al., 2005; Sengupta et al., 2011), to the hypnotic effects of ketamine, understanding the mechanism by which ketamine and propofol interact electroencephalographically will ultimately assist in the development of improved approaches to clinically monitor the hypnotic effects of combinations of these drugs. The combination of propofol and ketamine (often referred to as *ketofol*) is becoming increasingly important in the procedural sedation setting where rapid and effective sedation and analgesia, with minimal cardiorespiratory/hemodynamic compromise, is required (Hui et al., 1995; Frizelle et al., 1997; Sakai et al., 1999; Phillips et al., 2010).

## **MATERIALS AND METHODS**

#### **MODELING DRUG RESPONSE AND INTERACTIONS**

The simplest pharmacodynamic model of drug effect involving two or more agonists is that of competitive ligand-receptor binding. It is easily shown for two full agonists competing for the same receptor binding site, that the fractional receptor occupancy θ, as a function of the respective drug concentrations (*D*1, *D*2) is (Shafer et al., 2008)

$$\theta = \frac{k\_2 D\_1 + k\_1 D\_2}{k\_2 D\_1 + k\_1 D\_2 + k\_1 k\_2}, \text{ or } \tag{1}$$

$$\frac{\theta}{1-\theta} = \frac{D\_1}{k\_1} + \frac{D\_2}{k\_2}.\tag{2}$$

For *D*1 →∞ and/or *D*2 →∞, one then finds θ→1, i.e., full receptor occupancy. *k*1, *k*2 > 0 are the respective drugreceptor dissociation constants, which are equivalent to single drug concentrations that produce 50% receptor occupancy, i.e., θ = 1/2. In general a pharmacodynamic effect *E* is assumed to be some monotonic function of θ, i.e., *E* = *f*(θ) ≡ *g* (*D*1, *D*2). For a fixed effect *E* the locus of points (*D*1, *D*2) defines a *response isobole* and *E* = *g* (*D*1, *D*2) a *response surface*. Now consider the case of competitive binding and drug interaction (Greco et al., 1995)

$$\frac{\theta}{1-\theta} = \frac{D\_1}{k\_1} + \frac{D\_2}{k\_2} + \frac{\eta D\_1 D\_2}{k\_1 k\_2},\tag{3}$$

where η defines an interaction term. It can be easily demonstrated that

$$\begin{array}{rclcrcl}\frac{\theta-1}{\theta} & < & \eta & < & 0 & \Rightarrow & \text{infra-additivity / antagonism,} \\ & & \eta & = & 0 & \Rightarrow & \text{additivity,} \\ & & \eta & > & 0 & \Rightarrow & \text{synergy.} \end{array} \tag{4}$$

Inspired by these considerations, we chose here to describe the general pharmacodynamic effect of our two ligands by the following bilinear form

$$E = c\_1 D\_1 + c\_2 D\_2 + c\_{12} D\_1 D\_2.\tag{5}$$

This ansatz represents the simplest extension beyond the purely additive; and the sign of *c*12 then has the same interpretation as the sign of η in Eq. 4. We will use this bilinear form below to parameterize the dependence of the induced hyperpolarizations on normalized concentrations of propofol and ketamine, respectively.

One can however relate Eqs 3 and 5 more directly. Assume first that the pharmacodynamic effect is directly proportional to receptor occupancy, i.e., *E* ∝ θ. Then *k*1 and *k*2 become the respective "half maximum effective concentrations" (EC50s) at which 50% of the maximum response is observed for each drug applied alone. Furthermore, assume that the receptor occupancy remains relatively small θ ≈ θ/(1 − θ) = *D*1/*k*1 + *D*2/*k*2 + η*D*1*D*2/(*k*1*k*2) < 1/2, so that the effect *E* < *E*max/2. The half-maximal inhibition of HCN1 subunit-mediated ionic currents by racemic ketamine occurs at a concentration of approximately 16µM (Chen et al., 2009), which is significantly greater than the estimated minimum free plasma concentrations of 2.9µM required to produce surgical anesthesia in humans (Grant et al., 1983). Data for the half-maximal inhibition of HCN1-mediated ionic currents by propofol is to our knowledge not available. However, because HCN1-knockout mice are significantly less sensitive to the effects of propofol than wild-type ones, we can speculate that the ED50 (the "half maximum effective dose") for unresponsiveness with propofol in wild-type mice corresponds roughly to the half maximum of the neuronal changes (EC50). Chen et al. (2009) found this to be approximately 7 mg/kg. Using the volume of distribution of 1.38 l/kg (Cox et al., 1998) in the rat (murine values not available), EC50 is then about 5.1 mg/l or 29µM, which is significantly greater than the minimum free plasma concentration ∼8.5µM for surgical anesthesia. Thus *D*ketamine/*k*ketamine + *D*propofol/*k*propofol +η*D*ketamine*D*propofol/(*k*propofol*k*ketamine) < 1/2 is approximately satisfied as long as η < 1 2 . While the *E*max for ketamine-induced membrane hyperpolarization in murine pyramidal neurons is of the order of −4 mV, the actual value of *E*max will depend on the species and the recording conditions/preparation. In the absence of any information to the contrary one can assume that *E* < *E*max/2. Thus our ansatz Eq. 5 can be considered as following from Eq. 3 under a range of reasonable assumptions.

#### **LILEY MODEL AND EIGENSPECTRUM CALCULATION**

We base our investigation in this paper on the Liley et al. (2002) model, which is a typical neural field model (Deco et al., 2008; Coombes, 2010; Bressloff, 2012; Liley et al., 2012). In Bojak and Liley (2005) 73,454 different parameter sets, which produce biologically plausible resting state activity, were found for this model. We use here also the "eigenspectrum" approach introduced in Bojak and Liley (2005) to directly predict EEG power spectral densities (PSDs) from a given parameter set. In the following we will briefly review a few key features of the Liley et al. (2002) model and of the eigenspectrum approach that will play a role for the analysis in this paper, and refer the reader to the original reference for more detail. The Liley et al. (2002) model can be written concisely as follows:

$$\pi\_k \frac{\partial}{\partial t} h\_k \left( \mathbf{x}, t \right) = h\_k^r - h\_k \left( \mathbf{x}, t \right) + \sum\_{l=e, i} \frac{h\_{lk}^{eq} - h\_k \left( \mathbf{x}, t \right)}{\left| h\_{lk}^{eq} - h\_k^r \right|} I\_{lk} \left( \mathbf{x}, t \right), \tag{6}$$

$$\begin{split} & \left( \frac{1}{\gamma\_{lk}} \frac{\partial}{\partial t} + 1 \right)^{2} I\_{lk} \left( \mathbf{x}, t \right) \\ &= \frac{\Gamma\_{lk} e}{\gamma\_{lk}} \left\{ \frac{N\_{lk}^{\beta} S\_{l}^{\text{max}}}{1 + e^{\sqrt{2} \left[ h\_{l} (\mathbf{x}, t) - \mu\_{l} \right]} \Big/ \sigma\_{l}} + \Phi\_{lk} \left( \mathbf{x}, t \right) + p\_{lk} \left( \mathbf{x}, t \right)} \right\}, \end{split} \tag{7}$$

$$
\begin{split}
\left[\left(\frac{1}{\nu\_{\rm lk}\Lambda\_{\rm lk}}\frac{\partial}{\partial t}+1\right)^{2}-\frac{3}{2}\frac{1}{\Lambda\_{\rm lk}^{2}}\nabla^{2}\right]\Phi\_{\rm lk}(\mathbf{x},t) \\
=\frac{N\_{\rm lk}^{\alpha}S\_{l}^{\rm max}}{1+e^{\sqrt{2}\left[h\_{l}(\mathbf{x},t)-\mu\_{l}\right]/\sigma\_{l}}}.\end{split} \tag{8}
$$

In all these equations *l*, *k* = *e*, *i* serve as indices for excitatory and inhibitory neural populations, respectively, and **x** gives their position on a two-dimensional cortical sheet. The mean excitatory soma membrane potential *he*(**x**, *t*) of Eq. 6 is taken to predict the EEG. In the absence of postsynaptic inputs these potentials *hk*(**x**, *t*) decay to their resting values *h r k* . The inputs *Ilk*(**x**, *t*) correspond to postsynaptic potentials and are weighted by ionic driving forces *h* eq *lk* − *hk* (**x**, *t*), where the *h* eq *lk* are the respective Nernst potentials. These weights are normalized at rest to +1 (excitatory) and −1 (inhibitory), respectively. A postsynaptic input in Eq. 7 uses double indices to indicate source and target (for example, *Iei*(**x**, *t*) is excitatory input to an inhibitory neural population). Γ*lk* is the mean peak amplitude induced by a single presynaptic pulse δ(*t* − *t* p), and 1/γ*lk* the corresponding rise time to this peak of a postsynaptic "α form" response *I*(**x**, *t*) ∝ γ 2 *te*−γ*t* Θ(*t* − *t* p), where Θ is the Heaviside step function and δ the Dirac delta function. Extra-cortical input is given by *plk*(**x**, *t*), and is here assumed to be shaped noise (*pee*), static (*pei*), or absent (*pik*). The noise represents the average of uncorrelated input to the many neurons in the neural mass. For simplicity it is imposed only on the excitatory extra-cortical input to excitatory neurons, which is sufficient to generate the full dynamical range of the model. Finally, activity is propagated cortico-cortically via Eq. 8 with a standard damped wave equation (Jirsa and Haken, 1996; Robinson et al., 1997). The activity propagation through Φ*lk*(**x**, *t*) represents a synaptic footprint which falls off exponentially with characteristic distance scale Λ*lk*, and fibers having conduction velocity ν*lk*. Since there are no long-range inhibitory fibers, we can set Φ*ik* ≡ 0 in the following. Short range connectivity is both excitatory and inhibitory, and is represented by the first term in the curly brackets of Eq. 7. Note that Eq. 8 can be improved upon (Bojak and Liley, 2010), but its main role is in this case to include a larger variety of EEG wavelengths as will become apparent. Our main conclusions are not affected even for the radical choice of an entirely homogeneous cortex, i.e., upon removing all spatial dependence.

The eigenspectrum approach (Bojak and Liley, 2005) assumes that Eqs 6–8 have a "fixed point" solution for a homogeneous cortex with static *pee*. All variables are then linearly expanded around this solution, and auxiliary variables ˜*Ilk* = ∂*Ilk*/∂*t* and Φ˜ *ek* = ∂Φ*ek*/∂*t* are used to turn Eqs 6–8 into 14 first order ODEs. One can then Fourier-transform in space and time, and obtains an equation for the 14-dimensional state vector **s** in the form

$$i\boldsymbol{\alpha}\mathbf{s}\left(\boldsymbol{\omega},\mathbf{k}\right) = \mathbf{J}\left(\mathbf{k}^2\right) \cdot \mathbf{s}\left(\boldsymbol{\alpha},\mathbf{k}\right) + \mathbf{P}\left(\boldsymbol{\alpha},\mathbf{k}\right),\tag{9}$$

where **J** is the Jacobian matrix and **P**(ω, **k**) contains the remainder of the extra-cortical input, i.e., the variation of *pee* with subtracted mean. Note that the only spatial derivative here is the Laplacian in Eq. 8, hence the Fourier-transformed Jacobian is a function of the square of wavenumber **k**. One can then show (Bojak and Liley, 2005) that

$$|h\_{\mathfrak{e}}\left(\omega,k\right)|^{2} = \left|\sum\_{n=1}^{14} \frac{c\_{\mathfrak{n}}\left(k\right)}{i\omega - \lambda\_{\mathfrak{n}}\left(k\right)}\right|^{2},\tag{10}$$

where both the coefficients *cn* and the eigenvalues λ*n* can be obtained from a decomposition of the Jacobian in both left and right eigenmatrices, and *k* ≡ |**k**|.

Furthermore, if one makes the simplifying assumption that an EEG electrode aggregates the contributions of a disk-shaped part of the cortical sheet with radius *R*, then one can compute a prediction of the PSD as follows (Bojak and Liley, 2005)

$$PSD\left(f\right) = 2\pi R^2 \int\_0^\infty \frac{dk}{k} \mathbf{l}\_1^2\left(kR\right) \left| h\_\epsilon\left(\omega \equiv 2\pi f, k\right) \right|^2,\qquad(11)$$

where J1 is a Besselfunction of the first kind. In practice we evaluate the integral Eq. 11 numerically using a 64 point 0 < *ki* < 14.14/cm Gauss–Legendre quadrature, and hence need to evaluate Eq. 10 for all these *k* = *ki* . PSDs calculated in this manner from 10 parameter sets selected out of the 73,454 in Bojak and Liley (2005) are shown in **Figure 2**.We call a solution *stable* if for all 64 *k* = *ki* , as well as for homogeneous cortex *k* = 0 cm, the eigenvalues are such that ∀*n*: <λ*n*(*k*) < 0. Only for stable cases do all the approximations leading up to Eq. 11 make sense. The largest contributions to Eq. 10 arise when ω = =λ*n*(*k*); and if one disregards the *cn*(*k*), then the "least stable" eigenvalue with largest <λ*n*(*k*) < 0 will contribute most.

Consider now only those λ*m* that have non-zero frequencies *fm* ≡ =λ*m*(*k* = 0/cm)/(2π) 6= 0: due to the selection process (Bojak and Liley, 2005), the "least stable" λmax of these λ*m* will have 8 Hz ≤ fmax ≤ 13 Hz, i.e., a frequency in the α rhythm range. If we change from parameter set {**P**1} to {**P**2}, we can compute the resulting frequency shift ∆*f* ≡ *f* {**P**2} max − *f* {**P**1} max of this eigenvalue α frequency. We find that for the parameter changes considered below, this "theoretical" α frequency shift ∆*f* estimated directly from the eigenvalues provides a reasonable approximation for a more"experiment-like"calculation of the α peak shift. In an experiment one would typically seek the maxima of the measured PSDs in the 8–13 Hz range, and then compute their difference in frequency in order to determine an α frequency shift, cf. **Figure 1C** (Hayashi et al., 2007; Tsuda et al., 2007).We can do something similar here by evaluating the full PSDs with Eq. 11 for {**P**1} and {**P**2}, respectively, and then compute the difference of the maxima of these theoretical predictions. However, we use the "theoretical" ∆*f*

**Table 1 |** α **peak frequency shifts predicted from the leading eigenvalue (**∆**f ) and the full PSD, respectively, for the parameter sets of Figure 2.**

P and K are normalized propofol and ketamine concentrations, respectively.

in the following. It is much easier to compute, since it involves only one eigendecomposition for *k* = 0/cm compared to 64 for *k* = *ki* needed in the numerical PSD integration Eq. 11. Furthermore, the "theoretical"∆*f* separates the change of the α peak frequency from other changes to the spectrum. "Experiment-like" calculations of shifts directly from local maxima in the full spectrum can be confounded easily by other spectral changes, and a prior subtraction of the spectral "background" around these maxima would closely match our "theoretical" procedure. For the parameter sets shown in **Figure 2**, a comparison between "theoretical" and "experimentlike" frequency shifts is provided by **Table 1**. How these shifts are generated will be discussed in the following, but note for now that most results are quite similar. The big discrepancies for parameter set IX are caused precisely by a rise of the spectral "background," as discussed.

#### **DRUG EFFECT PARAMETERIZATION AND SELECTION OF SETS**

The effect of the action of both ketamine and propofol on HCN1 channels is to hyperpolarize the resting membrane potentials of pyramidal (excitatory) cells (Chen et al., 2009). Consider Eq. 6 in the absence of synaptic inputs*Ilk*(**x**,*t*) ≡ 0, then lim*t*→∞*hk* (**x**, *t*) = *h r k* . Thus *h r e* and *h r i* parameterize the excitatory and inhibitory resting membrane potentials, respectively. In the spirit of Eq. 5 we hence use the following ansatz:

$$\begin{split} \Delta h\_{\epsilon}^{r} \equiv \left. h\_{\epsilon}^{r} \right|\_{P,K} - \left. h\_{\epsilon}^{r} \right|\_{P=K=0} = -\left( a\_{1}P + a\_{2}K + a\_{12}PK \right) \\ &= -\Delta h \cos \theta, \end{split} \tag{12}$$

$$\begin{aligned} \Delta h\_i^r \equiv \left. h\_i^r \right|\_{P,K} - \left. h\_i^r \right|\_{P=K=0} = -\left( b\_1 P + b\_2 K + b\_{12} PK \right) \\ = -\Delta h \sin \theta, \end{aligned} \tag{13}$$

where *P*, *K* are normalized (dimensionless) concentrations of propofol and ketamine, respectively; and ∆*h r e* , ∆*h r i* are changes of the excitatory and inhibitory resting membrane potentials, respectively, due to these drugs. For convenience we have factored out the sign corresponding to hyperpolarization, and we have assumed that inhibitory neurons would react qualitatively like the pyramidal cells, i.e., ∆*h r e* , ∆*h r i* ≤ 0 mV in the considered ranges 0 ≤ *P* ≤ *P*max and 0 ≤ *K* ≤ *K* max, while quantitative differences are expressed by potentially different coefficients. Since the drugs applied individually lead to hyperpolarization, we must have coefficients *a*1, *a*2, *b*1, *b*2 > 0, whereas the sign of the interaction coefficients *a*12, *b*12 carries the same meaning as that of η in Eq. 4. In the following it often will be useful to express the "Cartesian" ∆*h r e* , ∆*h r i* ≤ 0 mV in the corresponding "polar coordinate" form as ∆*h* ≡ q ∆*h r e* 2 + ∆*h r i* 2 ≥ 0 mV and θ ≡ arctan ∆*h r i* ∆*h r e* ∈ [0 ◦ , 90◦ ].

As a first step, we have investigated which of the 73,454 human α rhythm sets from Bojak and Liley (2005) can be extended viably via Eq. 12 and Eq. 13. Chen et al. (2009) found for rat pyramidal neurons that ∆*h r e* = −4.0 mV for ketamine at 20µM concentration and ∆*h r e* = −3.7 mV for propofol at 5µM concentration. Assuming that in humans (and in inhibitory neurons) hyperpolarizations of similar sizes occur, we varied both *h r e* and *h r i* away from their original values in steps of −0.05 mV up to a hyperpolarization of −6 mV, while the remaining parameters were left unchanged. This leads to a grid of 121 × 121 hyperpolarization combinations ∆*h r e* , ∆*h r i* , for which we tested whether the changed parameter sets remain stable, i.e., we computed eigendecompositions for 64 + 1 values of *k* and made sure that all eigenvalues had negative real parts. We also calculated the resulting shift in the α peak frequency as compared to the original parameter set in the "theoretical" manner discussed above: ∆*f* (∆*h r e* , ∆*h r i* ) ≡ *f* {**P**(∆*h r e* ,∆*h r i* )} max − *f* {**P**(0,0)} max .

We find that of the 73,454 parameter sets only 1,627 remain stable for all 121 × 121 combinations of hyperpolarizations up to −6 mV. This does not mean that the other parameter sets are thereby rejected on biological or physiological grounds; rather their PSDs cannot be calculated with the eigenspectrum approximation used here, but would have to be estimated from explicit simulations with the fully non-linear Eqs 6–8. This ordinary numerical procedure is several orders of magnitude slower and hence not employed here. **Figure 3A** displays the average ∆*f* (∆*h r e* , ∆*h r i* ) over the 1,627 stable sets. The color bar indicates the corresponding frequency values.We can see that in this average there is little effect of ∆*h r i* , whereas decreasing ∆*h r e* (increasing the hyperpolarization of the pyramidal neurons) leads to an increasingly negative h∆*f*i. The lowest average value for the 1,627 sets in **Figure 3A** is h∆*f*(−6 mV, −6 mV)i = −2.03 Hz, whereas the highest is h∆*f*(0 mV, −6 mV)i = 0.336 Hz.

Since more substantial increases in frequency are expected for the interaction of ketamine and propofol (Hayashi et al., 2007; Tsuda et al., 2007),we introduce the following cut: a set will be kept only if for at least one of the 121 × 121 hyperpolarization combinations ∆*h r e* , ∆*h r i* we find ∆*f* > 1.6 Hz. Similarly, since ketamine on its own should introduce a decrease in ∆*f* (Schuttler et al., 1987; Kochs et al., 1996), we require that for at least one other hyperpolarization combination ∆*f* < −0.8 Hz. Finally, propofol on its own is assumed here to not change the α frequency significantly ∆*f* ≈ 0 Hz (Schwender et al., 1996; Kuizenga et al., 1998, 2001; Feshchenko et al., 2004; Breshears et al., 2010; Cimenser et al., 2011), at least not by a HCN1-mediated mechanism, as was discussed in the Introduction. It is more difficult to introduce a simple cut for this property, since for small hyperpolarizations by definition one finds ∆*f* ≈ 0 Hz. We orient ourselves here to ∆*h r e* = −3.7 mV for propofol from Chen et al. (2009), and require that at least for one combination with ∆*h r i* ≤ −4.3 mV one has |∆*f*(−3.7 mV, ≤4.3 mV)| < 0.4 Hz. Considered individually, the low frequency cut for ketamine eliminates only 80 parameter sets, whereas the high frequency cut for the interaction of ketamine and propofol leaves only 66 parameter sets. Combining these two cuts then leaves 64 parameter sets in total. Individually, the cut for propofol limiting the frequency shift leaves 149 parameter sets. Combined with the other two cuts, we arrive at 10 parameter sets. Their original PSDs are the ones that were displayed previously in **Figure 2**, and we display their parameter values in **Table A1** in the Appendix. We show the resulting h∆*f*i, now averaging over only the 10 selected sets, in **Figure 3B**. It is immediately apparent that there are now three zones: for small ∆*h r e* but large (negative) ∆*h r i* one sees large increases in frequency,for large (negative) ∆*h r e* but small ∆*h r i* large decreases in frequency, and in between there is a corridor with little change in frequency. This same basic structure is found in all 10 selected sets individually.

Now we can use this structure to determine the coefficients in Eq. 12 and Eq. 13. Starting with the case of giving ketamine only, we can write

$$P = 0 : \tan \theta\_K = \frac{b\_2}{a\_2}. \tag{14}$$

Thus the effect of increasing ketamine concentration in the plane of hyperpolarizations is to move out along a line through the origin with angle θ*K*. Ketamine on its own is supposed to deliver shifts to low frequencies, for which we have set a cut ∆*f* < −0.8 Hz above. We now determine for every hyperpolarization combination how many of the 10 selected parameter sets have ∆*f* ∆*h r e* , ∆*h r i* < −0.8 Hz. This leads to a 121 × 121 grid of values between 0 and 10. In **Figure 3C** this is shown by blue contour lines for 4, 7, and 10 sets fulfilling this cut. We choose the mean of all ∆*h r e* , ∆*h r i* in the "maximal fulfillment" (10 sets) region (tip of blue arrow) to determine tan θ*K* = 0.315. Given that we do not know the dependence of hyperpolarizations on ketamine concentrations in humans, we choose ∆*h r e* ≡ −4 mV at *K* = 1 and thus consequently *a*2 ≡ 4 mV. This implies an unknown normalization *K* ≡ *cK* /*c* ∗ *K* , so that at a ketamine concentration *c* ∗ *K* one finds ∆*h r e* = −4 mV in humans. Given this choice, we have *b*2 ≡ 1.26 mV from the ketamine angle tan θ*K* = 0.315. Since *a*2 > *b*2, we can now also find *K* max = (6 mV)/*a*2 = 1.5 as the largest value for the normalized ketamine concentration for which both hyperpolarizations remain below −6 mV.

In a similar manner we can deal with the case of propofol as the sole drug. Then we find the angular dependence:

$$K = 0 : \tan \theta\_P = \frac{b\_1}{a\_1}.\tag{15}$$

**Figure 3C** shows contour lines for 4, 7, and 10 parameter sets fulfilling the cut for an α frequency shift ∆*f* ∆*h r e* , ∆*h r i* < 0.4 Hz, this time in green color. Since the cut was evaluated for ∆*h r e* = −3.7 mV only to find these sets, we determine the mean of combinations ∆*h r e* ≡ −3.7 mV, ∆*h r i* that have "maximal fulfillment" (10 sets) in order to obtain tan θ*P* = 1.297, indicated by the tip of the green arrow. We choose ∆*h r e* ≡ −3.7mV at *P*≡ *cP* /*c* ∗ *P* = 1, so that *a*1 ≡ 3.7 mV and at an unknown propofol concentration *c* ∗ *P* one finds ∆*h r e* = −3.7 mV in humans. Then *b*1 ≡ 4.8 mV, and since *b*1 > *a*1 it follows that *P*max = (6 mV)/ *b*1 = 1.25.

Finally, **Figure 3C** shows red contour lines for 4, 7, and 10 parameter sets fulfilling ∆*f* (∆*h r e* , ∆*h r i* ) >1.6 Hz. Again we find the mean of "maximal fulfillment" (10 parameter sets), as indicated by the tip of the red arrow. These mean values are ∆*h r e* ≡ −0.177 mV and ∆*h r i* ≡ −5.903 mV in this case. We now extend to the −6 mV hyperpolarization limit by setting ∆*h* max *i* ≡ −6 mV and ∆*h* max *e* ≡ (−6 mV)(−0.177 mV)/(−5.903 mV) = −0.180 mV.We can now solve the following two equations

$$
\Delta \overline{h}\_e^{\text{max}} = - \left( a\_1 P\_{\text{max}} + a\_2 K\_{\text{max}} + a\_{12} P\_{\text{max}} K\_{\text{max}} \right), \tag{16}
$$

$$
\Delta \overline{h}\_i^{\text{max}} = - \left( b\_1 P\_{\text{max}} + b\_2 K\_{\text{max}} + b\_{12} P\_{\text{max}} K\_{\text{max}} \right) . \tag{17}
$$

This will then mean that our entire hyperpolarization grid −6 mV ≤ ∆*h r e* , ∆*h r i* ≤ 0 mV will be projected onto a rectangular area bounded by 0 ≤ *P* ≤ *P*max and 0 ≤ *K* ≤ *K* max, respectively. Solving Eq. 16 and Eq. 17 with our previous results yields *a*12 = −5.57 mV and *b*12 = −1.01 mV. **Figure 3D** shows the projected h∆*f*(*P*, *K*)i. Clearly the intended α frequency shifts are now achieved: negative ones for only ketamine, none for only propofol, and positives ones for propofol and ketamine together.

#### **RESULTS**

We have parameterized the HCN1-mediated hyperpolarizations of neuron membrane potentials in order to reproduce

the observed EEG effects of ketamine and propofol, and in particular of their interaction when concurrent. The coefficients that we have obtained for Eqs 12–13 afford the following interpretation: pyramidal neurons react similarly to ketamine and propofol (*a*1 = 0.925 × *a*2), whereas inhibitory neurons react much more strongly to propofol than to ketamine (*b*1 = 3.81 × *b*2). Furthermore, and perhaps most interestingly, there is an antagonism of ketamine and propofol (*a*12, *b*12 < 0), which leads to infra-additivity in the investigated effect of HCN1-mediated hyperpolarization, cf. Eq. 4. This antagonism is stronger in pyramidal neurons *a*12/(*a*1*P*max + *a*2*K* max) = 4.10 *b*12/(*b*1*P*max + *b*2*K* max), though the precise proportion depends on the given concentrations of the drugs. Intuitively it makes sense however that in inhibitory neurons, where one drug is much more effective than the other, the antagonism between the drugs is less pronounced.

To illustrate these results we look again at the "theoretical" estimates of the α peak frequency in **Figure 4**, where we compare now the effects of changing propofol and ketamine concentration on the 10 selected sets (red) with those computed for all the valid 1,627 sets (gray). Note that the 1,627 sets include the 10 selected ones. Quantile bands are computed to summarize the results for the individual parameter sets, as indicated by the legend. Starting from a baseline without drugs, four phases are being considered: first, propofol concentration is increased linearly; then propofol is maintained at maximum concentration and ketamine concentration is increased linearly; next propofol concentration is decreased linearly while ketamine is maintained at maximum concentration, and finally ketamine concentration is decreased linearly for a return to the baseline. It should be noted that no attempt at modeling the pharmacodynamics/pharmacokinetics of ketamine and propofol drug action beyond drug interaction has been made here. Furthermore, the eigenspectrum approach assumes that the system has reached equilibrium for the given parameters. Thus every single *f*max(*P*, *K*) predicted here, and consequently every single quantile band value, represents a "steady state" result for that particular drug concentration combination. Hence one can for example view **Figure 4** from right to left, beginning with an increase in ketamine concentration, followed by an increase in ketamine at maximum propofol concentration, and so forth.

Comparing now the red with the gray quantile bands, we see that our cuts selected sets that react particularly dramatically to the concurrence of propofol and ketamine (phases 2 and 3), while

**FIGURE 4 | Estimated** α **peak frequency shifts.** Shifts of the α peak frequency for normalized propofol P and ketamine K concentrations estimated as described below Eq. 11, using the hyperpolarizations in Eq. 12 and Eq. 13. Either all 1,627 (gray) or the 10 selected sets (red) are used to compute quantile bands, as indicated by the legend. The median value is shown by a thick black or red line, respectively. There are four phases of drug variation, as indicated by the titles and dotted lines, quantified by bars below the main panel: first, P = 0→1.2 linearly, while K = 0. Then K = 0→1.4 linearly, while P = 1.2. Next P = 1.2→0, while K = 1.4. Finally, K = 1.4→0, while P = 0. No pharmacodynamics has been modeled here, so every (P, K) combination yields an independent "steady state" result. Hence for example an increase of P at high K is shown by the third phase viewed from right to left.

being unresponsive to propofol alone (phase 1). Nevertheless, it is not the case that the results for the 1,627 sets show a totally divergent response pattern. In fact, the median rise of estimated α peak frequency in phase 2 (upon introducing ketamine at maximum propofol concentration) is comparable: 1.88 Hz for the selected sets (from 11.48 to 13.36 Hz) vs. 1.56 Hz for all sets (from 9.37 to 10.93 Hz). Thus the predicted boost of α peak frequencies due to the interaction between ketamine and propofol is a robust result for all sets given our drug effect parameterization, which is infra-additive concerning HCN1-mediated hyperpolarization. The main difference appears to be rather that the α peak frequencies of the selected sets do not react significantly to propofol, whereas they are similar to all other sets in the reaction to ketamine and the interaction between these drugs.

Turning to results for full PSDs from Eq. 11, we will consider the 10 selected sets only due to the higher computational demands. **Figure 5** shows results for one individual set (Set III of **Figure 2**) under three variations of drug concentration. In **Figure 5A** we see that as desired and estimated, the α peak frequency stays roughly the same during propofol anesthesia (Schwender et al., 1996;Kuizenga et al., 1998, 2001; Feshchenko et al., 2004;Breshears et al., 2010; Cimenser et al., 2011). The damping seen here would be more characteristic of occipital than frontal regions, though other processes in particular related to the GABAA agonism could

modify these results. The "theoretical" ∆*f* = 0.24 Hz at *P* = 1.2 is larger than the "experiment-like" shift of 0.13 Hz. In **Figure 5B** we can see the reaction to increasing ketamine concentration. As expected, the α peak gets shifted to lower frequencies. The "theoretical" ∆*f* = −1.37 Hz estimate at *K* = 1.4 is somewhat lower than the "experiment-like" shift of the local maxima of the PSDs of −1.69 Hz. While the α oscillations get dampened, they contribute to a net increase in the θ frequency range (Schuttler et al., 1987; Kochs et al., 1996) thanks to their downward frequency shift. But one sees also a general rise in power at lower frequencies. **Figure 5C** shows that adding ketamine in the presence of high doses of propofol leads to a shift of the α peak into the beta band, as observed by Hayashi et al. (2007) and Tsuda et al. (2007). The "theoretical" ∆*f* = 1.91 Hz at *P* = 1.2 and *K* = 1.4 is smaller than the "experiment-like" shift of 2.27 Hz. We predict here an increase in power, unlike the experiment, which observed a small but significant reduction, and our frequency shift of 2.27 Hz is less than half the observed 4.7 Hz. But this could be explained easily by the missing GABAA and NMDA mechanisms, or the precise circumstances of the experiment. A bolus of ketamine was given by intravenous injection in the experiment, whereas here we calculate "steady state" results. Furthermore, while the frequency shifts from our "theoretical" ∆*f* method are largely in agreement with those obtained from the local α maxima of the PSDs, they do not

**FIGURE 5 | Power spectral densities for Set III under drug variations. (A)** PSDs for increasing normalized propofol concentration from none (thinnest green line) to 1.2 (thickest green line). **(B)** PSDs for increasing normalized ketamine concentration from none (thinnest blue line) to 1.4 (thickest blue line). **(C)** PSDs for increasing normalized ketamine concentration

from none (thinnest red line) to 1.4 (thickest red line), while normalized propofol concentration is held constant at 1.2. **(D)** Comparison of the PSDs representing the highest normalized concentrations from **(A)** in green, **(B)** in blue, and **(C)** in red. The black curve is the PSD without drugs. In all four panels the dotted line represents the position of the α peak of this curve.

agree perfectly – in spite of both being derived using the same eigenspectrum technique. The difference is that by looking at the maxima the results are influenced by changes to the overall spectrum, which provide the "background" on which the α resonance sits. **Table 1** gives ∆*f* and "experiment-like"α peak frequency shifts for the 10 sets at the three highlighted drug concentrations.

Finally, in **Figure 6** we show similar results for all the 10 selected sets. We follow here the same scheme of changing drug concentrations as in **Figure 4**. We see that the α peaks of the full PSDs (here shown in decibels by color) of the individual sets indeed follow the "zigzag" shape we saw in the quantile bands of **Figure 4**. Panel III in **Figure 6** can be directly compared to **Figure 5**, which we have just discussed. For example, **Figure 5A** corresponds to the first phase in panel III here. Overall we see that while the sets clearly change in a similar way, they all have individual features that set them apart from the others. For example, parameter set IV shows particularly strong changes in the low frequency range, whereas parameter set VII reacts with a particularly strong lowering of the α peak frequency in the presence of ketamine. These variations can be considered as representing the variations that one can also observe in humans.

## **DISCUSSION**

We have shown that observed changes of the EEG α peak frequency induced by the presence of the anesthetic agents propofol and ketamine, but in particular also by their interaction when given concurrently, can be explained based on the modeling of HCN1-mediated hyperpolarizations alone, at least qualitatively. This is surprising, since the main mechanism of action of these drugs is supposed to be through NMDA antagonism (ketamine) and GABAA agonism (propofol), respectively. However, since HCN1-knockout mice are indeed less sensitive to the hypnotic effects of both drugs (Chen et al., 2009), this would indicate that the EEG remains useful as an indicator of anesthetic action. It is perhaps interesting to note that while ketamine is famous for its hallucinatory action (Wolff and Winstock, 2006), and hence is considered a dissociative anesthetic agent, propofol as a classic general anesthetic agent is also capable of inducing a range of hallucinatory phenomena (Balasubramaniam and Park, 2003). Hence it is possible that psychotropic HCN1-mediated effects are simply masked behaviorally more by propofol's GABAA agonism than by ketamine's NMDA antagonism, but that the EEG is particularly sensitive to these underlying changes.

Only a fraction of all considered parameter sets (1,627 of the 73,454 parameter sets from Bojak and Liley (2005) proved "stable" under the HCN1-mediated hyperpolarization changes up to −6 mV on both excitatory and inhibitory neurons. However, this is at least partly due to the computational methods used here: the eigenspectrum method (Bojak and Liley, 2005) can only be used for "fixed point" dynamics. Rejected sets could follow the

**under drug variation.** We use here the same four phases of drug variation as in **Figure 4**, as indicated by the dotted lines and quantified by bars below the main panels: first, P = 0→1.2 linearly, while K = 0. Then K = 0→1.4 linearly, while P = 1.2. Next P = 1.2→0, while K = 1.4. Finally, K = 1.4→0, while P = 0. Every panel corresponds to 1 of the 10 selected parameter sets, as indicated by a white roman numeral. The PSD for one

line corresponding to frequencies from 0 to 20 Hz. Colors here indicate decibels of the PSD, with dark red corresponding to large, green to medium and dark blue to small values. (The "jet" colormap of Matlab has been mapped for each panel individually, to the full range of PSD decibel values shown in the panel.) A white dashed line indicates the α peak frequency in the absence of drugs.

same kind of drug-induced changes, but be inaccessible with our method chosen for its computational speed. In principle it would be desirable to carry out fully non-linear calculations instead, but under extra-cortical noise input it takes about a minute to estimate a sufficiently detailed and accurate power spectrum on a regular PC. For the 121 × 121 hyperpolarization changes of **Figure 3** that would lead to about 10 days of calculations even for a single parameter set – but we investigate here 73,454 different ones. The eigenspectrum method is several orders of magnitude faster even for a single neural mass. Yet we include here also effects due to the spatial distribution of neural masses, see Eq. 11. The computational load scales roughly linearly with the number of integration points for the eigenspectrum method, but roughly as a square for an equivalent fully non-linear simulation based on spatial grids. This increases the difference in computation speed even further. Thus for the investigations carried out here only the eigenspectrum method proves practicable. In addition, limit cycle or chaotic dynamics are often more representative of seizures or other pathological brain states, which in clinical practice would lead to the termination of the pharmacological intervention that we intend to describe here. Other than by numerical simulation in every single case, we do not know how to determine the characteristics of the dynamics past the point of instability.

Furthermore, our current investigation does not include the NMDA and GABAA actions commonly assumed to be dominant in these drugs. We speculate that a more complete simulation could allow the use of a larger fraction of the 73,454 parameter sets, since these omitted actions can affect the required stability. A prolongation of the inhibitory postsynaptic potentials due to GABAA, for example, could suppress excessive excitation and thus stabilize a parameter set. These neglected stabilizing effects would increase also in due proportion to the agent concentration, just as the potentially destabilizing hyperpolarizations we have modeled here do. In order to obtain spectral changes that demonstrate clearly the expected frequency shifts, we introduced three further selection cuts, leaving us with only 10 parameter sets out of the 1,627. Again we speculate that NMDA and GABAA actions may ameliorate this reduction. If this is not the case, then this may point to underlying correlations between neural parameters or functional properties that were not considered in Bojak and Liley (2005), but which now prove crucial for a realistic description. Note that in terms of the hyperpolarizations, cf. **Figures 3B,C**, the non-reactivity to propofol means that for the 10 selected parameter sets an increase in the hyperpolarization of excitatory neurons can be compensated by an increase in the hyperpolarization of inhibitory neurons. This could suggest a particular intrinsic balance of excitation and inhibition maintaining functional stability against extrinsic disturbances.

As is apparent from **Figure 2**, most of our 10 selected sets have relatively high α peak frequencies. However, this simply reflects an underlying bias in the original 73,454 parameter sets, cf. Figure 8 in Bojak and Liley (2005); and we were able to lower base α frequencies through the adjustment of a few parameters (e.g., cortico-cortical connectivity *N* β *lk* ) within physiological limits, without thereby qualitatively changing the relative frequency shifts due to propofol and ketamine. Other spectral features that could be selected for are also unlikely to affect the α frequency

shifts here qualitatively. The constraints used by Bojak and Liley (2005) leave plenty of room for such adjustments: effectively only four out of 14 system eigenvalues were used to establish the "1/*f*" background and the α resonance. Adding a weak beta frequency resonance to enhance realism, for example, would only require the adjustment of two further system eigenvalues. Thus we expect that our results here would hold true if one were to redo the parameter space search of Bojak and Liley (2005) first with additional constraints on the base power spectra. Finally, for concurrent application of propofol and ketamine we predict increases of the α peak frequency of around 2 Hz, falling short of the 4.7 Hz seen experimentally (Hayashi et al., 2007; Tsuda et al., 2007). Yet our calculations were for "steady state" concentrations, and for rapid increases in dosage, as in the experimental injection of a single ketamine bolus here, one often finds a more complex response due the pharmacodynamics and transient neural responses. The so-called biphasic responses to anesthetic agents (Kuizenga et al., 1998, 2001) has received theoretical attention from several groups (Bojak and Liley, 2005; Wilson et al., 2006; Molaee-Ardekani et al., 2007; Hutt and Longtin, 2010; Steyn-Ross and Steyn-Ross, 2010), see also the review in Foster et al. (2008). Hints of such a biphasic response could be visible in **Figure 5** of Tsuda et al. (2007), which shows a significant drop of the α peak shift with time from the initial 4.7 Hz to values around 2 Hz. However, this would have to be disentangled from the decrease in ketamine concentration due to natural clearance after the bolus. We intend to investigate all the mentioned issues in future work.

We found that we could account for the heterogeneous effects of ketamine on the EEG if we assumed that propofol and ketamine interacted in an *infra-additive or antagonistic* manner in their inhibition of HCN1-mediated neuronal membrane hyperpolarization. While most anesthetic and sedative agents are reported to interact synergistically ketamine is well-known to be a major exception (Hendrickx et al., 2008). The interaction between ketamine and GABAA agonists (most sedative/anesthetic agents) to produce hypnosis is reported to range from additivity to infraadditivity/antagonism (Hendrickx et al., 2008). On the basis of the limited clinical data available it appears that the GABAA agonist propofol can interact infra-additively to produce a given hypnotic endpoint. For instance Hui et al. (1995), in a study involving 180 female patients presenting for minor gynecological surgery, calculated quantal dose-response curves for propofol and ketamine administered alone and in combination. On the basis of logarithmic regression of a response surface model, it was first suggested that the dose-response for the combination was best explained by additivity, but on reanalysis (Hendrickx et al., 2008), by infraadditivity (significantly so for immobility, as trend for hypnosis). In apparent contradiction, Sakai et al. (1999) claim a significantly additive interaction between propofol and ketamine for the endpoints of unresponsiveness to vocal command and the loss of eyelash reflex. However, the former study determined the probability of a given hypnotic endpoint in response to *single doses* of fixed proportions of ketamine and propofol, whereas the latter gave *continuous infusions* of ketamine and propofol until a given hypnotic endpoint was reached. Another study by Frizelle et al. (1997) used smaller doses of ketamine and targeted a lower level of sedation (rousable to verbal stimuli) with a bolus of propofol and ketamine followed by concomitant infusion. They found no statistical evidence that the addition of ketamine reduced the required amount of propofol to reach their intended sedation level, suggesting once more infra-additivity. Clearly, additional work regarding the pharmacodynamic interactions of ketamine and propofol are required, as it is difficult to reconcile the results of these studies.

The use of neural field/mass approaches to modeling drug action on the EEG is emerging as a powerful explanatory framework (Liley and Bojak, 2005; Foster et al., 2008; Hutt, 2011), which is able to retain meaningful connections to the brain's physiology despite its mesoscopic scale of description. Because neural field/mass models generally have much smaller parameter and state spaces than biophysically plausible neural network

# **REFERENCES**

brain states under general anesthesia by using global coherence analysis. *Proc. Natl. Acad. Sci. U.S.A.* 108, 8832–8837.

models, they are not only much easier to parameterize and simulate, but offer a simplerframeworkfrom which to make predictions and derive hypotheses that can be empirically tested. Our ability to account qualitatively for the effects that propofol and ketamine have on the EEG adds to a growing list of phenomena that are amenable to neural field/mass description (Deco et al., 2008; Coombes, 2010; Bressloff, 2012; Liley et al., 2012). Ultimately it is hoped that by accounting for large scale phenomena using neural field/mass models, genuine and enduring insights into brain function and physiology will emerge.

## **ACKNOWLEDGMENTS**

Harry C. Day is supported by a scholarship from the College of Life and Environmental Sciences at the University of Birmingham.

neural populations.*Cogn. Neurodyn.* 4, 37–59.

treatment-resistant unipolar depression: current evidence. *CNS Drugs* 26, 189–204.

of ketamine on clinical endpoints of hypnosis and EEG variables during propofol infusion. *Acta Anaesthesiol. Scand.* 43, 212–216.

*Acta Anaesthesiol. Scand.* 51, 472–481.

**Conflict of Interest Statement:** David T. J. Liley is Chief Scientific Officer of Cortical Dynamics Ltd., an unlisted subsidiary of Biopharmica Ltd. (Perth, Australia), which is a medical device company focused on developing an EEG based depth of anesthesia monitor. David T. J. Liley is an inventor on several patent applications filed by Cortical Dynamics Ltd. since 2004 that describe new approaches to monitoring depth of anesthesia using the EEG. None of the IP declared in any published, pending, or granted patent has been licensed.

*Received: 07 December 2012; accepted: 10 March 2013; published online: 05 April 2013.*

*Citation: Bojak I, Day HC and Liley DTJ (2013) Ketamine, propofol, and the EEG: a neural field analysis of HCN1-mediated interactions. Front. Comput. Neurosci. 7:22. doi: 10.3389/fncom.2013.00022 Copyright © 2013 Bojak, Day and Liley. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which per-*

*mits 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.*

# **APPENDIX**

**Table A1 | The 10 selected parameter sets, whose PSDs are shown in Figure 2.**

# Neural mass modeling of power-line magnetic fields effects on brain activity

# *J. Modolo1,2,3\*, A. W. Thomas 1,2,3 and A. Legros 1,2,3,4*

*1 Human Threshold Research Group, Lawson Health Research Institute, London, ON, Canada*

*2 Department of Medical Biophysics, Western University, London, ON, Canada*

*3 Department of Medical Imaging, Western University, London, ON, Canada*

*4 School of Kinesiology, Western University, London, ON, Canada*

#### *Edited by:*

*Peter B. Graben, Humboldt-Universität zu Berlin, Germany*

#### *Reviewed by:*

*Fabrice Wendling, INSERM, France Flavio Frohlich, University of North Carolina, USA*

#### *\*Correspondence:*

*J. Modolo, Human Threshold Research Group, Lawson Health Research Institute, 268 Grosvenor Street, London, ON N6A4V2, Canada.*

*e-mail<EMAIL_ADDRESS>
Neural mass models are an appropriate framework to study brain activity, combining a high degree of biological realism while being mathematically tractable. These models have been used, with a certain success, to simulate brain electric (electroencephalography, EEG) and metabolic (functional magnetic resonance imaging, fMRI) activity. However, concrete applications of neural mass models have remained limited to date. Motivated by experimental results obtained in humans, we propose in this paper a neural mass model designed to study the interaction between power-line magnetic fields (MFs) (60 Hz in North America) and brain activity. The model includes pyramidal cells; dendrite-projecting, slow GABAergic neurons; soma-projecting, fast GABAergic neurons; and glutamatergic interneurons. A simple phenomenological model of interaction between the induced electric field and neuron membranes is also considered, along with a model of post-synaptic calcium concentration and associated changes in synaptic weights Simulated EEG signals are produced in a simple protocol, both in the absence and presence of a 60 Hz MF. These results are discussed based on results obtained previously in humans. Notably, results highlight that (1) EEG alpha (8–12 Hz) power can be modulated by weak membrane depolarizations induced by the exposure; (2) the level of input noise has a significant impact on EEG power modulation; and (3) the threshold value in MF flux density resulting in a significant effect on the EEG depends on the type of neuronal populations modulated by the MF exposure. Results obtained from the model shed new light on the effects of power-line MFs on brain activity, and will provide guidance in future human experiments. This may represent a valuable contribution to international regulation agencies setting guidelines on MF values to which the general public and workers can be exposed.

**Keywords: neural mass models, power-line magnetic fields, electroencephalogram (EEG), synaptic plasticity, brain stimulation**

# **INTRODUCTION**

Since the pioneering work of Wilson and Cowan (1973) and Amari (1977), neural field models have been increasingly used, expanded and studied by a developing multidisciplinary community including: mathematicians, physicists, neuroscientists, medical imaging scientists *etcetera*. Neural field models provide a concise, yet insightful description of cortical activity. This theory has not only led to successful reproduction of numerous experimental results, but also to the prediction of a certain number of phenomena that have been observed *in vivo* (for a review of these phenomena, see Modolo et al., 2010). A popular simplification of neural field models consists in neglecting the role of space, consider neural populations present in cortical columns (such as pyramidal neurons) and to consider a connectivity matrix between the different neural masses considered (Wendling et al., 2002; Sotero and Trujillo-Barreto, 2008; Bojak et al., 2010). This approach has been used to build large-scale models of brain activity, including models simulating the electroencephalogram (EEG), or the blood oxygen level-dependent signal (BOLD) reflecting metabolic activity of brain tissue (Wendling et al., 2002; Sotero and Trujillo-Barreto, 2008; Bojak et al., 2010). Indeed, since neural mass models are tractable and make the link between local variables and observables, these models appear as an excellent comprise between biological realism and computational complexity. Therefore, here we present an application of neural mass model to a specific question that many teams over the world have been tackling for a number of years using various neuroimaging modalities: how do power-line frequency (60 Hz in North America) magnetic fields (MFs) interact with human brain activity (interaction mechanisms), and how does it translate into observable outcomes (neuroimaging data such as EEG, motor/cognitive performance)?

The effects of extremely low-frequency (categorized as being <300 Hz) MF such as power-line MF on human neurophysiology have been studied for several decades. Despite an impressive amount of experimental data *in vitro*, *in vivo*, and in humans, a complete understanding of the interaction mechanisms and associated effects is still to be achieved. One complication in comparing outcomes from these studies is the wide range of MF flux densities, MF exposure setups, and exposure protocols used (Crasson, 2003). First, significant work has been done has been done *in vivo* and *in vitro* in order to characterize the effects of electric fields in terms of membrane potential perturbation, excitability and neural network oscillations (Jefferys et al., 2003; Bikson et al., 2004; Deans et al., 2007; Fröhlich and McCormick, 2010; Reato et al., 2010). These studies have highlighted the role of neuronal morphology and orientation were critical in understanding the interaction with electric field, but also that membrane depolarization far below the firing threshold can influence the activity of neuronal networks. Second, among the reported effects in humans, let us mention the modulation of pain threshold (Ghione et al., 2005), effects on resting tremor (Legros and Beuter, 2005), modulations in functional brain activity as measured by BOLD (Legros et al., 2011), interference with learning in a short-term memory test (Corbacio et al., 2011). The most established interaction mechanism of MF exposure consists in induced currents, resulting from Faraday's law of induction, stating that a time-varying MF induces a time-varying electric field. This electric field will induce charge movement, creating a current. This is termed as the *induced current mechanism* (National Institute of Environmental Health Sciences of the National Institutes of Health, 1998). The induced current mechanism will be the mechanism considered in this paper. Interestingly, several studies have reported lasting effects associated with ELF MF exposure, i.e., an effect that is still detectable after cessation of the exposure. Using a specific pulsed MF, modulations of the EEG have been observed post-exposure (Cook et al., 2004, 2005). In the case of 60 Hz MF in the millitesla range (1.8 and 3 mT), modulations of the BOLD signal measured in humans during motor (finger tapping) and cognitive (mental rotation) tasks (Legros et al., 2010; Miller et al., 2010) have been found to be modulated post-exposure.

The interest for the interaction between power-line MF and human neurophysiology is twofold. First, there is a growing concern from the general public regarding the possible deleterious effects of power-line MF on human health, even if such negative effects remain to be demonstrated. Second, international regulation agencies such as ICNIRP (International Commission on Non-Ionizing Radiation Protection, http://www.icnirp.de) need results from the scientific literature to set their exposure guidelines (ICNIRP, 2010), aiming to protect the general public and workers, that are used by governments. One of the most long awaited data is the threshold in MF flux density at 60 Hz resulting in detectable effects in humans. One well-known effect of ELF MF exposure is magnetophosphenes, the perception of flickering lights in the visual field in the presence of a sufficiently strong MF. ICNIRP states that: "*Since the perception of magnetophosphenes constitutes the most reliable effect of MF exposure on human biology, this serves as a basis for the ICNIRP guidelines*" (ICNIRP, 2010). However, no data is available at 60 Hz regarding threshold values, therefore new approaches that could assist in the interpretation of existing experimental results, but also in the prediction of new results such as an estimation of threshold values that could be tested experimentally, would constitute significant advances in the field.

In order to shed light on the mechanisms involved in lasting effects of 60 Hz MF exposure on human neurophysiology, and also provide an estimate of the threshold value resulting in detectable changes in EEG caused by 60 Hz MF exposure, we present in this paper a neural mass model aiming to model brain tissue dynamics at different time scales, bridging biophysical mechanisms with changes in observables. Among the panel of brain tissue dynamics models available, the approach initiated by Jansen and Rit (1995), consisting in considering sub-populations of neurons synaptically connected; later extended and improved by Wendling et al. (2002), provides a meaningful and accurate description of cortical dynamics. Indeed, such models have been successfully applied to understand the transition between baseline EEG and epileptic activity (Wendling et al., 2002, 2005; Molaee-Ardekani et al., 2010). First, we present an extension of the neural mass model developed by Sotero and Trujillo-Barreto (2008), by including a population of fast GABAergic neurons as suggested by Wendling et al. (2002). Second, we include a simple biophysical model of interaction between the 60 Hz MF and neuronal activity, along with a model of synaptic plasticity changes related to postsynaptic calcium concentration levels (Shouval et al., 2002a,b). Third, we use this model to investigate *in silico* the effects of 60 Hz MF on cortical dynamics, notably to evaluate the threshold in MF flux density resulting in detectable changes in variables of interest. The role of synaptic input noise, neuronal populations modulated by the exposure, and synaptic plasticity are also explored. Finally, we discuss future directions of research using this modeling approach.

# **MATERIALS AND METHODS**

# **NEURAL MASS MODEL**

In order to develop a biologically grounded model to study the effect of 60 Hz MF on neuronal activity, we have extend the thalamo-cortical model proposed by Sotero and Trujillo-Barreto (2008) by including a population of soma-projecting, fast inhibitory γ-amino-butyric acid (GABA) interneurons to extend the possible dynamical repertoire of the model [as shown by Wendling et al. (2002)]. The proposed modification of the block diagram proposed by Sotero and Trujillo-Barreto (2008) used to describe the thalamo-cortical model is the following:

Model equations are obtained by using the fact that the synaptic response function (Green's function) for a type of synapse *i* (e.g., glutamatergic) writes as *Vi*(*t*) = *Ai*.*ai*.*t*. exp(−*ai*.*t*), where *Ai* is the response amplitude and *ai* the response time constant. Considering the temporal operator *L*ˆ = *d*2 *dt*2 + 2*ai d dt* + *a*2 *i* , and using the fact that the synaptic response is a Greens' function for the temporal operator, we can use *LV*ˆ *i*(*t*) = δ(*t*) to write the following neural mass equation:

$$\frac{d^2}{dt^2}V\_i(t) + 2.a. \frac{d}{dt}V\_i(t) + a^2. V\_i(t) = A.a. \upsilon\_i(t) \tag{1}$$

where υ*i*(*t*)is the incoming firing rate. It is often practical to write Equation (1) under the form of a system of two first-order differential equations:

$$\begin{aligned} \frac{d}{dt}V\_i(t) &= \wp\_i(t) \\ \frac{d}{dt}\wp\_i(t) &= A.a.\upsilon\_i(t) - 2.a.\frac{d}{dt}\wp\_i(t) + a^2.V\_i(t) \end{aligned} \tag{2}$$

**permission).** The inclusion of the new population of fast GABAergic neurons and its connectivity with other neuronal populations is highlighted in red.

Using this principle for each block of EPSP/IPSP presented in **Figure 1**, it is possible to formulate our extended thalamo-cortical model as a system of 22 differential equations (6 new equations corresponding to the new population of fast inhibitory interneurons, its feedback loop with pyramidal neurons, and its inhibitory input from slow inhibitory interneurons) presented below:

$$\begin{aligned} \dot{\mathcal{Y}}\_1^{\vec{n}\dagger}(t) &= \mathcal{Y}\_{12}^{\vec{n}\dagger}(t) \\ \dot{\mathcal{Y}}\_2^{\vec{n}\dagger}(t) &= \mathcal{Y}\_{13}^{\vec{n}\dagger}(t) \\ \dot{\mathcal{Y}}\_3^{\vec{n}\dagger}(t) &= \mathcal{Y}\_{14}^{\vec{n}\dagger}(t) \\ \dot{\mathcal{Y}}\_4^{\vec{n}\dagger}(t) &= \mathcal{Y}\_{15}^{\vec{n}\dagger}(t) \\ \dot{\mathcal{Y}}\_5^{\vec{n}\dagger}(t) &= \mathcal{Y}\_{16}^{\vec{n}\dagger}(t) \\ \dot{\mathcal{Y}}\_6^{\vec{n}\dagger}(t) &= \mathcal{Y}\_{17}^{\vec{n}\dagger}(t) \\ \dot{\mathcal{Y}}\_8^{\vec{n}\dagger}(t) &= \mathcal{Y}\_{19}^{\vec{n}\dagger}(t) \\ \dot{\mathcal{Y}}\_9^{\vec{n}\dagger}(t) &= \mathcal{Y}\_{20}^{\vec{n}\dagger}(t) \end{aligned}$$

$$\begin{aligned} \dot{\mathcal{Y}}\_{11}^{j}(t) &= \mathcal{Y}\_{21}^{j}(t) \\ \dot{\mathcal{Y}}\_{11}^{j}(t) &= \mathcal{Y}\_{22}^{j}(t) \\ \dot{\mathcal{Y}}\_{12}^{j}(t) &= A.a. \left\{ c\_{\mathcal{S}} \mathcal{S} \left[ \boldsymbol{\gamma}\_{1}^{\mathcal{N}}(t) - \boldsymbol{\gamma}\_{2}^{\mathcal{N}}(t) - \boldsymbol{\gamma}\_{3}^{\mathcal{N}}(t) \right] \\ &+ c\_{\mathcal{S}} \mathcal{S} \left[ \boldsymbol{\gamma}\_{3}^{\mathcal{N}}(t) \right] + K^{d,n} c\_{\mathcal{S}}^{i} \mathcal{S} \left[ \boldsymbol{\chi}\_{4}(t) \right] \right\} \\ &+ A.a. \left\{ \sum\_{m=1}^{M\_{\pi}} \left( k\_{c1}^{mj} c\_{\mathcal{T}} \mathcal{S} \left[ \boldsymbol{\gamma}\_{8}^{\mathcal{N}m}(t) \right] + k\_{c2}^{mj} c\_{\mathcal{S}} \mathcal{S} \left[ \boldsymbol{\gamma}\_{7}^{\mathcal{N}m}(t) \right] \right) \\ &\quad m = 1 \\ &+ \sum\_{i=1}^{N} K^{i,n} \sum\_{m=1}^{M\_{i}} c\_{\mathcal{S}} \mathcal{S} \left[ \boldsymbol{\gamma}\_{5}^{\mathcal{N}m}(t) \right] \\ &+ \sum\_{i=1}^{N} \left( \boldsymbol{\gamma}\_{i}^{\mathcal{N}}(t) - a^{2} \boldsymbol{\gamma}\_{1}^{\mathcal{N}}(t) + A.a. p^{\mathcal{T}}(t) \right) \\ &- 2.a. \mathcal{N}\_{12}^{m}(t) - a^{2} \boldsymbol{\gamma}\_{1}^{\mathcal{N}}(t) + A.a. p^{\mathcal{T}}(t) \end{aligned}$$

*y*˙ *nj* 13(*t*) = *B*.*b*. ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ *c*4.*S y nj* 4 (*t*) + *Mn m* = 1 *m* -= *j k mj i* .*c*9.*S ymn* 8 (*t*) ⎫ ⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎭ −2.*b*.*y nj* 13(*t*) − *b*2*y nj* 2 (*t*) *y*˙ *nj* 14(*t*) = *A*.*a*. - *c*1.*S y nj* 1 (*t*) − *y nj* 2 (*t*) − *y nj* 9 (*t*) + *Kth*,*n*.*c t* 4*S* [*x*5(*t*)] −2.*a*.*y nj* 14(*t*) − *a*2*y nj* 3 (*t*) *y*˙ *nj* 15(*t*) = *A*.*a*. - *c*3.*S y nj* 1 (*t*) − *y nj* 2 (*t*) − *y nj* 9 (*t*) + *Kth*,*n*.*c t* 5*S* [*x*6(*t*)] −2.*a*.*y nj* 15(*t*) − *a*2*y nj* 4 (*t*) *y*˙ *nj* 16(*t*) = *A*.*ad*1. - *S y nj* 1 (*t*) − *y nj* 2 (*t*) − *y nj* 9 (*t*) −2.*ad*1.*y nj* 16(*t*) − *a*2 *d*1*y nj* 5 (*t*) *y*˙ *nj* 17(*t*) = *A*.*ad*2. - *S y nj* 1 (*t*) − *y nj* 2 (*t*) − *y nj* 9 (*t*) −2.*ad*2.*y nj* 17(*t*) − *a*2 *d*2*y nj* 6 (*t*) *y*˙ *nj* 18(*t*) = *A*.*ad*3. - *S y nj* 3 (*t*) − 2.*ad*3.*y nj* 18(*t*) − *a*2 *d*3*y nj* 7 (*t*) *y*˙ *nj* 19(*t*) = *B*.*bd*4. - *S y nj* 4 (*t*) − 2.*bd*4.*y nj* 19(*t*) − *b*2 *d*4*y nj* 8 (*t*) *y*˙ *nj* 20(*t*) = *G*.*g*. - *c*12.*S y nj* 10(*t*) − *y nj* 11(*t*) + *Kth*,*n*.*c t* 5*S* [*x*6(*t*)] −2.*g*.*y nj* 20(*t*) − *g*2 *y nj* 9 (*t*) *y*˙ *nj* 21(*t*) = *A*.*a*. - *c*10.*S y nj* 1 (*t*) − *y nj* 2 (*t*) − *y nj* 9 (*t*) −2.*a*.*y nj* 21(*t*) − *a*2*y nj* 10(*t*) *y*˙ *nj* 22(*t*) = *B*.*b*. - *c*11.*S y nj* 4 (*t*) − 2.*b*.*y nj* 22(*t*) − *b*2*y nj* 11(*t*)

The equations describing the activity of the thalamus, composed of a population of thalamocortical cells and a population of reticular cells are [modified from Sotero and Trujillo-Barreto (2008)]:

$$\begin{aligned} \dot{\mathbf{x}}\_1(t) &= \mathbf{x}\_7(t) \\ \dot{\mathbf{x}}\_2(t) &= \mathbf{x}\_8(t) \\ \dot{\mathbf{x}}\_3(t) &= \mathbf{x}\_9(t) \\ \dot{\mathbf{x}}\_4(t) &= \mathbf{x}\_{10}(t) \\ \dot{\mathbf{x}}\_5(t) &= \mathbf{x}\_{11}(t) \\ \dot{\mathbf{x}}\_6(t) &= \mathbf{x}\_{12}(t) \\ \dot{\mathbf{x}}\_7(t) &= A\_{t}.a\_{t^\*} \left\{ \sum\_{i=1}^N K^{th,i} \sum\_{m=1}^{M\_i} c\_6.\mathcal{S} \left[ \boldsymbol{\nu}\_5^{im}(t) + p\_{lh}(t) \right] \right\} \\ &- 2.a\_{t^\*} \mathbf{x}\_7(t) - a\_t^2 \mathbf{x}\_1(t) \end{aligned}$$

$$
\dot{\mathbf{x}}\_8(t) = B\_t.b\_t.c\_{2t}.S\left[c\_{1t}\mathbf{x}\_3(t)\right] - 2.b\_t.\mathbf{x}\_8(t) - b\_t^2\mathbf{x}\_2(t)
$$

$$
\dot{\mathbf{x}}\_9(t) = A\_t.a\_t.S\left[\mathbf{x}\_1(t) - \mathbf{x}\_2(t)\right] - 2.a\_t.\mathbf{x}\_9(t) - a\_t^2\mathbf{x}\_3(t)
$$

$$
\dot{\mathbf{x}}\_{10}(t) = A\_t.a\_{d1t}.S\left[\mathbf{x}\_1(t) - \mathbf{x}\_2(t)\right] - 2.a\_{d1t}.\mathbf{x}\_{10}(t) - a\_{d1t}^2\mathbf{x}\_4(t)
$$

$$
\dot{\mathbf{x}}\_{11}(t) = A\_t.a\_{d2t}.S\left[\mathbf{x}\_1(t) - \mathbf{x}\_2(t)\right] - 2.a\_{d2t}.\mathbf{x}\_{11}(t) - a\_{d2t}^2\mathbf{x}\_5(t)
$$

$$
\dot{\mathbf{x}}\_{12}(t) = A\_t.a\_{d3t}.S\left[\mathbf{x}\_1(t) - \mathbf{x}\_2(t)\right] - 2.a\_{d3t}.\mathbf{x}\_{12}(t) - a\_{d3t}^2\mathbf{x}\_6(t)
$$

The physical meaning and values of model parameters are detailed in **Table A1** (Appendix section). For more details, the reader can refer to Sotero and Trujillo-Barreto (2008). Overall, the model is composed of 22 differential equations describing cortical dynamics of four different neuronal populations (pyramidal neurons, glutamatergic interneurons, fast/slow GABAergic neurons, and 12 differential equations defining thalamic activity). Therefore, this set of 34 differential equations describes the thalamocortical activity including time delays between cortical areas, connectivity parameters, synaptic responses derived from neurophysiology, and a biologically plausible (even if it is obviously simplified) circuitry between the neuronal populations considered. **Table A1** provided in Appendix summarizes the parameters used in the model, with new parameters added due to the population of fast inhibitory interneurons that have been highlighted.

## **MODEL OF INTERACTION BETWEEN 60 Hz EXPOSURE AND NEURON MEMBRANES**

In order to model the interaction between the electric field induced by 60 Hz MF exposure and neural tissue, we have used the λ · *E* model in order to simulate the modulation of neuron membrane polarization (Molaee-Ardekani et al., 2013). In this model, the membrane depolarization *dV* in the presence of an electric field **E** is a function of a constant **λ** termed "polarization length" (Radman et al., 2009). More precisely, the membrane depolarization is expressed as *dV =***λ** · *E*, where **λ** is a vector oriented along the neuron fibre, and *E* is the electric field vector. This expression is valid for a static electric field. In the case of a timevarying electric field (such as the electric field induced by 60 Hz MF), a frequency-dependent term needs to be included (Gianni et al., 2006), resulting in:

$$dV = \frac{\lambda \cdot E}{\sqrt{1 + \omega^2 \mathbf{r}^2}} \tag{5}$$

where ω = *2*π*f* , *f* being the frequency (in our case, *f* = 60 Hz), and τ is the polarization time constant. In order to use Equation (5), describing membrane depolarization induced by the induced electric field, at the level of a neuronal population, we made the following assumptions: (1) the induced 60 Hz electric field is homogeneous in space at the level of the neural mass (i.e., the MF flux density is constant at each point of the cortical column); (2) the MF-induced membrane depolarization is applied to pyramidal neurons only because of their large size compared to other types of neurons in the human cortex; (3) pyramidal neurons in a given neural mass all have the same spatial orientation. Taken

together, assumptions (1) and (3) result in identical *dV* values for all pyramidal neurons in a neural mass at a given time. These assumptions lead to the use Equation (3) in the context of a neural mass model.

This was achieved by modifying the expression of the total post-synaptic potential at the level of pyramidal neurons:

$$
\ddot{y}\_1(t) - \ddot{y}\_2(t) - \ddot{y}\_9(t) \to \ddot{y}\_1(t) - \ddot{y}\_2(t) - \ddot{y}\_9(t) + dV(t) \tag{6}
$$

As it is commonly calculated in neural mass model, the EEG signal was computed as the summation of excitatory and inhibitory post-synaptic potential at the level of pyramidal neurons:

$$EEG(t) = \wp\_1(t) - \wp\_2(t) - \wp\_9(t) + dV(t) \tag{7}$$

Let us mention that, even if the MF-induced depolarization is included in the model as a simple additive perturbation, it has the potential to induce non-linear effects. Indeed, the effective potential at the level of pyramidal neurons (7) is used as an input for other neuronal populations, and is transformed from a potential to a firing rate using a sigmoid function, which is fundamentally non-linear. In the "Results" section, we have used arbitrary values for the field-induced membrane depolarization *dV*, guided by preliminary simulation results. Based on *dV* values resulting in significant changes in the EEG alpha power with or without synaptic plasticity in the model, we will provide an estimate of the corresponding level of 60 Hz MF flux density. This will provide us with an order of magnitude of the 60 Hz MF flux density threshold value that should result in effects detectable experimentally in humans.

In our simulations, we focused specifically on the EEG alpha rhythm (8–12 Hz). The reason of this choice is twofold. First, as mentioned in the Introduction, there is converging evidence that extremely low-frequency MF in the millitesla range, such as 60 Hz MF, can induce EEG alpha activity modulation. Second, the model we have developed is basically an extension of the Jansen and Rit model, designed to model EEG alpha activity. It is possible to reproduce other types of EEG rhythms (e.g., beta −13 to 30 Hz), for example by introducing heterogeneity in the time constant of neural populations over different neural masses (Wendling et al., 2002), which exceeds the scope of this paper.

## **BIOPHYSICAL MODEL OF SYNAPTIC PLASTICITY**

In order to investigate the hypothesis that 60 Hz MF exposure might modulate with human neurophysiology by modulating synaptic plasticity, we have implemented a simplified model of synaptic plasticity, based on the biophysical model developed by Shouval et al. (2002a,b). It is now well accepted that the mechanisms of long-term synaptic potentiation and depression (LTP/LTD, respectively) involve changes in post-synaptic calcium concentration and the trafficking of α-amino-3-hydroxy-5-methyl-4-isoxazole-propionic acid (AMPA) glutamate receptors between the intracellular medium and the synapse site. Depending on the calcium concentration, AMPA receptors can either insert into the membrane at the level of the synaptic cleft, or undergo an endocytosis, which is termed *receptor trafficking* (Collingridge et al., 2004). An increase in the number of AMPA receptors at the synaptic level will increase the number of glutamate molecules that can bind on post-synaptic membranes, thereby increasing membrane depolarization during a synaptic event. Consequently, the number of post-synaptic AMPA receptors is directly proportional to the synaptic weight. The model proposed by Shouval et al. (2002a,b) has been a significant progress in the modeling of the biophysical processes at play during LTP/LTD. This model is based on the "calcium control hypothesis," according to which the level of post-synaptic calcium is the main factor regulating the exocytosis/endocytosis rate of AMPA receptors, and therefore the dynamics of synaptic plasticity changes.

We have adapted the model by Shouval et al. (2002a,b) to our neural mass model, and despite some simplifications with respects to the original model; our synaptic plasticity model captures some of its essential features. Based on the experimental literature on 60 Hz MF exposure effects on the EEG alpha rhythm, we assume that (1) no qualitative changes of EEG dynamics will occur due to 60 Hz MF exposure, changes will be purely quantitative (i.e., EEG alpha rhythm amplitude/spectral power changes, but no qualitative change in dynamical regime such as a transition toward high-amplitude, low-frequency spiking); (2) the coupling between the synaptic plasticity model and the 60 Hz MF is *via* the equation linking the EEG with the post-synaptic calcium concentration, occurring on long timescales (depending on the opening of N-methyl-D-aspartate (NMDA) glutamate receptors, not represented in the model). Therefore, the model offers the possibility to test the hypothesis that 60 Hz MF exposure can modulate synaptic plasticity by interfering with the calcium fluxes at the level of synapses. However, it does not take into account possible effects of 60 Hz MF exposure on spike timing (see the "Discussion" section).

Let us consider the average calcium post-synaptic concentration in a neural mass. The model proposed by Shouval et al. (2002a,b) links the calcium current at the level of NMDA receptors with the calcium concentration, and finally to a differential equation describing the dynamics of synaptic weight change as a function of two different calcium-dependent functions. The time constant of calcium concentration dynamics is long (on the order of minutes), and the calcium concentration increases with the membrane potential. Therefore, it appears reasonable to approximate calcium dynamics by a low-pass filtering of the mean potential of a given neural mass:

$$\text{tr}\_{\text{Ca}^{2+}} \frac{d}{dt} [\text{Ca}^{2+}] + [\text{Ca}^{2+}] = \chi (\wp\_1 - \wp\_2 - \wp\_9) \tag{8}$$

where *y*1 − *y*2 − *y*9 is the "EEG" signal at the level of a neural mass (e.g., summation of post-synaptic potentials at the level of pyramidal neurons as defined previously). Once the post-synaptic dendritic calcium concentration is obtained, it is possible to evaluate the calcium-dependent functions η and present in the Shouval et al. (2002a,b) model, used to express the dynamics of the synaptic weight *ci* (*i* denoting the type of synapse in the neural mass, e.g., afferent glutamatergic synapses on pyramidal neurons) at the level of a given neural mass:

$$\frac{d\mathbf{c}\_i(t)}{dt} = \eta(t) \cdot \left[\Omega\_{Ca^{2+}}(t) - \mathbf{c}\_i(t)\right] \tag{9}$$

The function was approximated by a combination of piecewiselinear and quadratic functions (see Appendix, **Figure A1** for details) similar to the function proposed in Shouval et al. (2002b). The function used in our model differs quantitatively from the one proposed in Shouval et al. (2002a,b), since the authors were linking with this function the level of post-synaptic calcium concentration with the relative change in synaptic weight, where we directly link the post-synaptic calcium concentration with the synaptic weight itself. Nevertheless, the used in this paper captures the most important qualitative properties proposed by Shouval et al. (2002a,b). In our simulations, we have assumed that was η(*t*) a constant, [*Ca*2+(*t*)] being bounded between 0 and 1μM. We assumed that the synapses modulated were the synapses terminating on pyramidal neurons, pooled in the constant c5, becoming the variable *c*5(*t*) in our model. Numerical implementation for the neural mass model was performed using Matlab 2010 (The Mathworks, USA) on a quad-core Apple iMac (2.66 GHz/CPU) with 8 GB of RAM. The simulation of a neural mass using the complete model during 2 h with a time step of *dt* = 1 ms took typically 8 min.

## **60 Hz MF EXPOSURE PROTOCOL**

In order to study the effects of 60 Hz MF on the simulated EEG, we used the following protocol: the neural mass was simulated during 2 h overall with a 1 ms resolution, which was decomposed as (1) 30 min without 60 Hz MF exposure (termed "sham," of sufficient duration to reach a steady state); (2) 60 min with 60 Hz MF exposure (sufficient to reach the new steady state); and (3) 30 min without 60 Hz MF exposure. Previous research in our team using fMRI to image the functional changes in brain activity due to 60 Hz MF exposure involved comparable durations (notably, a 60 min exposure period and fMRI acquisitions performed before and after, see Legros et al., 2010). Simulations were performed both (1) using the synaptic plasticity model, and (2) using a fixed synaptic weight value taken as the steady state value when synaptic plasticity was taken into account. By doing so, we aimed at decomposing the respective contribution of the 60 Hz sinusoidal perturbation in membrane potential one the one hand, and of possible calcium-related synaptic plasticity modulations on the other hand.

# **RESULTS**

# **EFFECT OF 60 Hz MF ON THE EEG IN THE MODEL**

As an example, we present in **Figure 2** an example of simulated EEG data, and associated mean post-synaptic calcium concentration and synaptic weight obtained by solving Equations (3, 4, 8, 9).

We have investigated the effect of increasing values (125, 250, 500, and 1000μV) for the MF-induced membrane depolarization on the EEG alpha power. EEG alpha spectral power was computed before, during and after the 1-h 60 Hz MF exposure. 10 runs of 7200 s were performed for each tested value of *dV*. The averaged EEG alpha power for each condition (before, during, and after exposure) is presented in **Figure 3**.

From the results presented in **Figure 3**, it appears that increasing values of *dV* gradually decreases EEG alpha power during

exposure. In order to test the significance of the amplitude of *dV* on the EEG alpha power during 60 Hz MF exposure, we conducted a statistical analysis of the results. We performed a 4 × 3 × 2 ANOVA for repeated measures (SPSS 21, IBM, USA), respectively testing for the effects of "*dV*" (125, 250, 500, and 1000μV), "time" (before, during and after), and "plasticity" (with/without synaptic plasticity). The standard *p*-value of 0.05 (Greenhouse-Geisser) was chosen as the threshold for significance, and *p*-values were corrected for multiple comparisons. The statistical results reveal a significant decrease of the EEG alpha power for *dV* = 500μV as compared to the other values of *dV* (*p* < 0.001). This indicates that the threshold for a significant decrease of EEG alpha power due to 60 Hz MF exposure lies between induced membrane depolarization values of 250–500μV. In the next section, we attempt to link the membrane depolarization values to the corresponding MF flux density at 60 Hz.

Due to the possibility that the weak 60 Hz membrane depolarization can be seen as an additive noise, we have tested the influence of the input noise level [*p*(*t*) in the model, see **Table A1** of the Appendix] variance on EEG alpha power modulation due to the exposure. Since the dynamics of the model itself depend critically on the input noise level, we have indeed investigated the possibility that the 60 Hz MF exposure has an effect of variable amplitude depending on input noise. The interest is that model predictions could be tested experimentally (e.g., in an experimental setting where different levels of visual input would be tested). Therefore, the objective was not to study the influence of the noise level on the model dynamics, but rather how the effects of the 60 Hz MF on model dynamics are dependent on the noise level. Four different values of noise variance were tested (σ = 120, 150, 180, and 210 spikes/s) for the same maximal *dV* value of 0.5 mV. 10 runs of 7200 s following the same protocol than previously were run for each noise level value (40 runs total). The influence of the input noise level on EEG alpha power modulation by the 60 Hz MF is presented in **Figure 4**.

The results presented in **Figure 4** highlight the importance of the input noise level of the model. If the input *p*(*t*), representing external noisy input to the neural mass, is too high; then the effect induced by the 60 Hz MF on EEG alpha power modulation decreases. This has an immediate consequence on threshold values of MF flux density resulting in detectable effects in brain activity: the MF flux density needed to elicit a response in brain tissue will be lower in the presence of a low level of noise. Interestingly, there is experimental evidence that the visual input can play a role on the effects of MF exposure in humans, with an higher effect when the eyes are closed (Legros et al., 2011). EEG alpha oscillations increase dramatically eyes closed, and decrease in the presence of a visual input, that increases the input noise to the occipital cortex. Therefore, even if there is a considerable gap between the model and human data, it is tempting to make a parallel between smaller effects of 60 Hz MF in the model with high levels of noise, and smaller effects of 60 Hz MF exposure eyes open with an increased input noise level. One advantage of using our neural mass over interpreting experimental results is the possibility to point at precise mechanisms by which the observed decrease in EEG alpha activity occurs due to the 60 Hz MF exposure. From

a physiological point of view, it is relevant to investigate which neuronal pathways are mainly modulated by the exposure. In the model, we observe an immediate decrease in the activity of the loop between pyramidal neurons and slow GABAergic neurons, likely increasing the effect of excitatory input. To complement this observation, it is relevant to note that, using a bifurcation theory analysis of the Jansen and Rit model (the core of our model), Grimbert and Faugeras (2006) have shown that increasing the input noise level at the level of pyramidal neurons in the alpha oscillations regime (corresponding to a Hopf bifurcation) had the effect to decrease the amplitude of alpha oscillations. Therefore, the 60 Hz MF stimulus used in our model seems to have a similar effect than an additive, positive constant membrane depolarization on pyramidal neurons. This results physiologically speaking from an efficiency decrease of the slow inhibitory GABAergic feedback at the pyramidal neurons level. In terms of dynamical systems theory, this seems to be the natural result of increased input level in a specific dynamical system on a Hopf cycle.

the weak 60 Hz membrane potential perturbation becomes "buried" in

In order to distinguish between the contribution of the MF-induced membrane polarization on the one hand, and changes in synaptic plasticity on the other hand, we ran the same simulations than previously for four different values of *dV* (125, 250, 500, and 1000μV), with a constant value for the synaptic weight *c*5. The objective was to identify if synaptic plasticity was affecting the direction (increase/decrease of EEG alpha spectral power) or amplitude of the effects. In the following, the value *c*5 of was chosen as the steady-state value in the case where synaptic plasticity was considered. 10 simulations of 7200 s were ran for each value of *dV*. The results are presented in **Figure 5**.

noise.

From the results in **Figure 5**, it appears that the modulation of post-synaptic calcium concentration and corresponding changes in synaptic weight plays a minimal role in EEG alpha power modulation due to the 60 Hz MF exposure, and does not impact qualitatively the result (the direction of the effects is the same, and the amplitude of the effects is minimally affected). Indeed, the results from the ANOVA shows no significant interaction effect between synaptic plasticity and *dV* values (*p* = 0.253). This indicates that the presence of the synaptic plasticity mechanisms included in the model does not significantly change the effect of the membrane depolarization. It seems however to induce a non-significant increase the amplitude of the 60 Hz MF exposure effect. Different choices for the function linking the post-synaptic calcium concentration level with the updated synaptic weight - lead to similar results (not shown). Therefore, it seems that, if a modulation of synaptic plasticity explains lasting effects of 60 Hz MF exposure, it does not occur primarily by the modulation of post-synaptic calcium currents. However, it is still plausible that receptor trafficking and synaptic plasticity could be impacted by a perturbation of spike timing due to the 60 Hz MF exposure, a mechanism not included in the present model, which we discuss later.

In the model develop by Molaee-Ardekani et al. (2013), investigating the effects of transcranial direct current stimulation (tDCS), the neurons being modulated by the induced field were pyramidal neurons and inhibitory interneurons. Since the simulated EEG could result in different outcomes due to the exposure depending on the neuronal populations simulated, possibly assisting in discriminating between different interaction mechanisms; we have also simulated the EEG in the case of a 60 Hz MF exposure modulating the activity of different populations of inhibitory interneurons. To do so, in a similar fashion to Equation (7) describing the membrane depolarization of pyramidal neurons, we simulated different scenarios: (1) slow and fast GABAergic interneurons are involved, (2) slow GABAergic interneurons are involved, and (3) fast inhibitory neurons are involved. Depending on the scenario, we also added the variable *dV*(*t*) to *y*4(*t*) (slow GABAergic interneurons), or to *y*10(*t*)-*y*11(*t*) (fast GABAergic interneurons), and to both of these quantities for scenario (1). We have simulated the EEG for a similar protocol than previously (30 min without exposure, 1 h of exposure, 30 min without exposure), with a maximal value of *dV* = 1000μV. The resulting EEG alpha power before, during and after exposure for each of these scenarios is presented in **Figure 6**.

value of the MF-induced membrane depolarization was *dV* = 1 mV. Pyramidal neurons are considered to be modulated in each scenario. Pyr, pyramidal neurons only; slow+fast, pyramidal neurons, slow and fast GABAergic interneurons; slow, pyramidal neurons, slow GABAergic neurons; fast, pyramidal neurons, fast GABAergic neurons.

The results presented in **Figure 6** show that, if the fast GABAergic interneurons are modulated by the 60 Hz MF in addition to pyramidal neurons, the difference in the EEG is minimal. However, if the slow GABAergic interneurons are modulated, the decrease in EEG alpha power is much smaller, dropping from 17% (pyramidal neurons only) to 5%. Therefore, if the slow GABAergic interneurons are also modulated by the induced electric field due to the 60 Hz MF exposure, the threshold leading to a systematic EEG alpha power modulation will be higher. From the interpretation on the decrease in EEG alpha power when the pyramidal neurons alone are stimulated, where a decrease in the efficiency of the loop between pyramidal neurons and slow inhibitory interneurons is observed, we can speculate that the effect of modulating the slow GABAergic interneurons in addition to pyramidal neurons has an opposite effect of increasing the activity of this loop. Consequently, if slow GABAergic interneurons are modulated by the exposure to the 60 Hz MF exposure, the model suggests that it would result in a compensation of the modulation of pyramidal neurons' activity alone, thereby increasing the threshold in MF flux density leading to a systematic decrease in alpha. In other words, the modulation of pyramidal neurons and slow GABAergic interneurons activity would have competing effects regarding the decrease in EEG alpha power.

# **THRESHOLD OF 60 Hz MF FLUX DENSITY RESULTING IN DETECTABLE EFFECTS**