3. WC choose sigmoid response functions such that *SE*(0) = *SI*(0) = 0 in order to force (0,0) to be a critical (equilibrium) point where the variable time derivatives equal zero. Such critical points may be either state or unstable. If a critical point is stable, any brain dynamic state coming sufficiently close to this point will become fixed (forever static or "brain dead"). Of more interest to us are unstable critical points associated with on-going oscillations (limit cycles), possibly underlying EEG. In WC, the tissue response functions can become negative, a physiological impossibility; thus,in our following modified analysis, the WC conditions are replaced by

$$\mathcal{S}\_E[P(t), \left< E(t) \right>, \left] \ge 0$$

$$\mathcal{S}\_I[\left< E(t) \right>, \left] \ge 0 \tag{A3}$$

4. WC interpret the WC parameter µ as a single membrane time constant. Based on the classic solution of the cable equation and the distribution of excitatory synapses on dendrites with inhibitory synapses typically near cell bodies and our increased appreciation of wide ranges of excitatory and inhibitory time constants τ*E*, τ*I* , we consider casesτ*E* 6= τ*I* .

### **A MODIFIED VERSION OF THE WC ANALYSIS**

Define the non-dimensional time *t*1 = *t* τ*E* and time constant ratio *A* = τ*I* τ*E* so that with rescaled variables Eq. A1 become

$$\frac{d\langle E\rangle}{dt\_1} = -\left\langle E\right\rangle + \left(1 - \left\langle E\right\rangle\right) S\_E \left(P + \left.\omega\_{EE}\left\langle E \right\rangle - \left.\omega\_{IE}\left\langle I\right\rangle\right)\right.\tag{1}$$

$$A\frac{d\left\langle I\right\rangle}{dt\_1} = -\left\langle I\right\rangle + \left(1 - \left\langle I\right\rangle\right) S\_I \left(\left.\omega\_{EI}\left\langle E \right\rangle - \left.\omega\_{II}\left\langle I\right\rangle\right)\right.\tag{A4}$$

Here the non-dimensional parameter *A* may range from somewhat less than one to as high as perhaps 5. For convenience we drop the subscript 1 on the non-dimensional time variable. We introduce the new dependent variables*XE*(*t*),*XI*(*t*),which provide perturbations about the critical point (*E*0, *I* 0)

$$
\begin{aligned}
\langle E(t) \rangle &= E\_0 + X\_E(t) \\
\langle I(t) \rangle &= I\_0 + X\_I(t)
\end{aligned}
\tag{A5}
$$

We assume that the input function *P*(*t*) is exclusively excitatory such that *P*(*t*) ≥ 0 at all times. The following conditions follow from the variable definitions

$$0 \le E\_0 + X\_E(t) \le 1$$

$$0 \le I\_0 + X\_I(t) \le 1\tag{A6}$$

Equations (A4) then yield

$$\frac{dX\_E}{dt} = -\left[E\_0 - X\_E + (1 - E\_0 - X\_E)\,\mathrm{S}\_E\,\left[P + \omega\_{EE}(E\_0 + X\_E)\right]\right]$$

$$-\omega\_{IE}(I\_0 + X\_I)]$$

$$A\frac{dX\_I}{dt} = -\left[I\_0 - X\_I + (1 - I\_0 - X\_I)\,\mathrm{S}\_I\,\left[\omega\_{EI}(E\_0 + X\_E)\right]\right]$$

$$-\omega\_{II}(I\_0 + X\_I)\,\mathrm{I}\,\tag{A7}$$

The sigmoid functions *SE*, *SI* are chosen here such that (*E*0, *I* 0) is a critical point when *P*(*t*) = 0. Thus, we choose the following alternate sigmoid tissue response functions.

$$S\_I = \frac{1}{1 + \exp\left[-\mathsf{w}\underline{\mu}(E\mathsf{o} + X\_E) + \mathsf{w}\underline{\mu}(\mathsf{lo} + X\mathsf{i}) + K\mathsf{i}\right]}\tag{A8}$$

$$S\_E = \frac{1}{1 + \exp\left[-\omega\_{EE}(E\_0 + X\_E) + \omega\_{ie}(I\_0 + X\_I) + K\_E\right]} \quad \text{(A9)}$$

The constants *KE*, *KI* are added to the original WC sigmoid functions, thereby determining the response range during oscillations about fixed points. This choice of the forms of the sigmoid response functions insures that (1) (*E*0, *I* 0) or (*X,Y* ) = (0,0) is a critical point and (2) 0 ≤ *SE* ≤ 1 and 0 ≤ *SI* ≤ 1. Substitution of Eqs A8 and A9 into Eq. A7 and setting *XE*(*t*) = *XI*(*t*) = 0 yields the sigmoid constants in terms of the critical point

$$K\_I = \left. \boldsymbol{w}\_{EI} \boldsymbol{E}\_0 - \left. \boldsymbol{w}\_{II} \boldsymbol{I}\_0 + \text{Log} \left[ \frac{1}{\boldsymbol{I}\_0} - 2 \right] \right|\_{\text{m}} \qquad (\text{A10})$$

$$K\_E = \left. \varkappa\_{EE} E\_0 - \left. \varkappa\_{IE} I\_0 + \log \left[ \frac{1}{E\_0} - 2 \right] \right| \tag{A11}$$

Substitution of Eqs A10 and A11 into Eqs A8 and A9 yields

$$\begin{aligned} S\_I \left[ \left. \boldsymbol{w}\_{EI} (\boldsymbol{E}\_0 + \boldsymbol{X}\_E) - \boldsymbol{w}\_{II} (\boldsymbol{I}\_0 + \boldsymbol{X}\_I) \right] \right] \\ &= \frac{1}{1 + \left( \frac{1}{I\_0} - 2 \right) \exp \left( -\boldsymbol{w}\_{EI} \boldsymbol{X}\_E + \boldsymbol{w}\_{II} \boldsymbol{X}\_I \right)} \quad I\_0 < \frac{1}{2} \end{aligned} \tag{A12}$$

$$\begin{aligned} \text{SE}\left[\mathsf{w}\_{\text{EE}}(\text{Eo} + X\_{\text{E}}) - \mathsf{w}\_{\text{IE}}(\text{I}\_{0} + X\_{\text{I}})\right] \\ = \frac{1}{1 + \left(\frac{1}{E\_{0}} - 2\right) \exp\left(-\mathsf{w}\_{\text{EE}}X\_{\text{E}} + \mathsf{w}\_{\text{IE}}X\_{\text{I}}\right)} \quad E\_{0} < \frac{1}{2} \quad \text{(A13)} \end{aligned} $$

#### **PHASE-PLANE ANALYSIS**

We assume *wII* ∼= 0 in all of the following analyses based on our preliminary studies: In simulations with non-zero *wII* the effect is only to set the activity level of the system and has no significant influence on the dynamics. The first step in the analysis of Eq. A7 is to find the nature of the critical point (0,0). To accomplish this we expand the functions *F*(*XE*, *XI*) and *G*(*XE*, *XI*) about (0,0), where these functions are the expressions on the right sides of Eq. A7, that is

$$F(X\_E, X\_I) = -E\_0 - X\_E + (1 - E\_0 - X\_E) \text{ S}\_\varepsilon$$

$$G(X\_E, X\_I) = [-I\_0 - X\_I + (1 - I\_0 - X\_I) \text{ S}\_i]/A \tag{A14}$$

Taylor expansion about the critical point (0,0) yields equations of the general form

$$\begin{split} \frac{dX\_E}{dt} &= F(X\_E, X\_I) \cong F(0, 0) + \left(\frac{\partial F}{\partial X\_E}\right)\_0 X\_E + \left(\frac{\partial F}{\partial X\_I}\right)\_0 X\_I \\ &= aX\_E + bX\_I \end{split} $$

$$\begin{split} \frac{dX\_I}{dt} &= G(X\_E, X\_I) \cong G(0, 0) + \left(\frac{\partial G}{\partial X\_E}\right)\_0 X\_E + \left(\frac{\partial G}{\partial X\_I}\right)\_0 X\_I \\ &= cX\_E + dX\_I \end{split} \tag{A15}$$

Here we have forced *F*(0,0) = *G*(0,0) = 0 by proper choice of the constants *KE*, *KI* in Eqs A10 and A11. The partial derivatives are evaluated at (0,0) yielding the parameters (*a*, *b*, *c*, *d*). Eq. A15 then consist of two first order linear equations governing the dynamic behavior of the non-linear system close to (0,0). Define the parameters β = *a* + *d* and γ = *ad* − *bc*; the eigenvalues λ1,2 of the linear system satisfy λ 2 − βλ + γ =0 with solution

$$
\lambda\_{1,2} = \frac{1}{2} \left( \mathfrak{k} \pm \sqrt{\mathfrak{k}^2 - 4\gamma} \right) \tag{A16}
$$

We are mainly interested in oscillatory solutions about the critical point (0,0); that is, stable limit cycle solutions. These are expected when (0,0) is an *unstable* spiral, which occurs when β > 0 and β 2 < 4γ. By contrast, stable spirals result in damped oscillations. Saddle points and nodes result in non-oscillatory solutions (stable or unstable) that are of minimal interest here.

Consider the following example with *E*0 = *I* 0. The critical point (*E*0, *I* 0) is unstable if the following condition is met

$$\omega\_{EE} > \frac{A+1}{AE\_0 \left(1 - 2E\_0\right)} \quad 0 < E\_0 < \frac{1}{2} \tag{A17}$$

Note that *A* = τ*I* τ*E* so if inhibitory time constants are much shorter than excitatory time constants, larger values of *wEE* are required to produce linear instability. For the physiologically interesting range 0.1 < *E*0 < 0.4 and *A* = 1, all *wEE* > 25 cause the fixed point to be unstable.

For the case *E*0 = *I* 0 and *A* = 1, the fixed point is a spiral if

$$
\omega\_{EE} < 2\sqrt{\mathbf{w}\_{IE}\mathbf{w}\_{EI}}\tag{A18}
$$

By combining Eqs A17 and A18 we find the necessary (but possibly not sufficient) conditions for a stable limit cycle about (*E*0, *E*0) when *A* = 1, that is

$$\frac{2}{E\_0 \left(1 - 2E\_0\right)} < \omega\_{EE} < 2\sqrt{\omega\_{IE}\omega\_{EI}}\tag{A19}$$

If this condition is met, the corresponding spiral frequency is

$$
\omega\_{spiral} = N(E\_0)\sqrt{4\omega\_{IE}\omega\_{EI} - \omega\_{EE}^2} \tag{A20}
$$

Here the numerical factor lies in the range 0.00790 < *N*(*E*0) < 0.0294 if 0.1 < *E*0 < 0.4. An *unstable* spiral point at (*E*0, *I* 0) suggests a likely *stable limit cycle*, but the limit cycle frequency will not generally equal the (linear) spiral frequency. If *wEE* exceeds the upper limit in Eq. A19 an unstable node or saddle point will occur. In this case the solutions *XE*(*t*), *XI*(*t*) are likely to grow beyond physiologically realistic ranges, implying that the basic WC equations are no longer valid.

# Cortical information flow in Parkinson's disease: a composite network/field model

# *Cliff C. Kerr 1,2,3\*, Sacha J. Van Albada4, Samuel A. Neymotin1,5, George L. Chadderdon1, P. A. Robinson2,3 and William W. Lytton1,6*

*1 Department of Physiology and Pharmacology, State University of New York Downstate Medical Center, Brooklyn, NY, USA*

*2 School of Physics, University of Sydney, NSW, Australia*

*3 Brain Dynamics Centre, Westmead Millennium Institute, Westmead, NSW, Australia*

*4 Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6), Jülich Research Centre and JARA, Jülich, Germany*

*5 Department of Neurobiology, Yale University, New Haven, CT, USA*

*6 Department of Neurology, Kings County Hospital, Brooklyn, NY, USA*

#### *Edited by:*

*Dimitris Pinotsis, University College London, UK*

#### *Reviewed by:*

*Petra Ritter, Charité University Medicine Berlin, Germany Stephan Van Gils, Universiteit Twente, Netherlands*

#### *\*Correspondence:*

*Cliff C. Kerr, Department of Physiology and Pharmacology, State University of New York Downstate Medical Center, 450 Clarkson Avenue, Po Box 31, Brooklyn, NY 11203, USA. e-mail<EMAIL_ADDRESS>
The basal ganglia play a crucial role in the execution of movements, as demonstrated by the severe motor deficits that accompany Parkinson's disease (PD). Since motor commands originate in the cortex, an important question is how the basal ganglia influence cortical information flow, and how this influence becomes pathological in PD. To explore this, we developed a composite neuronal network/neural field model. The network model consisted of 4950 spiking neurons, divided into 15 excitatory and inhibitory cell populations in the thalamus and cortex. The field model consisted of the cortex, thalamus, striatum, subthalamic nucleus, and globus pallidus. Both models have been separately validated in previous work. Three field models were used: one with basal ganglia parameters based on data from healthy individuals, one based on data from individuals with PD, and one purely thalamocortical model. Spikes generated by these field models were then used to drive the network model. Compared to the network driven by the healthy model, the PD-driven network had lower firing rates, a shift in spectral power toward lower frequencies, and higher probability of bursting; each of these findings is consistent with empirical data on PD. In the healthy model, we found strong Granger causality between cortical layers in the beta and low gamma frequency bands, but this causality was largely absent in the PD model. In particular, the reduction in Granger causality from the main "input" layer of the cortex (layer 4) to the main "output" layer (layer 5) was pronounced. This may account for symptoms of PD that seem to reflect deficits in information flow, such as bradykinesia. In general, these results demonstrate that the brain's large-scale oscillatory environment, represented here by the field model, strongly influences the information processing that occurs within its subnetworks. Hence, it may be preferable to drive spiking network models with physiologically realistic inputs rather than pure white noise.

**Keywords: neural field model, spiking neural networks, Parkinsons's disease, thalamus, cortex, basal ganglia, Granger causality, interlaminar processing**

# **1. INTRODUCTION**

Parkinson's disease (PD) is a multiscale phenomenon, encompassing pathology at the level of single neurons, local networks, large neuronal ganglia, and the complex interactions between these ganglia and the cortex. PD is caused by the degeneration of dopaminergic neurons in the substantia nigra pars compacta, with the damage later spreading to dopaminergic neurons in the ventral tegmental area (Cools, 2006). The loss of dopaminergic input alters the dynamics of the striatum, which then affects the dynamics of large portions of the thalamus and cortex, which in turn affects the spinal cord and muscles (Bolam et al., 2002). Striatal dynamics are crucial to several large-scale projection pathways, including the well-characterized direct and indirect pathways. Dopaminergic input to the striatum increases transmission in D1-expressing striatal neurons involved in the direct pathway. These neurons inhibit the globus pallidus internal segment (GPi). Dopaminergic input also decreases input to D2-expressing striatal neurons involved in the indirect pathway. These neurons inhibit the globus pallidus external segment (GPe), which in turn inhibits the GPi. Thus, alterations to the direct and indirect pathways in PD are both thought to increase the firing rate of the GPi, which in turn inhibits the thalamus. There is also a hyperdirect pathway from the cortex to the GPi via the subthalamic nucleus (STN), as well as other lesser pathways (**Figure 1**).

Numerous models of PD and the basal ganglia have been proposed, using either field or network approaches. Van Albada and Robinson (2009) developed a field-based model of the basal

ganglia/thalamocortical system. This model was shown to reproduce realistic firing rates of each neuronal population in both healthy and PD states. One early network model was that of Terman et al. (2002), which represented a small network of neurons in the GPe and STN. A considerably larger and more complex (non-spiking) network model was developed by Leblois et al. (2006). This model explored both basal ganglia and thalamocortical cell populations, looking at competition between the direct and hyperdirect pathways. They suggested that PD disrupted this competitive balance, resulting in loss of the network's ability to select motor programs. Another network model focusing on motor-selection abilities was developed by Humphries et al. (2006), who also found that decreased dopamine interfered with the basal ganglia's capacity for selecting actions. Network models have also been used to analyze and predict the effects of deep brain stimulation on basal ganglia nuclei (Hahn and McIntyre, 2010; Guo and Rubin, 2011; Dovzhenok et al., 2013).

Previous neuronal network models of PD have either not included a cortex at all (Terman et al., 2002; Rubchinsky et al., 2003; Park et al., 2011), approximated it as a random Poisson process (Humphries et al., 2006), or considered it as a single layer with a single cell type (Leblois et al., 2006). The thalamus has also either been omitted or treated as a single population. In this work, we sought to fill this gap by exploring the interactions of the largescale dynamics of basal ganglia, represented by a field model, with a far smaller but more spatially detailed network model of the thalamus and six-layered cortex.

# **1.1. COMPOSITE MODEL**

The primary aim of this paper is to determine how the largescale dynamics of the brain affect the information flow in small networks of neurons. Most previous brain modeling efforts have been directed at one of these two scales, rather than their interaction. These efforts have consisted of either (1) neural field models that describe the dynamics of the whole brain, without explicitly modeling the activity of individual neurons (Nunez, 1974; Jirsa and Haken, 1996; Robinson et al., 1997; Destexhe and Sejnowski, 2009), or (2) spiking neuronal network models that capture individual neurons' dynamics, but are many orders of magnitude smaller than the brains of even the simplest vertebrates (Lumer et al., 1997; Neymotin et al., 2011b). Several large network models have also been published that have roughly as many "neurons" as the full mammalian brain (Izhikevich and Edelman, 2008; Ananthanarayanan et al., 2009). However, these models have not yet reproduced large-scale dynamics with the same degree of fidelity as neural field models. For example, the model of Izhikevich and Edelman (2008) showed simultaneous peaks in the delta and alpha bands, whereas experimentally these peaks are characteristic of sleep and wakefulness, respectively, and are hence rarely observed simultaneously (Niedermeyer and Lopes da Silva, 1999). Such infidelity may be because the enormous computational resources required to run these models makes it impractical to constrain their parameters by fitting their dynamics to experimental data.

Recently, both Deco and Jirsa (2012) and Wilson et al. (2012) described approximations that allow small and large spatial scales to be spanned at a mesoscopic level of description, allowing largescale dynamics (e.g., BOLD signals) to be related to small-scale network properties (e.g., criticality). Robinson and Kim (2012) took a different approach: they described the theoretical basis of combining spiking network and neural field components into a single model. The fundamental challenge in combining these two modeling approaches is to create a common representation of neuronal activity, since individual spikes are used in network models, while field models use average firing rates. Converting individual spikes into an average firing rate is a straightforward reduction of dimensionality: one simply needs to average over multiple neurons in the model. In contrast, converting an average population firing rate into individual spikes in multiple neurons requires an increase in dimensionality. This is a degenerate problem, so additional assumptions must be made. One approach, described in Robinson and Kim (2012), is to treat each neuron as a phase oscillator. The average firing rate then represents the instantaneous rate of phase change, with a given neuron firing whenever its phase advances by 2π radians. However, here we used an alternative approach, in which the average firing rate is taken as the instantaneous rate for an ensemble of Poisson processes. These are then used to generate individual spike times (Dayan and Abbott, 2001; Leblois et al., 2006; Chadderdon et al., 2012). This approach produces variability in spike timings even with a constant average firing rate, as is seen in real neuronal populations.

## **2. METHODS**

The model we used consisted of a network of spiking neurons that was "embedded" in a neural field model. The embedding consisted of having the field model generate spikes (via an ensemble of Poisson processes) that were used to drive the network model. Except where otherwise noted, all analyses were performed on the network model. The complete model is publicly available via ModelDB: https://senselab.med.yale.edu/modeldb/ ShowModel.asp?model=147366.

#### **2.1. NEURAL FIELD MODEL**

The neural field model was based on the work of Van Albada and Robinson (2009) and Van Albada et al. (2009). The neuronal populations and connections that constitute this model are shown in **Figures 1** and **3A** respectively. The basal ganglia nuclei modeled were the striatum, internal and external pallidal segments, and STN. The internal pallidal population can be thought of as including the substantia nigra pars reticulata, which has very similar connections and properties. The substantia nigra pars compacta was not explicitly modeled, except through its effects on the other nuclei. The thalamus was modeled as two populations: the inhibitory thalamic reticular nucleus (TRN) and the excitatory thalamocortical relay nuclei (TCR). The cortex was also modeled as two populations, representing inhibitory interneurons and excitatory pyramidal neurons. Since together these neuronal populations comprise a large portion of the brain, a network formulation would be computationally intractable. Except for a unitless normalization constant, all parameter values were based on anatomical and physiological data, as listed in Table 2 of Van Albada and Robinson (2009).

In neural field models, neuronal properties are spatially averaged. The dynamics are then governed by a set of equations relating the mean firing rates of populations of neurons to changes in mean cell-body potential, which are in turn triggered by mean rates of incoming spikes. The neural field model used here was based on a previously published model of the electrophysiology of the thalamocortical system (Robinson et al., 1997, 2001, 2002, 2005; Rennie et al., 1999), which in turn was based on earlier field models (Wilson and Cowan, 1973; Nunez, 1974; Freeman, 1975; Steriade et al., 1990; Wright and Liley, 1996).

The first component of the model is the description of the average response of populations of neurons to changes in mean cell-body potential. The mean firing rate *Qa* of each population *a* is the maximum attainable firing rate *Q*max *a* times the proportion of neurons with a membrane potential *Va* above the mean threshold potential θ*a*. This can be approximated by the sigmoid function

$$Q\_a(\mathbf{r}, t) = \frac{Q\_a^{\text{max}}}{1 + e^{-\{V\_a(\mathbf{r}, t) - \theta\_a\}/\sigma'}},\tag{1}$$

where **r** is the spatial coordinate, *t* is time, and σ is √3/π times the standard deviation of the distribution of firing thresholds (Wright and Liley, 1995). This function increases smoothly from 0 to *Q*max *a* as *Va* changes from −∞ to ∞.

The change in the mean cell-body potential due to afferent activity depends on the mean number of synapses *Nab* from neurons of population *b* to neurons of population *a* (note that the direction of projection *b* → *a* follows the conventions of control theory and matrix multiplication). The change in potential also depends on *sab*, the time-integrated change in cell-body potential per incoming spike. Defining ν*ab* = *Nabsab*, the change in the mean cell-body potential in neurons of population *a* is (Robinson et al., 2004).

$$D\_{a\emptyset}(t)V\_a(t) = \sum\_b \upsilon\_{ab} \phi\_b(t - \tau\_{ab}),\tag{2}$$

$$D\_{\alpha\beta}(t) = \frac{1}{\alpha\beta} \frac{d^2}{dt^2} + \left(\frac{1}{\alpha} + \frac{1}{\beta}\right) \frac{d}{dt} + 1. \tag{3}$$

Here, φ*b*(*t* − τ*ab*) is the incoming firing rate, τ*ab* represents the axonal time delay for signals traveling from population *b* to population *a* neurons, and α and β are the decay and rise rates of mean cell-body potential. The differential operator *D*αβ(*t*) represents dendritic and synaptic integration of incoming signals (Robinson et al., 1997; Rennie et al., 2000). The synapses and dendrites form an effective low-pass filter with a cut-off frequency between 1/α and 1/β.

In this model, neuronal activity spreads along the cortex in a wavelike fashion. This reflects previous models (Nunez, 1995; Jirsa and Haken, 1996; Bressloff, 2001) as well as experimental observations of such waves following cortical stimulation (Burns, 1951; Nunez, 1974; Rubino et al., 2006). Estimates of characteristic axonal ranges and propagation speeds suggest that these waves are significantly damped on the scale of the human cortex (Robinson et al., 2001, 2004; Wright and Liley, 1995). Assuming that the range distribution of corticocortical fibers decays exponentially at large distances, activity propagates according to a 2D damped-wave equation of the form (Robinson et al., 1997)

$$Q\_a(\mathbf{r}, t) = \left[ \frac{1}{\chi\_a^2} \frac{\partial^2}{\partial t^2} + \frac{2}{\chi\_a} \frac{\partial}{\partial t} + 1 - r\_a^2 \nabla^2 \right] \phi\_a(\mathbf{r}, t), \qquad (4)$$

where γ*a* = *va*/*ra* is the damping rate, consisting of the average axonal transmission speed *va* (10 m·s −1) and the characteristic axonal range *ra*. In practice, most types of axons are short enough to justify setting γ*a* = ∞, which has been termed the local interaction approximation (Robinson et al., 2004). We therefore take only γ*e*, the damping rate of cortical pyramidal neurons, to be finite. This turns all wave equations except the cortical one into delayed one-to-one mappings. The model was implemented on a 5 × 5 grid of nodes with coupling to nearest-neighbor nodes via this damped-wave equation.

# **2.2. SPIKING NETWORK MODEL**

The spiking network was based on several previous models developed by our group (Lytton and Stewart, 2005; Lytton et al., 2008b; Neymotin et al., 2011b; Kerr et al., 2012; Song et al., 2013). It consisted of 4950 event-driven integrate-and-fire neurons. These were divided into three types (excitatory pyramidal cells *E*, fast-spiking inhibitory interneurons *I*, and low-threshold spiking interneurons *IL*), which were in turn distributed across the six layers of the cortex, plus two thalamic cell populations (excitatory thalamocortical relay *TCR* and inhibitory thalamic reticular *TRN*), for 15 distinct neuronal populations in total. The numbers and locations of each neuronal population are illustrated in **Figure 2**, and were as follows: E2 (i.e., excitatory pyramidal neurons of layer 2/3), 1500; I2, 250; IL2, 150; E4, 300; I4, 200; IL4, 150; E5R, 650; E5B, 150; I5, 250; IL5, 150; E6, 600; I6, 250; IL6, 150; TCR, 100; and TRN, 100. The pyramidal neurons

in layer 5 are divided into two populations, *R* (regular firing) and *B* (bursting), since these have different cellular properties and connectivity patterns.

Connectivity (shown in **Figure 3B**) and the relative numbers of neurons per layer were based on published models (Traub et al., 2005; Neymotin et al., 2011a,b) and anatomical studies (Thomson et al., 2002; Binzegger et al., 2004; Song et al., 2005; Lefort et al., 2009; Adesnik and Scanziani, 2010). Connectivity was strongest between populations within a given layer, as seen from the four clusters visible along the diagonal of **Figure 3B**. Overall, excitatory neurons had more projections than inhibitory ones, but inhibitory projections were typically stronger. This balanced excitation and inhibition such that the overall gain of the system (the number of additional output spikes for every additional input spike) was close to unity. Such balance is necessary for avoiding the stable but undesirable states of seizure

**FIGURE 3 | Connectivity of the models.** Color shows normalized effective connectivity (probability × weight) from rows to columns, with red denoting excitation and blue denoting inhibition. **(A)** Connections in the field model (CE, cortical excitatory; CI, cortical inhibitory; TCR, thalamocortical relay; TRN, thalamic reticular nucleus; SD1, striatal D1; SD2, striatal D2; GPi, internal globus pallidus; GPe, external globus pallidus; STN, subthalamic nucleus). **(B)** Connections in the network model. Approximate diagonal symmetry shows that most connections are reciprocal; relatively strong connections along the diagonal indicate high intralaminar connectivity.

(pathologically high firing) and quiescence (pathologically low firing).

Individual neurons were modeled as event-driven, rule-based units. Since computing resources are finite, a tradeoff must be made between the complexity of neurons vs. the complexity of the network. The neuron model used was complex enough to replicate key features found in real neurons, including adaptation, bursting, depolarization blockade, and voltage-sensitive NMDA conductance (Lytton and Stewart, 2005, 2006; Lytton and Omurtag, 2007; Lytton et al., 2008a,b; Neymotin et al., 2011b), yet was simple enough to connect into large (103 −106 neuron) networks.

Each neuron had a membrane voltage state variable (*Vm*) with a baseline value determined by a resting membrane potential parameter (*V*RMP, set at −65 mV for pyramidal neurons and low-threshold-spiking interneurons, and at −63 mV for fastspiking interneurons). This membrane voltage was updated by one of three events: synaptic input, threshold spike generation, and refractory period. These events are described briefly below; further detail can be found in the papers and code cited above.

# *2.2.1. Synaptic input*

The response of the membrane voltage to synaptic input was modeled as an instantaneous rise and exponential decay: *Vn*(*t*) = *Vn*(*t*0) + *ws*(1 − *Vn*(*t*0)/*Ei*)*e* − *t* − *t*0 τ*i* , where *Vn* is the membrane voltage of neuron *n*; *t*0 is the synaptic event time (i.e., *t* − *t*0 is the time since the event); *ws* is the weight of synaptic connection *s*; *Ei* is the reversal potential of ion channel *i*, relative to resting membrane potential (where *i* = AMPA, NMDA, or GABAA; and *E*AMPA = 65 mV, *E*NMDA = 90 mV, and *E*GABAA = −15 mV); and τ*i* is the receptor time constant for ion channel *i* (where τAMPA = 20 ms; τNMDA = 30 ms; and τGABAA = 10 or 20 ms for somatic and dendritic GABAA, respectively).

## *2.2.2. Action potentials*

A neuron fires an action potential at time *t* if *Vn*(*t*) > *Tn*(*t*) and *Vn*(*t*) < *Bn*, where *Vn*, *Tn*, and *Bn* are the membrane voltage, threshold voltage (−40 mV for pyramidal neurons and fast-spiking interneurons, −47 mV for low-threshold-spiking interneurons), and blockade voltage (−10 mV for interneurons and −25 mV for pyramidal neurons), respectively, for neuron *n*. Action potentials arrive at target neurons at time *t*2 = *t*1 + *l*(*n*1, *n*2)/*v* + τ*s*, where *t*1 is the time the first neuron fired, τ*s* is the delay due to synaptic conduction effects, *l*(*n*1, *n*2) is the axon length between neurons *n*1 and *n*2, and *v* is the axonal conduction velocity (1 m·s −1, which is smaller than in the field model, since long-range fibers tend to be more heavily myelinated).

# *2.2.3. Refractory period*

After firing, a neuron cannot fire during the absolute refractory period, τ*A* (10 ms for interneurons and 50 ms for pyramidal neurons). Firing is reduced during the relative refractory period by two effects: first, an increase in threshold potential, *Tn*(*t*) = 1 + *Re*− *t* − *t*0 τ*R Tn*(*t*0), where *R* is the fractional increase in threshold voltage due to the relative refractory period (0.25 for interneurons and 0.75 for pyramidal neurons) and τ*R* is its time constant (1.5 ms for interneurons and 8 ms for pyramidal neurons); and second, by hyperpolarization, *Vn*(*t*) = *Vn*(*t*0) − *He*− *t* − *t*0 τ*H* , where *H* is the amount of hyperpolarization (0.5 mV for interneurons and 1 mV for pyramidal neurons) and τ*H* is its time constant (50 ms for interneurons and 400 ms for pyramidal neurons).

Local field potentials (LFPs) were computed for each cortical layer as the average membrane voltage across all neurons in that layer; after baseline removal and normalization, this approach is roughly equivalent to summing over all synaptic currents (Mazzoni et al., 2010). While this approach does not take into consideration synaptic and dendritic geometry, this is not possible in the event-driven point-neuron model used here.

Simulations were run in NEURON 7.3 (Hines and Carnevale, 2001; Carnevale and Hines, 2006) on a Linux workstation with an Intel Xeon 2.7 GHz CPU; each 20 s simulation took approximately 10 min to run on a single core. To avoid edge effects, the first and last 2 s of simulated data were discarded. All analyses were performed on the remaining 16 s of simulated data. Since the model is at steady-state and does not incorporate plasticity effects, longer runs produced similar results (data not shown). Model parameters were tuned manually (within physiological limits) to match experimentally observed firing rates, dynamics, and information-theoretic properties, as described in Song et al. (2013).

# **2.3. INPUT DRIVE**

The composite model consisted of the spiking network model being driven by ("embedded in") the activity of the field model. Since the field model represents a brain region much larger than the network model, the field causally influences the network, but not vice versa. The key methodological novelty of this work is that the spiking network model is thus embedded in an environment with physiologically realistic dynamics (as provided by the field model), rather than the white noise environment such models are typically embedded in.

To obtain realistic firing rates in the network model, the input spiking rate each neuron receives must be bounded. Hence, the firing rate from each neuronal population in the field model was normalized so that the minimum and maximum input spiking rates were 225 and 1125 s−1 for excitatory neurons and 30% lower for inhibitory neurons. The input drive was obtained by treating each of these normalized instantaneous firing rates as the rate of an ensemble of Poisson processes for generating spikes. These spikes were then used to drive each population of spiking neurons, using the same connections as used in the field model itself (e.g., excitatory cortical neurons in the network model received input from the excitatory cortical field, the inhibitory cortical field, and the thalamic field); relative connection weights were also set to match those of the field model. Thus, each neuron belonging to a given population in the network model receives the same average rate of input from the field model, but from a separate Poisson process, thereby avoiding artificial correlations in input spike times between neurons.

Four different inputs were explored in this work. First, all neurons in the network were driven by spikes drawn from a spectrally white distribution ("WN", the white noise model). This represents the control condition, and is identical to the approach used in previous work with the network model (Neymotin et al., 2011b). Second, neurons were driven by the thalamocortical version of the field model ("TC", the thalamocortical model); i.e., connection strengths to and from the basal ganglia neuronal populations were set to zero. Third, neurons were driven by the full basal ganglia/thalamocortical model described above ("BG", the healthy basal ganglia model). Finally, neurons were driven by the full basal ganglia/thalamocortical model, using parameter values shown by Van Albada and Robinson (2009) to best match the electrophysiological changes associated with the degeneration of dopaminergic projections to the striatum ("PD", the Parkinson's disease model).

# **2.4. ANALYSIS**

# *2.4.1. Fano factor*

The Fano factor, a common measure of spiking variability (Churchland et al., 2010), is the ratio of the variance to the mean of the spike rate:

$$F = \frac{\sigma^2}{\mu},\tag{5}$$

where σ is the standard deviation and μ is the mean of the time series of binned spiking activity across all neurons. To explore spiking variability on a range of different time scales, the time bin size was varied from 1 ms (resulting in 16,000 bins, with an average of roughly 10 spikes per bin) to 8 s (resulting in 2 bins, with roughly 80,000 spikes per bin).

# *2.4.2. Population burst probability*

A population burst (Benayoun et al., 2010) was defined as ≥2 neurons firing within a given 10 ms time bin. The probability of a burst of size *N* was defined as the number of time bins with *N* cells firing divided by the total number of time bins. The relative burst probability was calculated by dividing the observed number of bursts of each size by the number of bursts of that size expected from uncorrelated activity, which in turn was determined via the observed firing rate (averaged over the entire simulation) and the binomial probability distribution.

# *2.4.3. Spectral granger causality*

Information flow was quantified in terms of spectral Granger causality, also called the directed transfer function (Kaminski et al., 2001). Although many alternative tools for inferring causality exist, such as directed transfer entropy (Lizier et al., 2011), no others allow the spectral properties of the signals to be analyzed in detail.

As in standard Granger causality analysis, spectral Granger causality of α(*f*) → β(*f*) is non-zero if prior knowledge of variable α at frequency *f* reduces error in the prediction of β at frequency *f* . The directionality of the causation arises from the fact that Granger causality quantifies how much the history of time series α can be used to predict the future of time series β: if α has a strong causal influence on β, then the prediction error will be reduced.

Spectral Granger causality is calculated by Fourier transforming the multivariate autoregressive model used in standard Granger causality. Hence, the spectral Granger causality from time series α(*t*) to time series β(*t*) is defined as (Cui et al., 2008)

$$G\_{\alpha \to \beta}(f) = -\log\left(1 - \frac{\left(N\_{\alpha a} - \frac{N\_{\beta a}^2}{N\_{\alpha a}}\right) |H\_{\beta, a}(f)|^2}{S\_{\beta, \beta}(f)}\right),\tag{6}$$

where *N* is the noise covariance, *H*(*f*) is the transfer function, and *S*(*f*) the spectral matrix, as derived from the bivariate autoregressive model of α(*t*) and β(*t*). This analysis was performed in Matlab 2012a using code based on the BSMART toolbox, available via http://www.brain-smart.org.

# **3. RESULTS**

The neural field model results were similar to those reported previously (Van Albada and Robinson, 2009; Van Albada et al., 2009), and are briefly presented here for completeness. We then present the overall dynamics of the spiking network model (Kerr et al., 2012), comparing its dynamics for each of the four drives (white noise, the thalamocortical model, the healthy basal ganglia model, and the PD model). Finally, we focus more closely on the alterations that occur in the PD-driven model and their implications. We have split the results into these sections in order to better accomplish our dual goals of (1) presenting the new composite model, and (2) applying this model to help understand the pathophysiology of PD.

# **3.1. FIELD MODEL DYNAMICS**

Firing rates in each neuronal population were similar to those reported previously (Van Albada and Robinson, 2009; Van Albada et al., 2009). Because the drive from the field model to the network model was normalized to a range that provided realistic firing rates in the latter, tonic firing rates had negligible effect on the simulations.

Changes in coherence are a commonly reported finding in PD. In the PD model, coherence between the GPe and the GPi was lost, and high frequency power (>10 Hz) in the GPi was reduced (**Figure 4**). In the healthy state, activity in the GPe and GPi is strongly correlated (*r*2 = 0.9). Following dopamine loss, this correlation is substantially reduced (*r*2 = 0.3). This is because the GPe and GPi are both mostly influenced by the striatum in the healthy state, whereas the GPi is strongly driven by the STN in the parkinsonian state, resulting in strong coherence between the GPi and STN in the PD model. Increased STN-GPi coherence at frequencies up to about 35 Hz has indeed been found in PD off levodopa compared to the on-levodopa condition (Brown et al., 2001). Since the GPi is the only nucleus of the basal ganglia that projects to the thalamus or cortex (**Figure 1**), all changes observed in the network model in the healthy versus PD cases are due to the altered dynamics of the GPi.

To characterize the overall dynamics of the different field models, we looked at their power spectra. In the absence of the basal ganglia, cortical excitatory neurons had a strong alpha peak (10 Hz), and a weaker harmonic in the beta range (20 Hz), as shown in **Figure 4B**. Cortical inhibitory neurons were driven strongly by thalamocortical cells, evident both in the phase locking between the two populations (**Figure 4A**), and in the similarity of their power spectra below 70 Hz (**Figure 4B**). The addition of the basal ganglia (**Figure 4B**, middle panel) reduced the strength of the alpha peak in cortical excitatory neurons and reduced the slope of the power law spectral fall-off at high frequencies; in cortical excitatory neurons, this slope changed from *P*(*f*) ∝ *f* −5.3 to *P*(*f*) ∝ *f* −4.3. Reduced dopamine corresponding to PD reduced the power of higher frequencies (>10 Hz) relative to lower frequencies (<10 Hz) in the cortical, thalamic, and GPi spectra. For example, the GPi showed a 2% decrease in power at 10 Hz compared to a 76% decrease at 20 Hz. In contrast, reduced dopamine increased power in the STN at frequencies >10 Hz (e.g., 2.2 times larger at 20 Hz), a result also reported experimentally (Brown et al., 2001; Cassidy et al., 2002; Priori et al., 2004).

## **3.2. NETWORK MODEL DYNAMICS**

The field drive into the network model strongly modulated its spiking activity (**Figure 5**). Firing rates varied from near zero during the troughs of input activity to >10 Hz during the peaks (**Figure 5A**). The temporal structure of the spiking activity depended strongly on the type of input drive used (**Figure 5B**). As a control, white noise produced no consistent temporal structure. The TC-driven model input produced some structure, with a characteristic time scale below 500 ms. The BG-driven model added some features on longer time scales (of order 1 s) to the activity produced by the TC-only field model. Variability in firing rate, as measured by the Fano factor, was lowest in the WN-driven model (**Figure 5D**)—as would be expected since the white noise had the lowest variability of the four inputs. On time scales <1 s, the PD-driven model had the greatest variability, while the BG-driven model had the greatest variability on scales >1 s.

The power spectra of the network model, shown in **Figure 6A**, were broadly similar to those of the input drives, but with several interesting differences. The basic filter properties of the network model are apparent from the shape of spectrum of the WN-driven model; to a first approximation, the network acts like a low-pass filter, with *P*(*f*) ∝ *f* −4.0 for *f* > 20 Hz. However, actual afferent activity in the brain is already low-pass filtered due to dendritic properties, so a more realistic input (the thalamocortical drive) results in even greater low-pass filtering. For example, the WNdriven model predicts 5.4 times more power at 10 Hz than the TCdriven model. Both BG- and PD-driven models differed markedly from the TC-driven model in the 20–30 Hz band, where many basal ganglia nuclei have their peak power. Interestingly, this peak was much sharper in the network model than in the input drive, demonstrating a resonance effect (compare **Figure 4A** with **Figure 6B**).

To quantify synchrony in the model on a population level, we used population burst size (Benayoun et al., 2010). All of the fielddriven models showed substantially higher population bursting than the WN-driven model (**Figure 6B**). This is because the field drive applies a global modulatory signal to the network, which organizes the firing of its neurons into up and down states (as evident from the bands of spikes in **Figure 5A**); in contrast, the WN-driven model has a constant, intermediate level of activation.

## **3.3. DYNAMICAL CHANGES IN THE PARKINSON'S DISEASE MODEL**

The Parkinson's disease model (PD-driven model) showed a number of changes that suggest possible mechanisms underlying the clinical dysfunctions of the disease. Compared to the healthy control (BG-driven model), the PD-driven model showed a shift in the LFP spectrum toward lower frequencies, with higher delta power and a lower beta peak frequency (**Figure 6**), consistent with clinical findings (Stoffers et al., 2007). These changes were also readily apparent looking at the LFP time series, which showed a flattening of activity between the slow, high-amplitude features (**Figure 5B**). Soikkeli et al. (1991) noted such slowing in 10 out of 18 non-demented PD patients, as well as in all 18 demented PD patients studied [see Figure 1 in Soikkeli et al. (1991)].

The PD-driven model showed an 18 ± 2% decrease in firing rates compared to the healthy model (**Figure 5C**), consistent with changes in fMRI indicators of activity (Monchi et al., 2007). The PD-driven model also showed greater firing variability than the healthy model on most time scales. For example, with a bin size of 1 ms, the Fano factor was 41% higher in the PD-driven model (**Figure 5D**). However, it showed less variability on very long time scales: with a bin size of 8 s, the Fano factor was 2.4 times higher in the BG-driven model. The increased variability in the PD-driven model on all but the longest timescales is consistent with the enhanced oscillations and synchrony associated with PD (Goldberg et al., 2002). Note that maximal dynamical richness does not necessarily correspond to maximal variability in firing rates: for example, tonic firing will have low dynamical richness

and low variability on all time scales, while strong, seizure-like oscillations will also have low dynamical richness, despite very

activity. **(C)** Mean firing rates by cell type (averaged over both cortical and

high variability (at least on the time scale of the oscillation). The concentration of activity in large population bursts was a prominent feature of the PD-driven model. For example, bursts consisting of 40 neurons were 60% more common in the PDdriven model than in the healthy model, while 70-neuron bursts were three orders of magnitude more common. (Population bursts smaller than 30 neurons were more common in the healthy model, a result of its higher firing rate.) Although it is tempting to consider these large population events in the context of parkinsonian tremor, we did not note a clear periodicity in their occurrence.

A crucial question in PD is the mechanism by which information flow is disrupted from higher cortical areas (e.g., those involved in motor planning) to primary areas (e.g., those involved in motor execution). Although information flow between cortical layers is bidirectional, a dominant direction of information flow is suggested by both anatomical and functional studies (Bollimunta et al., 2008). This dominant information pathway is believed to

on all scales.

stream from thalamic inputs to layer 4 (or upper layer 5 in agranular motor cortices), up to layer 2/3 for processing, and thence to layer 5, which in turn produces outputs to multiple sites including the thalamus, basal ganglia, and brainstem. We hypothesized that damage to this dominant pathway would represent a pathology with major functional consequences. We therefore used Granger causality to quantify information flow between the cortical layers that comprise this pathway.

Overall, interlaminar spectral Granger causality was highest in the BG-driven model, and lowest in the WN-driven model (**Figure 7**). Most notably, the BG-driven model showed a prominent peak in causality in the high-beta/low-gamma band (20–35 Hz). This peak was almost entirely absent in the PD-driven model; for example, peak causality from layer 4 to layer 5 in this frequency range was only half that of the BG-driven model (0.23 and 0.45 for PD- and BG-driven models, respectively), even though these models had similar spectral power (**Figure 6A**). As shown in **Figure 7**, similar results were seen in other layer pairs (e.g., 4→2/3, 2/3→5, and 6→2/3).

# **4. DISCUSSION**

We have explored the effects of driving a spiking network model with several different types of input, including those corresponding to the healthy brain and to PD. Many of the differences between the healthy and PD models accord with prior experimental findings. For example, we found a modest but consistent reduction in firing rates of cortical neurons in PD. Although there are no direct studies of cortical firing rates during PD in humans, several indirect measures from functional imaging suggest such a decrease (Jenkins et al., 2004; Monchi et al., 2004, 2007). We also found a shift toward lower LFP frequencies, a finding consistent with PD electroencephalography (Soikkeli et al., 1991; Bosboom et al., 2006; Stoffers et al., 2007). We found increased synchrony between neurons in our PD model, as measured by population burst size and probability; increased synchrony among basal ganglia neurons is a commonly reported finding in PD (Raz et al., 1996), and increased synchrony among cortical neurons has also been reported (Goldberg et al., 2002).

Our major finding was the loss of Granger causality between cortical layers in the high-beta/low-gamma band. The Granger causality for the PD-driven model was more similar to the TCand WN-driven models than to the BG-driven model, suggesting that the dynamical properties of the basal ganglia that facilitate cortical information flow are almost entirely lost in PD. The frequency range of this disrupted information flow is thought to be crucial for encoding motor commands, especially limb movements (Van Der Werf et al., 2008; Muthukumaraswamy, 2010). Gamma has also been implicated in many cognitive processes (Fries et al., 2007), including the perceptual binding underlying sensorimotor coordination (Lee et al., 2003) and consciousness (Llinas et al., 1998). Hence, our observation of disrupted causality might also partially account for some of the cognitive symptoms of PD, including bradyphrenia and planning deficits (Morris et al., 1988; Chaudhuri and Schapira, 2009).

The fact that Granger causality was disrupted in the PDdriven model (**Figure 7**) while the power spectrum was nearly unchanged in the same frequency band (**Figure 6A**) shows that the changed input drive has reorganized the dynamics of the network in complex ways. Since the GPi does not project directly to the cortex, these changes are entirely mediated by the thalamus; indeed, thalamic lesions alone are sufficient for producing parkinsonian symptoms in rats (Oehrn et al., 2007). Since the thalamus projects differentially to the different layers of the cortex, a major change in thalamic input is sufficient explanation for why the causality would shift so dramatically. Specifically, the thalamus normally projects strongly to layer 4; the peak in causality at 20–35 Hz is consistent with thalamic modulation by the GPi. In PD, inhibition to the thalamus is increased, which results in weaker drive to the cortex and thus a loss of information flow.

layer 4 to 2/3 in the delta (<5 Hz) and high-beta/low-gamma (20–35 Hz) bands; causality in the latter band is almost entirely lost in Parkinson's disease. **(B)** The causality from layer 2/3 to layer 5 is slightly reduced in this band in Parkinson's disease. **(C)** These two effects combine to

layer 6 to layer 2/3. In each case, the high-beta/low-gamma band Granger causality is significantly higher in the BG-driven model than in any of the other models.

Our findings suggest that therapeutic interventions, such as deep brain stimulation (Deuschl et al., 2006), may be more effective if they restore both the dynamics and the tonic level of activity of the GPi, rather than just the latter.

Several of our findings are qualitatively consistent with experimental results pointing to a loss of complexity in EEG time series from patients with a variety of cognitive disorders, including PD (Stam et al., 1994, 1995; Vaillancourt and Newell, 2002). For example, in the healthy model, the slope of the Fano factor increases roughly linearly on time scales from 1 ms to 10 s, indicating dynamical structure across a wide range of time scales (**Figure 5D**). This result can be seen qualitatively in the LFP time series of the healthy model, which appeared to show meaningful structure over a broader range of time scales than any of the other models (**Figure 5B**). We speculate that these properties may reflect the number of possible states that the network can assume, which may in turn be related to the number of different motor programs that can be implemented by the network. This principle is closely related to the concept of φ, defined as "the repertoire of causal states available to a system as a whole" (Balduzzi and Tononi, 2008). While φ cannot be easily computed for moderately large networks such as ours, we expect that it will be manifested in terms of the network's ability to perform real motor tasks—a topic we will explore in future work. Specifically, we predict that the BG-driven model will perform better on simulated reaching tasks than the WN-, TC-, or PD-driven models.

Beta-band activity (15–30 Hz) was predominantly generated by the thalamic and inhibitory cortical neuronal populations in our model (**Figure 4B**, top panel), in agreement with previous experimental and modeling studies (Brown and Williams, 2005; Hahn and McIntyre, 2010). Most empirical studies of beta activity in PD have focused on the basal ganglia nuclei, with increased power in the STN being a commonly reported finding (Brown and Williams, 2005; Kühn et al., 2005; Weinberger et al., 2006). In our model, we found that beta power in the STN was indeed enhanced in PD (**Figure 4B**, middle and bottom panels), which may reflect an idling or antikinetic state (Brown and Williams, 2005; Engel and Fries, 2010).

## **4.1. LIMITATIONS**

Several experimentally observed features of PD, such as increased coherence among neurons in the STN, can only be explicitly represented using a neuronal network model of the basal ganglia (Terman et al., 2002)—a major benefit of that modeling approach. However, it is not known whether these phenomena are causally linked to parkinsonian symptoms. Hence, in the present context, the benefits of using a neural field model for the basal ganglia outweigh the drawbacks of this approach. In future, a spiking network model of the basal ganglia would be desirable in order to account for these and other phenomena, such as reinforcement learning. An explicit representation of dopamine in such a mode—rather than the implicit representation used here would also allow the effects of pharmacological interventions to be modeled directly.

Due to the eloquence of the motor system, movement disorders are the most obvious symptoms of PD. Yet the pathophysiology of the disease is widespread; even the retina is affected (Hajee et al., 2009). We stress that the spiking network model used here was designed as a model of association cortex, not primary motor cortex; for example, our model includes layer 4 cells, which are absent from the latter. However, since the thalamus and striatum have broad projections to the cortex, we expect the dynamical and information-theoretic changes in PD (such as increased synchrony and reduced complexity) to extend to motor areas as well. In the future, we will explore the effects of PD in a model of primary motor cortex controlling a virtual arm (Chadderdon et al., 2012), with the aim of directly demonstrating classical parkinsonian motor symptoms. By incorporating sensory feedback into this model, the white noise that was used to drive the neural field component can be replaced with more realistic input, thereby addressing another obvious limitation of the method used here.

# **4.2. MULTISCALE DYNAMICS IN A COMPOSITE MODEL**

To our knowledge, this work represents the first composite spiking network/neural field model of the brain. This is a multiscale model that spans spatial scales from 10μm to 30 cm and temporal scales from 1 ms to tens of seconds. The composite method provides a way of linking two types of models that provide access to different spatial scales—a network model than spans scales from individual neurons (10μm) to a cortical column (600μm), and a field model encompassing the whole diencephalon (30 cm).

# **REFERENCES**

Temporally, both network and field models are valid over many orders of magnitude (approximately 10−3 −104 s).

The mechanism used here to couple the field and network models is just one of several alternatives (Wilson et al., 2012). In the present case, the coupling was unidirectional; the network model did not affect the dynamics of the field model. While this can be easily justified in terms of the effective size of each model, an alternative approach generates the neural field based on the dynamics of the network model, using the new neuronin-cell approach of Robinson and Kim (2012). Because spiking network models are still limited in their capacity to generate accurate dynamics on a large scale, this approach cannot yet be used in place of neural field models. However, this may change if scientific advances and improved computing facilities enable the development of larger and more realistic spiking network models.

Many spiking network models that are too small to show selfsustaining activity are driven by white noise (Hill and Tononi, 2005; Vogels and Abbott, 2005; Oswald et al., 2009; McDonnell et al., 2011; Volman et al., 2011; Kerr et al., 2012; Muller and Destexhe, 2012; Vijayan and Kopell, 2012). Here we demonstrated that using physiologically realistic input instead of white noise has a major impact on multiple measures of network activity, including power spectra, spiking variability, burst probability, and Granger causality. Thus, white-noise-driven spiking network models are an abstraction away from the physiological environment, and should perhaps be considered as being analogous to artificially driven slice preparations rather than *in vivo* activity.

## **ACKNOWLEDGMENTS**

Cliff C. Kerr, Samuel A. Neymotin, George L. Chadderdon, and William W. Lytton were supported by DARPA grant N66001-10- C-2008; P. A. Robinson was supported by the Australian Research Council and the Westmead Millennium Institute; and Sacha J. van Albada was supported by EU FP7 Grant 269921 (BrainScaleS) and the Helmholtz Alliance on Systems Biology. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

human. *Clin. Neurophysiol.* 116, 2510–2519.

cortex. *PLoS ONE* 7:e47251. doi: 10.1371/journal.pone.0047251

nucleus diminishes thalamocortical relay errors in a biophysical network model. *Neural Netw.* 24, 602–616.

D. (2006). Competition between feedback loops underlies normal and pathological dynamics in the basal ganglia. *J. Neurosci.* 26, 3567–3583.

topology, delays and short-term plasticity. *Brain Res.* 1434, 162–177.

dynamics in parkinsonian brain: the boundary between synchronized and nonsynchronized dynamics. *Phys. Rev. E* 83:042901. doi: 10.1103/PhysRevE.83.042901

characteristic of Parkinson's disease without dementia. *Brain* 130, 1847–1860.

trauma determines the threshold for epileptic activity in a model of cortical deafferentation. *Proc. Natl. Acad. Sci. U.S.A.* 108, 15402–15407.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 22 January 2013; accepted: 02 April 2013; published online: 25 April 2013.*

*Citation: Kerr CC, Van Albada SJ, Neymotin SA, Chadderdon GL, Robinson PA and Lytton WW (2013) Cortical information flow in Parkinson's disease: a composite network/field model. Front. Comput. Neurosci. 7:39. doi: 10.3389/fncom.2013.00039*

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

# The mesoscopic modeling of burst suppression during anesthesia

# **David T. J. Liley \* and MatthewWalsh**

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

#### **Edited by:**

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

#### **Reviewed by:**

Axel Hutt, INRIA CR Nancy, France Moira Steyn-Ross, University of Waikato, New Zealand

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

David T. J. Liley, Brain and Psychological Sciences Research Centre, Faculty of Life and Social Sciences, Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia. e-mail<EMAIL_ADDRESS>
The burst-suppression pattern is well recognized as a distinct feature of the mammalian electroencephalogram (EEG) waveform. Consisting of alternating periods of high amplitude oscillatory and isoelectric activity, it can be induced in health by deep anesthesia as well as being evoked by a range of pathophysiological processes that include coma and anoxia. While the electroencephalographic phenomenon and clinical implications of burst suppression have been studied extensively, the physiological mechanisms underlying its emergence remain unresolved and obscure. Because electroencephalographic bursting phenomenologically resembles the bursting observed in single neurons, it would be reasonable to assume that the theoretical insights developed to understand bursting at the cellular ("microscopic") level would enable insights into the dynamical genesis of bursting at the level of the whole brain ("macroscopic"). In general action potential bursting is the result of the interplay of two time scales: a fast time scale responsible for spiking, and a slow time scale that modulates such activity. We therefore hypothesize that such fastslow systems dynamically underpin electroencephalographic bursting. Here we show that a well-known mean field dynamical model of the electroencephalogram, the Liley model, while unable to produce burst suppression unmodified, is able to give rise to a wide variety of burst-like activity by the addition of one or more slow systems modulating model parameters speculated to be major "targets" for anesthetic action. The development of a physiologically plausible theoretical framework to account for burst suppression will lead to a more complete physiological understanding of the EEG and the mechanisms that serve to modify ongoing brain activity necessary for purposeful behavior and consciousness.

**Keywords: burst suppression, anesthesia, electroencephalogram, mean field model, neuronal hyperexcitability**

# **1. INTRODUCTION**

Prior to the development of the modern intensive care unit in the early 1960s, that featured intubation, artificial respiration, and comprehensive physiological monitoring, reports of the electroencephalographic pattern of burst suppression (BS) were confined to animal studies involving deep anesthesia and the occasional case of psychosurgery (Niedermeyer, 2009). Since then the burstsuppression pattern has become well recognized as a major diagnostic feature of the EEG waveform that is encountered in a range of encephalopathic conditions, in addition to its appearance in health during deep anesthesia. Typically the BS pattern consists of bursts of high amplitude slow, sharp, or spiking electroencephalographic activity separated by periods of electroencephalographic suppression (isoelectricity). The oscillatory features of the bursts, together with their duration and the duration of suppressed periods show a high degree of variability (see **Figure 1** for examples) that presumably reflects its myriad of initiating causes. First identified during deep anesthesia with tribromoethanol in cats (Derbyshire et al., 1936), labeled burst-suppression pattern by Swank and Watson (1949) during barbiturate and ether anesthesia in dogs, it is now associated with cortical deafferentation (Henry and Scoville, 1952), cerebral anoxia and hypoxia, various types of intracortical lesions (Fischer-Williams and Cooper, 1963), deep

coma, various infantile encephalopathies, the final stages of deteriorated status epilepticus (Treiman et al., 1990), hypothermia, and high levels of many sedative and anesthetic agents (Schwartz et al., 1989; Akrawi et al., 1996).

Burst suppression in the absence of anesthesia is in general associated with a very poor prognosis. For example in neonates (Grigg-Damberger et al., 1989) the appearance of BS, even if transient,is a portent of death or severe neurodevelopmental disability. In contrast, in adult populations while an anoxic/hypoxic BS pattern signals a serious pathophysiological event the outcome is not necessarily fatal and recovery with or without severe neurological damage is possible (Niedermeyer, 2009). Consistent with this are results of experimental work with EEG monitoring in rats revealing that animals with greater rates of high amplitude bursts have a better survival and neurological outcome compared to those with lower rates of low amplitude bursts (Geocadin et al., 2002).

While the electroencephalographic phenomenon and clinical implications of BS have been studied extensively (Brenner, 1985; Niedermeyer, 2009) the physiological mechanisms underlying its emergence remain in general unresolved and obscure.

Burst suppression is typically thought to be spatially homogeneous with burst onset and termination reported to occur near simultaneously across the entire scalp (Brenner, 1985; An et al.,

(1996)]. **(B)** Electroencephalogram recorded in acute anoxia showing a clear burst-suppression pattern with grouped spikes [figure reproduced with permission from Hockaday et al. (1965)]. **(C)** BS pattern during closed loop target controlled propofol infusion at a target level of

infusion [figure reproduced with permission from Lukatch and MacIver

amplitude spikes interspersed with isoelectric periods of variable length [figure reproduced with permission from Voss et al. (2006)]. **(E)** Electroencephalogram recorded from a 3-month-old infant suffering from infantile myoclonic encephalopathy [reproduced with permission from Niedermeyer (2005)].

1996; Ching et al., 2012), indicating that low level subcortical mechanisms may be playing a decisive role. However arguing against this is the fact that this pattern persists subsequent to cortical deafferentation (Lukatch and MacIver, 1996), indicating that it probably represents an intrinsic, though physiologically abnormal, dynamical mode of cortex. Indeed the phenomenal resemblance of the patterns of BS to disorders of neuronal hyperexcitability suggests the involvement of similar physiological mechanisms. For example the bursting during burst suppression is often associated with myoclonic jerks resembling those seen during epileptic fits. Like generalized epileptiform activity, bursts are recorded simultaneously at multiple electrode derivations, implying the wide synchronization of neuronal activity.

At the cellular level a commonly reported finding is that hyperpolarization of the membrane potential of cortical neurons reliably precedes any overt electroencephalographic activity of BS (Steriade et al., 1994). Such hyperpolarization, which has been attributed to an increase in neuronal membrane potassium conductance (Steriade et al., 1994), has been hypothesized to play a major role in the induction of BS. This implied importance of inhibition in the genesis of BS is further supported by results involving rat neocortical brain slice micro-EEG preparations in which the application of a direct acting GABA*A* agonist, muscimol, readily induces BS. However contradicting this result is work reporting that inhibition is diminished during isoflurane-induced BS, in an *in vivo* feline preparation, as evidenced by increases in cortical neuronal input resistance and extracellular chloride concentration (Ferron et al., 2009). Of course it may be that slow periodic modulations in inhibition, rather than singular increases or decreases in inhibition, underpin BS. In support of this view is the recent model of Ching et al. (2012), in which alterations in brain metabolism, due to the effects of hypoxia or anesthesia, parametrically regulate an activity dependent slow modulation of an adenosine triphosphate-gated potassium channel conductance to give rise to BS. However modulations in inhibitory activity alone may not be sufficient to account for BS and more consideration might need to be given to other mechanisms. For example Kroeger and Amzica (2007) present empirical evidence suggesting that modulations in excitatory synaptic efficiency, due to the progressive depletion of interstitial calcium during the periods of high amplitude electroencephalographic activity and its recovery during isoelectric periods, might account for BS. Consistent with this are reports involving laboratory slice preparations in which burst suppression induced by thiopental, propofol, and isoflurane is abolished by the application of glutamate receptor antagonists (Lukatch and MacIver, 1996).Whatever the pathophysiology of BS is it is reasonably clear that it is unlikely to be accounted for by a unitary physiological perturbation. That the physiological factors identified to date in BS all lead to a single well defined state suggests the possibility of an unifying dynamical mechanism. Thus the best hope for progress in understanding the phenomena of BS may be theoretical.

How might we theoretically approach BS? The well studied dynamical mechanisms of action potential bursting (Izhikevich, 2007) may be able to provide vital insights into the mechanisms of bursting in the EEG. In general the dynamical mechanisms underlying bursting can be divided into two broad classes (i) fast-slow

bursters in which there is a clear separation of the underlying time scales, with a fast system responsible for the fast spiking, and a slow system its slow modulation, and (ii) "hedgehog" bursters (Izhikevich, 2000) in which there is no clear separation of time scales. In terms of developing a theory of BS the former might represent the preferable starting point as the little empirical evidence that is available (Ching et al., 2012), at least in humans, suggests that alphoid activity, indicative of normal resting EEG, is preserved during the bursts of BS. Thus a theoretical starting point to understanding BS might be to consider the slow modulation of a dynamical system developed to describe the resting EEG.

One such dynamical system is the mesoscopic electrocortical model of Liley et al. (Liley et al., 1999, 2002, 2011; Bojak and Liley, 2005; Frascoli et al., 2011). This model is capable of accounting for a range of resting electroencephalographic phenomena that includes the alpha rhythm (Liley et al., 2002), the modulation of resting activity by sedative and anesthetic action (Bojak and Liley, 2005) as well as the proconvulsant properties of the latter (Liley and Bojak, 2005), all within a physiologically plausible/admissible parameter space. This model is therefore well suited as a foundationfrom which to explore the physiological and dynamical genesis of BS. However, because in this model rhythmogenesis emerges from a strong coupling between cortical excitatory and inhibitory population activity, in its present form it has a restricted ability to exhibit BS through the parametric separation of time scales, either through the simulated actions of anesthetics or through other parametric routes. Here we show that BS can however emerge in this model by the addition of a slow system driven by one or more of the originally defined mean fields. We speculate that such a slow system represents a mathematical ansatz for the slow neuromodulation of activity by a variety of intracortical, inter-cortical, and subcortical systems that include thalamus and the ascending neurotransmitter modulatory systems.

# **2. MATERIALS AND METHODS**

# **2.1. MESOSCOPIC MEAN FIELD MODELING OF ELECTROCORTICAL ACTIVITY**

The electroencephalogram and electrocorticogram arise out of the cooperative activity of many thousands of neurons. A single electroencephalographic electrode records the synaptically induced currents of well over a 100,000 neurons (Nunez and Srinivasan, 2005) and thus detailing each neurons contribution to this summed activity would appear superfluous. For this reason it is preferable instead to model the activity of populations of neurons. One general way of achieving this, in which known stochastic fluctuations can be included, is to dynamically evolve the probability distributions associated with the states of the neuronal ensemble. While in principle providing a rigorous way forward the formulation of such stochastic equations of motion entails a great deal of physiological uncertainty. For this and other reasons (Deco et al., 2008) a more resolute path is to dynamically evolve some average quantity such as the mean soma membrane potential or the mean firing rate of some suitably defined neuronal ensemble. In this manner a mesoscopic level model can be developed which acts as a bridge between cellular (or microscopic) level activity and whole brain (or macroscopic) level behavior. While the current mathematical approach for formulating the equations of motion for

the activity of neuronal populations or masses, stems principally from the works of Wilson and Cowan (1972, 1973), Nunez (1974), Freeman (1975), and Amari (1975), they are not particularly successful in articulating the genesis of rhythmic activity in the EEG and its modulation by pharmacological agents, due to a range of mathematical simplifications that cannot be justified by an appeal to physiology. For this reason a range of biologically more refined neuronal populations models have been developed (Deco et al., 2008) that form a more appropriate basis from which to model the physiological genesis of the electroencephalogram. The model of Liley et al. (Liley et al., 1999, 2002, 2011; Bojak and Liley, 2005; Frascoli et al., 2011) is one such model and will be the focus of our attempts to account for the phenomenon of BS.

The model of Liley et al. is able to account for the major features of the mammalian electroencephalogram in the context of a parameterization that accords with physiological measurement and experiment. Because it aims to provide a dynamical description of the electroencephalogram the model is spatially coarse-grained over the approximate extent of a cortical macrocolumn, reflecting the general physiological wisdom that such columniation signifies populations of neurons having similar functional behavior and anatomical connectivity. The multiple interactions between individual neurons are replaced by effectivefeed-forward andfeedback interactions between the mean activity (or mean fields) of excitatory and inhibitory populations of neurons over the spatial extent of the column. In this way mammalian neocortex is conceived as consisting of localized populations of excitatory and inhibitory neurons interacting by all possible combinations of feed-forward and feed-back connectivity interacting globally via long-range excitatory connections (**Figure 2**).

Thus the response of the mean soma membrane potential *hk* (*k* = *e, i*) at position *r* on a two-dimensional cortical sheet, in response to induced post synaptic activity *Ilk* (*l* = source, *k* = target population) is given by

$$\pi\_k \frac{\partial h\_k(r, t)}{\partial t} = h\_k^r - h\_k(r, t) + \sum\_{l=e, i} \frac{h\_{lk}^{eq} - h\_k(r, t)}{|h\_{lk}^{eq} - h\_k^r|} I\_{lk}(r, t) \quad (1)$$

The postsynaptic response to a single pre-synaptic action potential (at *t* = 0) is modeled by the well-known synaptic alpha function of cable theory as 0*lk*γ *lkt* exp(1 − γ *lkt*)2(*t*) where 0*lk* is peak amplitude (occurring at *t* = *t* peak = 1/γ *lk*) of the respective excitatory (*l* = *e*) or inhibitory (*l* = *i*) postsynaptic potential (PSP), and 2(*t*) is the Heaviside step function. Thus we assume that the time course of the synaptically induced excitatory and inhibitory currents is described by a critically damped oscillator driven respectively by the mean rate of incoming excitatory and inhibitory axonal pulses:

$$\left(\frac{\partial}{\partial t} + \mathcal{y}\_{\mathbb{lk}}\right)^2 I\_{\mathbb{lk}}(r, t) = \exp(1) \Gamma\_{\mathbb{lk}} \mathcal{y}\_{\mathbb{lk}} A\_{\mathbb{lk}}(r, t), \qquad (2)$$

with

$$A\_{ek}(r,t) = N\_{ek}^{\beta} S\_t[h\_e(r,t)] + \phi\_{ek} + p\_{ek}(r,t), \tag{3}$$

and

$$A\_{ik}(r,t) = N\_{ik}^{\beta} \mathcal{S}\_i[h\_i(r,t)] + p\_{ik}(r,t), \tag{4}$$

where *Alk* comprises the different sources of pre-synaptic spikes: *N* β *lk Si* (input from local cortical neuronal populations), φ*ek* (input from long-range excitatory cortico-cortical fibers), and *plk* (extracortical sources).While the present consensus is that extra-cortical sources (thalamo-cortical afferents) are purely excitatory in nature and thus *pik* = 0, we choose to retain these terms as when time independent they can be utilized to include the effects of tonic inhibition that are known to be induced by anesthetic action. The time courses of the synaptically induced currents, 1/γ *lk* are taken to describe the time course of "fast" excitatory [*l* = *e*: α-amino-3 hydroxl-5-methyl-4-isoxazole-propionate (AMPA) and kainate] and inhibitory [*l* = *i*: γ -amino-butyric-acid type A (GABAA)] neurotransmitter kinetics. Thus each type of PSP (excitatory, inhibitory) is described by two parameters 0*lk*, γ *lk*. However, as we will describe later, a parametrically more flexible description of the PSP is required to meaningfully model the effects of anesthetics in which we can independently vary peak amplitude, rise (*t* peak) and decay times. Mean neuronal population firing rates, *Sl* , are assumed to instantaneous sigmoid functions of the mean soma membrane potential i.e.,

$$S\_l[h\_l(r,t)] = S\_l^{\max} / \{1 + \exp\{\sqrt{2}(h\_l(r,t) - \bar{\mu}\_l/\sigma\_l)\}\} \quad (5)$$

The axonal pulses φ*ek* propagated by the exclusively excitatory long-range cortico-cortical fiber system is in the simplest case described by the following two-dimensional telegraph equation,

$$\begin{split} \left(\frac{\partial}{\partial t} + \nu\_{\rm ek} \Lambda\_{\rm ek}\right)^2 \phi\_{\rm ek}(r, t) - \frac{3}{2} \nu\_{\rm ek}^2 \nabla^2 \phi\_{\rm ek}(r, t) \\ = \nu\_{\rm ek}^2 \Lambda\_{\rm ek}^2 N\_{\rm ek}^\alpha \mathcal{S}\_{\rm e}[h\_{\rm e}(r, t)] \end{split} \tag{6}$$

where *N* α *ek* is the total number of excitatory connections formed by long-range cortico-cortical axons on long on local population *k*, and assumes a single axonal conduction velocity *vek* and an exponential fall off with distance (characteristic scale = 1/3*ek*) of the strength of cortico-cortical connectivity. For simplicity, and given thefact that at least in anesthesia BS appears to have a degree of spatial uniformity, we chose to only study the spatially homogeneous case, i.e., 52 = 0.

Equations (1)–(4) and (6) represent a system of 8 coupled non-linear partial differential equations that typically defines the Liley model of electrocortical rhythmogenesis, which is capable of reproducing the main features of spontaneous human electroencephalogram (alpha resonance,"1/*f*" activity). **Table 1** summarizes all model parameters, their definitions, and approximate ranges.

## **2.2. MODEL PARAMETERIZATION: GENERATION OF NORMATIVE PARAMETER SETS**

Because BS activity (at least that induced by anesthetic and sedative action) is assumed to ultimately arise out of a background of normal electroencephalographic activity it is important to define parametrically normative states in order to study how they may be perturbed during health and disease. We therefore chose to

**Table 1 | List of spatially averaged parameters for different types k, l** ∈ **{e, i} of neuronal target populations in the electrocortical model of Liley et al. (Liley et al., 1999, 2002, 2011; Bojak and Liley, 2005; Frascoli et al., 2011), with typical ranges that are assumed to be**

**physiologically admissible.**

Table adapted from Liley et al. (2011). EPSP, excitatory PSP; IPSP, inhibitory PSP.

utilize previously defined parameter sets (Bojak and Liley, 2005) that have the following properties: (i) are confined to the physiologically admissible parameter space (see **Table 1**), (ii) give rise to electroencephalographically and physiologically plausible activity ("1/*f*" decay at low frequencies plus a relatively sharp peak at alpha frequencies, 8–13 Hz; mean excitatory/inhibitory neuronal firing rates <20 s−1 ) and (iii) that exhibit transient increases in total power and monotonic reductions in mean frequency with respect to modeled anesthetic action (see below). In general such sets can be found by randomly searching the high dimensional physiologically admissible (and plausible) parameter space. For further details see Bojak and Liley (2005).

## **2.3. MODELING ANESTHETIC ACTION**

The range of molecular and cellular targets identified to date as sites of anesthetic action is so varied that a unitary biological mechanism for anesthetic effect seems unlikely. Nevertheless, at least functionally, at the level of cortex anesthetics seem to act principally by enhancing the actions of inhibitory activity (Liley et al., 2011). Indeed from the perspective of the mean field model we have described many of its parameters can be related in a fairly straightforward way to these alterations in inhibitory activity and other identified sites of anesthetic action in cortex (see **Table 2**). However a parametrically more flexible description of the PSP, than is presently incorporated, is required to meaningfully model the effects of anesthetics in which we can independently vary peak amplitude, rise (*t* peak) and decay times. For example isoflurane, a volatile halogenated anesthetic, has been shown to prolong the decay time of the unitary IPSP without altering its time to peak. Fortunately a simple modification of the equation describing the dynamics of the PSP enables independent adjustment of the peak amplitude, rise (*t* peak) and decay times. By defining *Ilk* to satisfy

$$
\left[\frac{\partial}{\partial t} + \mathcal{y}\_{\mathbb{lk}}(\varepsilon\_{lk})\right] \left[\frac{\partial}{\partial t} + \tilde{\mathcal{y}}\_{\mathbb{lk}}(\varepsilon\_{lk})\right] I\_{\mathbb{lk}}(r, t)
$$

$$
= \tilde{\mathcal{y}}\_{\mathbb{lk}}(\varepsilon\_{lk}) \exp[\mathcal{y}\_{\mathbb{lk}}(\varepsilon\_{lk})/\mathcal{y}\_{lk}^{0}] \Gamma\_{\mathbb{lk}} A\_{lk}(r, t),
\tag{7}
$$

$$
\gamma\_{\rm lk}(\varepsilon\_{\rm lk}) = \varepsilon\_{\rm lk} \gamma\_{\rm lk}^0 / (e^{\varepsilon\_{\rm lk}} - 1), \quad \tilde{\gamma}\_{\rm lk} = \gamma\_{\rm lk}(\varepsilon) e^{\varepsilon\_{\rm lk}} \tag{8}
$$

**Table 2 | Relationship between major experimentally identified sites of cortical anesthetic action and parameters of the electrocortical model of Liley et al. (Liley et al., 1999, 2002; Bojak and Liley, 2005).**

\*Parameters will depend on membrane potential in this case.

† Effect demonstrated in periphery, speculative in cortex (Swindale, 2003). where 1/γ 0 *lk* defines the time to peak, we can control the decay of the unitary PSP by altering ε*lk* > 0. Increasing ε*lk* will monotonically increase the decay time of the tail of the unitary PSP (see lower left panel, **Figure 6**). Empirically it is found that increasing the aqueous concentration of a range of GABAergic anesthetic agents leads to a progressive increase in the decay time of the unitary inhibitory PSP (e.g., Banks and Pearce, 1999) and thus ε*ik* will be a monotonic function of anesthetic concentration *c*, i.e., ε*ik*(*c*). Liley et al. (2011), based on a range of empirical evidence, have numerically estimated ε*ik*(*c*) for the volatile anesthetic agent isoflurane. However because ε*ik*(*c*) is not currently known for other GABAergic anesthetic agents we will assume that in general ε*ik* ∝ *c*.

It is worth noting that equation (7) reduces to equation (2) as ε*lk* →0. Further details regarding this formulation can be found in Bojak and Liley (2005).

## **2.4. DEFINING A THEORETICAL BASIS FOR BURSTING**

We call a bursting system *fast-slow* if it can be written in the following form

$$
\dot{\mathbf{x}} = f(\mathbf{x}, \mathbf{y}) \qquad \text{(fast oscillatory system)} \tag{9}
$$

$$
\dot{\chi} = \mu \lg(\mathfrak{x}, \mathfrak{y}) \quad (slow \text{ } modulator \text{ } system) \tag{10}
$$

where *x* ∈ R *m* describes the m-dimensional system responsible for the fast oscillatory (spiking in single neuron models) dynamics and *y* ∈ R *n* the n-dimensional slow system that modulates the fast oscillations (or spiking behavior). The parameter µ represents the ratio of the time scales between the slow and fast system. It is typically assumed that µ 1. Because µ can be made as small as we like equations (9) and (10) represent a singularly perturbed system.

We will assume that equations (1)–(4) and (6), which putatively describe the genesis of the "fast oscillatory" resting EEG, correspond to the m-dimensional fast system. To motivate the slow n-dimensional slow modulatory system we will make a plausible appeal to the biophysics of synaptic resource depletion and recovery during periods of sustained neuronal population activity. During periods of high firing neural activity a variety of factors come into play to diminish synaptic efficiency. The most important of these are receptor desensitization and synaptic vesicle depletion. Tsodyks and Markram (1997) developed a model to account for the biophysics of such activity dependent short term synaptic depression estimating that its onset is rapid, of the order of milliseconds, but that its recovery is quite slow, of the order of 800 ms. Given that such a time scale is approximately at least an order of magnitude greater than the characteristic time scales associated with resting EEG activity, this may represent a candidate slow EEG modulatory system. On this basis we choose to include this activity dependent short term synaptic depression using the following two different formulations, referred respectively to as SS1 and SS2,

$$
\dot{\Gamma}\_{lk} = \mu\_l [\theta\_l - k\_l \mathbb{S}\_l(h\eta)], \quad \text{(SS1)} \tag{11}
$$

$$\dot{\Gamma}\_{lk} = \mu\_l \| \Gamma\_{lk}^0 / (1 + \exp[\kappa\_l(h\_l - \xi\_l)]) - \Gamma\_{lk} \|, \quad (\text{SS2}) \tag{12}$$

where 1/µ*l* is the characteristic time scale of the respective slow modulatory system. Based on Tsodyks and Markram (1997) we will fix 1/µ*l* to 1000 ms. The advantage of the first formulation is that the rates of synaptic recovery (µ*l*θ*l*) and depletion (µ*lkl*) can be separately specified. The advantage of the second formulation is that 0*lk* remains bounded between 0 0 *lk* (the resting value) and zero decreasing monotonically with increasing mean soma membrane potential *hl* , with 0 0 *lk* remaining as a free parameter. For low levels of the respective neuronal activity (*hl*) there is very little decrease in 0*lk* until a threshold ξ *l* is reached, with the parameter κ*l* controlling the sensitivity of the change at this threshold to variations in neuronal activity. This formulation has previously been used by Tabak and Rinzel (2005) in their mean field model for spontaneous electrical bursting activity in embryonic chick spinal cord.

#### **2.5. COMPUTATIONAL METHODS**

All numerical integrations and one-dimensional dynamical continuations were performed using the XPPAUT package (Ermentrout, 2002). A 4-th order Runge-Kutta scheme with a time step of 0.1 ms was used to numerically integrate the differential equations. Because of the multiple time scales our differential system may suffer from stiffness and numerical solutions may not converge. In these cases we have used the recommended "stiff" integrator CVODE (Cohen and Hindmarsh, 1994) as implemented in XPPAUT.

# **3. RESULTS**

All numerical simulations were performed using a single model parameter set having a physiologically plausible white noise fluctuation spectrum (see top left panel **Figure 3**) and a single stable fixed point. This parameter set was found using the methods described in section 2.2. The parameters used, all within the physiologically admissible domain, can be found in the **Table 3**.

**Figure 3** shows the effects of the activity dependent modulation of 0*ii* on simulated mean field EEG activity using SS1 [equation (11)]. Bursts emerge periodically, with intervening near isoelectric intervals, apparently driven by slow variations in 0*ii* (red line, top right panel). A spectrogram of a sufficiently long simulated time series reveals that the frequency of the model EEG activity decreases from low beta (∼= 15 Hz) to high alpha (∼= 12 Hz) over the period of the bursts. Such intra-burst "chirping" is a common feature of many of the topologically identified single neuronal bursters (Izhikevich, 2007). Of interest are the multiple harmonics of this dominant oscillatory activity.

**Figure 4** illustrates how we might dynamically account for the transition to, and cessation from, bursting and follows the now standard method of the dissection of neural bursting pioneered by Rinzel (1985). Here we have set µ*i* = 0 and consider how the "fast" EEG system responds. **Figure 4** shows a one-dimensional bifurcation diagram of this "fast" EEG system with 0*ii* as the bifurcation parameter. Thick black lines show the fixed points as a function of 0*ii*. For small values of 0*ii* there is a single stable fixed point. As 0*ii* is increased this fixed point loses stability by a super-critical Hopf bifurcation, thus signaling the onset of limit cycle activity. Periodic continuations of this low amplitude activity reveals that it, and a stable fixed point, co-exist with a higher amplitude limit

cycle, thus suggesting that an activity dependent hysteresis drives the system between a fixed point and a high amplitude oscillation, thus giving rise to the bursting activity observed. This can be better seen by superimposing on this diagram the trajectory of a single burst (thin solid black line). Here we can see that the burst terminates through a fold-limit cycle bifurcation. At this stage it is not clear what bifurcation accounts for the emergence of the burst.

A well described feature of anesthetic action is the reduction in cerebral blood flow and metabolism (Kaisti et al., 2003). Therefore during anesthetic action it would be reasonable to assume that the recovery of pre-synaptic neurotransmitter levels will be impaired. In particular as the anesthetic level increases then the rate of synaptic recovery should decrease. **Figure 5** shows the effects of systematic reductions in the synaptic recovery rate for the SS1 model. As θ*i* (synaptic recovery) is decreased the burst duration decreases and the period of non-oscillatory isoelectric activity increases until the model EEG becomes fully isoelectric. This trend is also observed clinically during anesthesia.

Short term synaptic depression would be expected to affect all synapses, though the depression would not be expected to be uniform. So far we have assumed that the synaptic depression would principally affect inhibitory synapses between inhibitory neurons. Will such bursting survive when all types of "fast" synaptic activity is subject to the biological forces of short term synaptic depression? **Figure 6** reveals that bursting does occur when both excitatory and inhibitory synaptic activity undergoes activity dependent short term synaptic depression. Further, the bursting that emerges is strongly modulated by parameters of the "fast" **Table 3 | Model parameter set used in simulations of Figures 3–8.**

Parameter set taken from Bojak and Liley (2005). Because parameters were obtained as part of a numerical search their full precision had been detailed, however their sensitivity to perturbation is much less than the precision reported (Bojak and Liley, 2005).

EEG system (see **Table 2**) that have been identified as targets for the action of anesthetic agents. Prolonging the decay of the unitary IPSP and reducing subcortical input are both found to significantly modulate modeled EEG bursting. In particular it is found that reducing excitatory extra-cortical input (*pee*), which

**FIGURE 4 | One-dimensional continuation of the "fast" EEG system in he suggests how we might account for the dynamics of bursting.** Thick black lines (solid = stable, dashed = unstable) shows the fixed points (fp) as a function of 0ii. Solid black circle indicates fixed point for the default parameter values for "fast" EEG system (parameters as for **Figure 3**). Solid blue lines indicate the locus of maximum and minimum amplitudes respectively of stable limit cycle (lc) activity, whereas solid red lines correspond similarly to locus of unstable limit activity. Transitions between stable and unstable lc activity is predominantly through fold-limit cycle bifurcations, except for middle limit cycle branch where transition is through a torus bifurcation. Superimposed on this one-dimensional bifurcation diagram is the trajectory of a single burst.

presumably dominantly arises from thalamus, leads to very long quiescent (isoelectric) periods.

One of the limitations in using SS1 is that the PSPs are, in principle, free to take any value whereas physiology would dictate that they should remain bounded. To explore the effects of this restriction we chose to define an alternative slow modulating system [SS2, equation (12)]. **Figure 7** shows the effect of utilizing this system to provide a slow activity dependent modulation of excitatory and inhibitory synaptic efficacy of our "fast" EEG system. In the left panel of this figure we have plotted *ce* ≡ 0*el*/0 0 *el* as a function of time. The interesting thing to note is that in contrast to **Figure 3** excitatory synaptic efficacy decreases during the quiescent interburst period and increases during the burst. The left hand panel however shows that there is a significant phase difference between the normalized excitatory synaptic efficacy *ce* ≡ 0*el*/0 0 *el* and the

normalized inhibitory synaptic efficacy *ci* ≡ 0*il*/0 0 *il* . This suggests that there is an important dynamical interplay between excitatory and inhibitory synaptic efficacy to regulate neuronal population excitability such that bursting occurs. This may explain why there is confusion in the empirical literature regarding the role alterations in synaptic efficiency have in the genesis of BS.

solutions will undergo a bifurcation to a non-oscillatory state. All parameters

An important difference between SS2 and SS1 is that parameters hypothesized to be important targets of anesthetic action *pee* and ε*il* are able not only to parametrically regulate bursting but appear also able to switch bursting on (presumably through a bifurcation from a stable fixed point). **Figure 8** illustrates this. If the inhibitory neuronal IPSP decay time is not long enough then a single fixed point dominates which has an associated white noise fluctuation spectrum. But as the IPSP decay time increases (beyond ε*ii* > 1.8 for the parameter set chosen) then bursting emerges. However if *pee* is decreased, as we would expect during anesthesia, then the isoelectric period is prolonged until at some critical value of *pee* bursting is extinguished to be replaced by an infinitely long quiescent/isoelectric period.

# **4. DISCUSSION**

as per **Figure 4**.

We have described here how a well-known model of the "fast" dynamics of the EEG can be modulated by a number of slow systems to produce bursting activity that bears some resemblance to BS seen clinically. The slow systems that we used were all based on some form of activity dependent short term plasticity that has been empirically observed, and used successfully in other models of macroscopic level bursting (Tabak and Rinzel, 2005). While we were able to clearly show the existence of bursting, because we did not include any additive or multiplicative noise sources, we were unable to account for the quasi-periodicity of BS. Thus all our busting arises from purely deterministic processes, presumably involving a range of well described bifurcations (Izhikevich,

**FIGURE 6 | The effects of short term synaptic depression of both excitatory and inhibitory cortical synapses in the genesis of burst suppression and its modulation by variations in extra-cortical input (pee)** **and the IPSP decay time (**ε**ii,** ε**ie) (see section 2.3 for further details)**. Parameters: θ e = 0.1818 mV, θi = 0.07 mV, ke = 14 mV s. All other parameters as for **Figure 3**.

2007). However because our system clearly exhibited bistability (see **Figure 4**) it is almost certain that our system will be able to exhibit some form of burst excitability in response to stochastic forcing. Such burst excitability has been described experimentally. For example during BS induced by various halogenated anesthetic agents, bursts can be readily evoked by auditory, visual, or somatosensory stimuli (Hartikainen et al., 1995). Nevertheless while burst onset and duration may be random variables it would seem that the bursts themselves should reveal a high degree of determinism (weak non-linearity) when compared to EEG in which bursting or epileptiform activity is not evident.

Because the parameter space of the underlying "fast" EEG model is potentially extremely large it is not possible to systematically explore its dynamical repertoire and it may be possible that this system, not augmented with one of the slow systems described, is able to burst. Nevertheless, on the basis of our results, and what is known regarding the dynamical mechanisms of bursting, it would seem likely that multiple pathways to BS exist through a variety of activity dependent slow modulatory systems.

Further we might hypothesize that such slow modulatory systems might span a number of functional scales in the brain. **Figure 9** diagrammatically illustrates some possible candidate systems. An obvious activity driven slow modulatory system would be that associated with thalamus and the corresponding thalamocortical feed-back. Mean field models of the EEG that incorporate thalamo-cortical feed-back have been developed (Rennie et al., 2002) and it will be interesting to see if they are structurally configured to support BS. In addition to synaptic fatigue another obvious cortical system that might be marshaled to provide slow activity dependent modulation is the cortico-cortical conduction system. Although at this stage there is currently little evidence to suggest anesthetics slow conduction velocities, it is widely documented that axonal conduction velocities are significantly decreased in

hypothermia. Decreases of up to 5% per ˚C for conduction velocity have been reported (Waxman, 1980). As BS has been observed to occur in hypothermia (Schwartz et al., 1989;Akrawi et al., 1996) we can conclude that a slow system emerging in the long-range

## **REFERENCES**

Akrawi, W. P., Drummond, J. C., Kalkman, C. J., and Patel, P. M. (1996). A comparison of the electrophysiologic characteristics of eeg burstsuppression as produced by isoflurane, thiopental, etomidate, and propofol.*J. Neurosurg. Anesthesiol.* 8, 40–46.

coupling via a slowing of axonal conduction velocity is a possible route to BS.

Our attempts to account for the dynamical pathogenesis of burst suppression differ from other approaches, most notably Ching et al. (2012). Ching et al. described the scalp EEG in terms of the activity of a small scale, biophysically detailed, computational model of interacting populations of cortical and thalamic neurons. Burst suppression was modeled as arising from metabolically induced alterations in an ATP-gated slow neuronal membrane potassium current (*IKATP* ) on the basis that the reduction in cerebral metabolic rate (CMRO2) induced by anesthetic agents and hypoxia was associated with the depletion of ATP, and hence membrane hyperpolarization. While on this basis they claim to have accounted for a number of defining features of BS that included (i) the spatial synchrony of burst onset (ii) the parametric variability of burst duration/isoelectricity and (iii) the characteristically long time scales associated with bursting/isoelectricity compared to resting EEG, some caution needs to be exercised.

Firstly their model of resting/spontaneous EEG is constructed on the basis of the activity of no more than 20 model neurons. Because EEG is a distributed large scale phenomena such a model is unlikely to meaningfully account for resting/spontaneous activity particularly given the absence of any long-range excitatory cortico-cortical coupling. This has important implications for the propagation of burst activity particularly given that the onset of bursts, when examined at fine temporal scales, is probably not truly spatially homogeneous.

Secondly while the relationship between CMRO2 and ATP production cannot be reasonably disputed, not all anesthetic agents that produce reductions in CMRO2 produce BS. For example the noble gas xenon has been reported to reduce CMRO2 by up to 33% in human participants (Rex et al., 2006) yet is not associated with any discernable BS.

Thirdly the approach they have taken to producing BS essentially depends on the slow modulation of a faster system, the approach we have adopted here. The modeled time scales of *IKATP* variability are very long, of the order of tens of seconds.

For clarity and tractability the current investigations have focused on the spatially homogeneous case for the model of Liley et al. [i.e., 52 = 0 in equation (6)]. Clearly the emergence of BS in the spatially extended case will need to be investigated through the appropriate numerical solution of the defining partial differential equations. Because the cortical phase synchrony (Hartikainen et al., 1995) of burst suppression has not, as far as we are aware, been explicitly investigated it will be crucial to empirically determine the spatiotemporal emergence of bursts in order to assess the importance of excitatory cortico-cortical connectivity in the emergence and modulation of BS as implied by the model and as we have suggested. This will require recording high density EEG during anesthesia in which BS is present.

of burst-suppression activity. *Neurophysiol. Clin.* 26, 329–334.

of consciousness. *Neurol. Clin.* 3, 615–631.

V. (1995). Cortical reactivity during isoflurane burst-suppression anesthesia. *Anesth. Analg.* 81, 1223–1228.

*the Nervous System*, eds S. Coombes and P. C. Bressloff (Singapore:World Scientific), 273–301.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 06 December 2012; accepted: 05 April 2013; published online: 30 April 2013.*

*Citation: Liley DTJ and Walsh M (2013) The mesoscopic modeling of burst suppression during anesthesia. Front. Comput. Neurosci. 7:46. doi: 10.3389/fncom.2013.00046*

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

# Neural masses and fields in dynamic causal modeling

# *Rosalyn Moran1,2,3\*†, Dimitris A. Pinotsis 1† and Karl Friston1*

*1 Wellcome Trust Centre for Neuroimaging, Institute of Neurology, University College London, London, UK*

*2 Virginia Tech Carilion Research Institute, Virginia Tech, Roanoke, VA, USA*

*3 Bradley Department of Electrical and Computer Engineering, Virginia Tech, Blacksburg, VA, USA*

#### *Edited by:*

*Peter Robinson, University of Sydney, Australia*

#### *Reviewed by:*

*Peter Robinson, University of Sydney, Australia James Roberts, Queensland Institute of Medical Research, Australia*

#### *\*Correspondence:*

*Rosalyn Moran, Virginia Tech Carilion Research Institute, Virginia Tech, 2 Riverside Drive, Roanoke, 24016 VA, USA e-mail<EMAIL_ADDRESS>
*†These authors have contributed equally to this work.*

Dynamic causal modeling (DCM) provides a framework for the analysis of effective connectivity among neuronal subpopulations that subtend invasive (electrocorticograms and local field potentials) and non-invasive (electroencephalography and magnetoencephalography) electrophysiological responses. This paper reviews the suite of neuronal population models including neural masses, fields and conductance-based models that are used in DCM. These models are expressed in terms of sets of differential equations that allow one to model the synaptic underpinnings of connectivity. We describe early developments using neural mass models, where convolution-based dynamics are used to generate responses in laminar-specific populations of excitatory and inhibitory cells. We show that these models, though resting on only two simple transforms, can recapitulate the characteristics of both evoked and spectral responses observed empirically. Using an identical neuronal architecture, we show that a set of conductance based models—that consider the dynamics of specific ion-channels—present a richer space of responses; owing to non-linear interactions between conductances and membrane potentials. We propose that conductance-based models may be more appropriate when spectra present with multiple resonances. Finally, we outline a third class of models, where each neuronal subpopulation is treated as a field; in other words, as a manifold on the cortical surface. By explicitly accounting for the spatial propagation of cortical activity through partial differential equations (PDEs), we show that the topology of connectivity—through local lateral interactions among cortical layers—may be inferred, even in the absence of spatially resolved data. We also show that these models allow for a detailed analysis of structure–function relationships in the cortex. Our review highlights the relationship among these models and how the hypothesis asked of empirical data suggests an appropriate model class.

**Keywords: dynamic causal modeling, electroencephalography, magnetoencephalography (MEG), local field potential (LFP), neural mass models**

# **INTRODUCTION**

Over the past two decades, BOLD neuroimaging techniques have been successfully applied in human studies to identify regions of functional specialization, to within a scale of a few millimeters (Ashburner, 2012). Electrophysiological recordings including magneto- and electro-encephalography (M/EEG) offer an aggregate measure of neuronal activity (in the order of hundreds of thousands of neurons) but at a millisecond timescale (Baillet et al., 2001; Nunez and Srinivasan, 2006). Though localizing activity is mathematically ill-posed in these electromagnetic modalities, the wealth of spatial information from fMRI studies can now support M/EEG as a powerful modality for the analysis of functional integration in the human brain.

Dynamic causal modeling (DCM) is designed to probe the mechanisms of effective connectivity (the influence of one brain region on another) that underlie multi-region network responses in neuroimaging (fMRI, M/EEG) data. The approach uses a neurobiologically motivated model which is inverted or fitted to empirical observations using Bayesian techniques (Daunizeau et al., 2011). These comprise separate generative processes at the neuronal level and at the observation level. For M/EEG, neural mass and neural field models in particular are used, to support this analysis by quantifying the temporal and spatiotemporal evolution of macroscopic brain activity, using physiologically plausible dynamics. In DCM; as implemented in the SPM software (Litvak et al., 2011), neural mass and field models are used as generative models to infer the synaptic parameters and effective connectivity that constitute active brain networks. More recently, these models have been applied in (single-region) DCMs as a "mathematical microscope"—to test synaptic hypotheses at the level of specific laminae and receptors (Moran et al., 2011b), and disambiguate between structural and functional hypotheses; for example, explaining intersubject variations in gamma oscillations (Pinotsis et al., 2013b). In this review, we summarize the state-of-the-art in modeling such population-based activity and demonstrate their use in the context of DCM. The development of these modeling approaches has been underpinned by pioneering developments several decades ago that produced generative models of EEG data characteristics based on neural masses (Wilson and Cowan, 1972; Nunez, 1974; Freeman, 1975, 1987; Jansen and Rit, 1995; Valdes et al., 1999; Wendling et al., 2000).  The types of models we will review here have played a direct implementational role in DCM and do not represent an exhaustive overview—a non DCM-centric treatment is provided in Deco et al. (2008). Here we consider, in particular, how DCM treats ensemble neuronal activity as a point process (neural mass models) or explicitly incorporate a spatial dimension (neural field models). Both types describe so-called "mesoscopic" properties of neural activity, employing statistical mechanics to transform single unit activity into population activity—where appropriate composites can be used to generate macroscopic data. Since this mesoscale is hidden from direct observation, we demonstrate how these models rely upon and exploit knowledge about synaptic and cell physiology, as well as neuroanatomy.

Two distinct biological perspectives have informed the development of neural models, leading to a taxonomy of "convolution" or "conductance-based" models. These distinctions arose from the consideration of cortical mesocolumns—convolution models (Freeman, 1975); and separately from the consideration of a single cell's electrophysiological properties—conductance models (Hodgkin and Huxley, 1952). Early work by Wilson and Cowan (1973) derived a sigmoidal relationship for transforming population membrane potential to an average population firing rate. These models consisted of sigmoidal and convolution-based operators and were refined on the basis of empirical observations by Freeman, Wendling and others (Freeman, 1987; Jansen and Rit, 1995; Wendling et al., 2000). In contrast, conductance based models were formulated as an equivalent circuit model of an excitable cell membrane: Hodgkin–Huxley's original description of the giant squid axon is the classical example of this sort of model—and was reduced to a two dimensional form by Morris and Lecar (1981). Their reduced circuit has been scaled up for M/EEG analysis in DCM using the Fokker-Planck formalism (Breakspear et al., 2010) to describe the evolution of population densities (Marreiros et al., 2009). In this setting, when only the first order statistics (e.g., mean) of the population density are considered, the model describes a neural mass (where the population density can be regarded as a point of mass). When higher order statistics are considered, we obtain a "mean-field" model (where the full density of one population depends on the mean of another). We note that the terminology here is rather specific to DCM; the focus of this review. In other settings, and in other treatments of neuronal activity (Coombes, 2010; Buice and Chow, 2013), mean-field models often refer to population dynamics with interacting means only. However, here we use the term *neural-mass* to refer to an interaction in population means and *mean-field* to higherorder interactions to remain consistent with the DCM literature and to acknowledge the early neural-mass nomenclature developed by Valdez-Sosa and other pioneering work in this field (Valdes et al., 1999). Both neural mass and mean field formulations can be applied to convolution and conductance based models: The choice of either convolution or conductance based model depends on the type of inference required (when applying the model to real data), with the latter offering a richer and biologically more realistic parameterization of synaptic currents.

The deployment of neural mass (or mean field) models of populations in DCM entails further neurobiological plausibility, through a laminar specification of cell types and their interconnectivity. For neocortical studies, a laminar architecture is populated with neuronal ensembles, so that forward (e.g., thalamo-cortical), backward or lateral (e.g., inter-hemispheric) extrinsic connections impinge upon pyramidal, spiny stellate or inhibitory interneurons (David et al., 2006). This construction is motivated by tracing studies in the macaque (Felleman and Van Essen, 1991) and demonstrates the first constraint under which these models were developed for DCM. Namely; that they conform to known physiological and anatomical principles. A second constraint is that they must be able to generate stereotypical features of empirical macroscopic measurements; for example, dominant alpha rhythms (David and Friston, 2003) or late potentials in evoked transients (Garrido et al., 2007a). In this sense, none of the models are "right" or "wrong"—but can be usefully compared to test a particular hypothesis (Box, 1976).

In addition to the distinction between neural mass and mean field formulations of either convolution or conductance based models, we also have to consider the distinction between models based upon ordinary differential equations and partial differential equations (PDEs) that endow neuronal populations with spatial attributes: incorporating the spatial domain into DCM was motivated by the advent of spatially resolved population recording modalities (Pinotsis et al., 2012). This use of neural fields, was proposed as a semi-quantitative treatment of electromagnetic brain activity by Jirsa and Haken (1996, 1997) and Robinson (2006). Crucially neural fields enable local axonal arborization to be modeled directly and can generate topological data features. These may be particularly resolved in high-density subdural grid electrodes (electrocorticography) and optical imaging techniques and also contribute to the topographical distribution of sensor/scalp space measurements in M/EEG.

In this review, we hope to provide a didactic treatment of the neural mass and neural field models available in DCM and highlight application studies that exemplify their use. This complements more general treatments of neural population modeling (Deco et al., 2008). The first section considers convolution-based neural mass models. We will demonstrate their use in inferring causal interactions among multiple brain regions and highlight the minimal assumptions needed to form—and test—competing hypotheses. In this section, we will also introduce the important distinction between different models and different data features; noting that the same models can be used for (and indeed should be capable of generating) different data features. We will focus on the distinction between time and frequency domain responses—highlighting the use of identical neural mass models when modeling evoked and steady state responses. In the second section, we examine conductance-based models and how new currents can be added to enhance physiological detail at the synaptic level. We also examine the impact of second-order interactions among neuronal ensembles in mean field models, particularly in the spectral domain. The third section introduces the spatial parameterization in the form of partial differential or neural field equations—and how these have been applied to test alternative explanations for gamma activity in the visual cortex. In this example, we reconsider lateral connections and the role of distinct pyramidal cell populations. This final model recapitulates "a canonical microcircuit" and provides a framework for investigating differences in directed oscillations. This development was motivated by theoretical considerations about message passing in the brain; namely predictive coding and implications for spectral asymmetries in laminar specific communication (Friston, 2005, 2009). Important asymmetries of the sort are evident in several recent empirical observations (Maier et al., 2010; Buffalo et al., 2011; Bastos et al., 2012) on the laminar specificity of oscillation frequencies. In principle, these sorts of observations can be used with DCM, to address key questions about reciprocal message passing in the brain and its hierarchical architecture.

# **CONVOLUTION BASED NEURAL MASS MODEL**

# **GENERATIVE MODELS OF EVOKED RESPONSES: THE ERP MODEL**

In DCM, event related potentials are modeled as the response of a dynamic input–output system to exogenous (experimental) inputs (David et al., 2006; Kiebel et al., 2006; Garrido et al., 2007b). The DCM generates a predicted ERP as the response of a network of coupled sources to sensory (thalamic) input with the form of a narrow (16 ms) Gaussian impulse function that accounts for some temporal smoothing in thalamic volleys. Each source is modeled as a point source (c.f., equivalent current dipole) comprising three subpopulations, each assigned to a particular cortical layer. For simplicity, we place an inhibitory subpopulation in the supragranular layer. This receives inputs from excitatory deep pyramidal cells in an infra-granular layer which are, in turn, driven by excitatory spiny cells in the granular layer; layer IV. These three subpopulations are connected with intrinsic coupling parameters as shown in **Figure 1**. Though these models operate as a point process, by specifying different layers we can call on anatomical rules of extrinsic (region to region) connectivity (Felleman and Van Essen, 1991). Specifically,

link subpopulations in each source. Neuronal states include currents, *i*, and membrane potentials *v*. Extrinsic connections enter at specific cortical

(*i*) synapses.

average post-synaptic depolarizations (*He/i*) at excitatory (*e*) and inhibitory

the "ERP" and "LFP" convolution based models can be assembled within a generative network using three distinct types of connections. Forward connections correspond to afferent pyramidal axons and synapse on layer IV stellate cells, while backward afferents impinge upon pyramidal and inhibitory interneurons outside of layer IV. Lateral, inter-hemispheric connections are modeled with a postsynaptic response that is elicited in all layers.

The evolution of neuronal activity in this anatomical architecture is controlled by two simple operations following Jansen and Rit (1995). The first is a convolution operation (**Figure 1**) that lends the model its name and models the average membrane depolarization response as a low-pass impulse response. This transforms the average density of pre-synaptic firing arriving at the population into the average postsynaptic membrane potential (PSP). This response is specified by two biologically informed parameters; one tunes the maximum amplitude of PSPs and represents the receptor density and the second is a lumped representation of the sum of synaptic rate constants (of passive membrane and other spatially distributed delays in the dendritic tree). The output operator (**Figure 1**) then transforms this average membrane potential into the average rate of action potentials fired by the population. This transformation is assumed to be instantaneous and is described by a sigmoid function with parameters that determine its shape and location. These parameters model the voltage sensitivity or gain of the subpopulation and its average threshold. It is this function that endows the model with non-linear behaviors that are crucial for phenomena like phase-resetting of the M/EEG. The sigmoid form for these activation functions was originally motivated as arising from a unimodal distribution of threshold potentials within a population of Heaviside response units (Wilson and Cowan, 1973). More recent formulations that connect directly to full mean field (population density) treatments consider the sigmoid form to arise from the distribution of depolarizations within a population, under a fixed threshold. This output of this sigmoid function (presynaptic input) is scaled by intrinsic and extrinsic connectivity parameters from subpopulations within a source or from pyramidal cell afferents that arise from other sources in the network. Thus activity promulgates and reverberates throughout the network (**Figure 1**). Delays along these connections are also parameterized with values that correspond to the time taken for axonal propagation between layers (∼2 ms) and regions (∼16 ms).

A pair of ordinary differential equations completely describes the dynamics of each subpopulation within a source (**Figure 1**). These are deterministic and—for a DCM of ERPs—a thalamic impulse timed to correspond to some experimental stimulus perturbs the sources. The spatial arrangement of pyramidal cell dendrites (perpendicular to the cortical surface) renders them the prominent sources of measurable electromagnetic signals and are thus harvested from each source in the network and passed through a lead field, to produce the spatiotemporal patterns observed in M/EEG sensors.

Early applications of this model led to a series of EEG-based investigations into oddball effects by Garrido et al. (2007a, 2008, 2009). The mismatch negativity (MMN) is a negative change in the auditory evoked potential that occurs after an unpredictable change in the acoustic stream. It is elicited, for example, when deviant frequency tones are embedded in a stream of repeated tones. The MMN has a fronto-temporal topology and is thought to reflect the updating of an internal model of the sensorium, where by a sensory prediction error is registered and a new prediction formed (Näätänen et al., 2005). Competing "ERP" neural mass models were compared (in terms of their model evidence, using standard Bayesian techniques) to probe the type of extrinsic connection changes that mediate the MMN. Indicating both a bottom-up sensory prediction error and a topdown change in predictions, a stimulus specific modulation of both forward and backward connections among hierarchically deployed sources exhibited trial specific (deviant compared to standard) modulation (Garrido et al., 2007b). This network and paradigm has since been investigated in pathological settings. A striking example by Boly et al. (2011) used the same model comparison procedure to distinguish between top-down and bottom-up extrinsic connections and their changes with levels of consciousness.

# **GENERATIVE MODELS OF EVOKED AND SPECTRAL RESPONSES: THE LFP MODEL**

A second convolution-based model, named the "LFP" model was developed from the "ERP" model. The model was augmented to address the neurotransmitter basis of changes in intracortical local field potentials from rat prefrontal cortex (Moran et al., 2007)—and now also serves as a generative model for noninvasive EEG and MEG studies (Boly et al., 2012). The model differed from the "ERP" model in two ways: first, based on biophysical models by Whittington et al. (1995) of gamma oscillations in the hippocampus, the role of inhibitory interneurons was augmented with recurrent self-connections (**Figure 1**). This subtle addition was important from the perspective of a new set of questions. Here our goal was to develop a generative model of spectral responses can exhibit high-frequency oscillations. A second extension refined the neurophysiological input-output transforms; whereby spike-rate adaptation was modeled at the input stellate cell population. This involves the addition of currents based on the phenomenological model described by Benda and Herz (2003) and combined several ionic currents modulating spike generation—including voltage-gated potassium currents (M-type currents), the interplay of calcium currents and intracellular calcium dynamics with calcium-gated potassium channels (AHP-type currents) and the slow recovery from inactivation of the fast sodium current.

Practically these additions lead to a larger dynamic repertoire using an identical connectivity architecture (three subpopulations within a source and extrinsic forward, backward and lateral connections with laminar specificity) as that described above. The focus of the LFP model was the reproduction of fast synchronous activity as summarized with the steady-state spectral density (Fourier transform) of time series data. In order to generate steady-state spectral responses, we linearized the model's differential equations around an equilibrium point. This equilibrium or operating point is obtained by integrating the system over a protracted time window. The linearity assumption will

accommodate parameter spaces in the region of fixed points and local bifurcations (Friston et al., 2012), known to emerge from this sort of model (Grimbert and Faugeras, 2006). In DCM, the neural masses are treated as a system that is perturbed by white and pink noise; which is explicitly parameterized. This provides a compact summary of the system, where the system's spectral responses can be obtained from its transfer functions, which depend on the physiological parameters of the model (and neural noise) (Nunez, 1974; Steyn-Ross et al., 1999; Robinson et al., 2001; Moran et al., 2007). In other words, the transfer function links unobserved physiological processes to measured spectral responses and is an essential part of forward or generative models of spectral measures. This "LFP" neural mass model was first used in a single-region DCM analysis to demonstrate how one can make inferences about synaptic function at the neuronal level, using macroscopic electrophysiological measurements. This proof of principle used microdialysis measures of extracellular neurotransmitter for validation (Moran et al., 2008) and extended the breadth of applications of DCM—in this case by inferring condition specific modulations of synaptic parameters.

# **CONDUCTANCE BASED NEURAL MASS MODELS THE NMM AND MFM**