PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS One. 2010; 5(9): e12695.
Published online 2010 September 22. doi:  10.1371/journal.pone.0012695
PMCID: PMC2943909

Mechanisms Explaining Transitions between Tonic and Phasic Firing in Neuronal Populations as Predicted by a Low Dimensional Firing Rate Model

Paul L. Gribble, Editor

Abstract

Several firing patterns experimentally observed in neural populations have been successfully correlated to animal behavior. Population bursting, hereby regarded as a period of high firing rate followed by a period of quiescence, is typically observed in groups of neurons during behavior. Biophysical membrane-potential models of single cell bursting involve at least three equations. Extending such models to study the collective behavior of neural populations involves thousands of equations and can be very expensive computationally. For this reason, low dimensional population models that capture biophysical aspects of networks are needed. The present paper uses a firing-rate model to study mechanisms that trigger and stop transitions between tonic and phasic population firing. These mechanisms are captured through a two-dimensional system, which can potentially be extended to include interactions between different areas of the nervous system with a small number of equations. The typical behavior of midbrain dopaminergic neurons in the rodent is used as an example to illustrate and interpret our results. The model presented here can be used as a building block to study interactions between networks of neurons. This theoretical approach may help contextualize and understand the factors involved in regulating burst firing in populations and how it may modulate distinct aspects of behavior.

Introduction

Different populations of cells in the nervous system of many organisms display sudden, organized, and collective changes in spiking activity. Such changes in population firing involve possibly many thousands of cells. A population burst occurs when the population firing rate suddenly increases and then goes back to the basal rate. Population bursts are produced during normal behavior, but also in pathological situations [1] and are displayed in a variety of central regions of the nervous system in vertebrates (e.g., midbrain, thalamus, subiculum, hippocampus, olfactory bulb, and spinal cord) and invertebrates (ventral cord and antennal lobe in insects, stomagogastric ganglia in lobster and other crustaceans). In addition, population bursts are believed to underlie different aspects of normal and pathological function [2] in the nervous system. For instance, periodic bursting in the respiratory groups of the mammalian brainstem occurs at fixed phase lags [3], [4]. These oscillations in population firing are also present in networks of motor neurons that control locomotion and other rhythmic activities [5], [6]. Oscillations in population activity are also important in sensory processing. For instance, olfactory projection neurons in the antennal lobe of many insects such as moths [7], flies [8], locust [9], and honeybees [10] display short-lasting responses to short-lasting olfactory stimuli. The different populations involved in these olfactory responses also display oscillatory firing for long-lasting stimuli [11], [12]. Population bursts are also believed to contribute to processes related to learning and memory. For instance, pyramidal cell bursts in the hippocampus are believed to underlie the initial representation and further transference of memory traces from short term to long term storage [13], [14].

There has been a considerable search for methods to appropriately study population activity, especially among neurocomputation studies related to perceptual decision making [15][18], central pattern generator [19][21] and synchronization [15]. The focus of this paper is to construct a computationally efficient model to study macroscopic biophysical mechanisms underlying transitions between different kinds of population firing. The model presented here was created with the idea of studying large circuits formed by different regions of the nervous system (e.g. the hippocampus-nucleus accumbens-pallidum-VTA loop [22]). One requirement for the construction of the model was that the same general formulation should be used as a template to model different populations of neurons, perhaps only differing in the choice of parameter spaces.

A variety of the currently existing network models are based on single cell activity. Some of these models include phenomenological population density formulations based on integrate and fire neurons [23][26], Poisson processes [27], and generalized linear point processes [28]. Among other limitations, these models do not include possibly important dependencies on physiologically relevant phenomena such as different sources of input with different time scales for excitation or inhibition. In comparison, biophysical single cell models require either two or three equations [29] (but see [30] for an interesting hybrid approach) and typically at least 4 parameters per ionic current. That is, biophysical single cell models are often too complex to be directly used as building blocks for a larger neural network. One problem is that the number of equations in a network model with biophysical cells is at least two times the number of neurons, but possibly much larger. Another potential problem is that the dependence of the model dynamics on the parameters can become intractable depending on the level of heterogeneity of the cells in the model. Examples of network simulations based of several biophysical point neurons or complex multiple compartment neurons can be found, respectively in [31], [32], or [33]. The study of small networks of synaptically coupled cells is thus computationally expensive when the size of a network grows to a few thousand cells, even if homogeneity in the parameters is assumed.

For these reasons, we decided to construct a macroscopic model of population activity such that each of the parameters of the model represents an experimentally measurable quantity. That is, we required the model to be macroscopic but biophysical. Importantly, the model is flexible enough to potentially represent several distinct neural populations with the same general formulation. The general formulation used here can be potentially used as a building block for the study of large collections of interacting populations, thereby capturing interactions between different areas of the nervous system with a small number of equations.

Our approach

We view the factors that bear upon cell populations and produce collective increases or decreases in the population firing rates as similar to those that produce spiking in a single cell. In single cells, ions that cross the membrane produce currents that change the membrane potential. Some ionic currents increase the membrane potential, whereas others contribute with a negative feedback that repolarizes the membrane. Analogously, there are factors that contribute to increase/decrease the firing rate in a population of cells. The above analogy can be observed in the firing rate model of Wilson and Cowan [34]. In addition, the phase space of the Wilson-Cowan model resembles the phase space of dynamical systems that describe the spiking activity of excitable cells [29], [35]. In mathematical terms, such a similarity suggests topological equivalence between dynamical systems that represent single cells and populations. In view of the above remarks, we hypothesized that it should be possible to construct population models with similar trajectories and overall qualitative behaviors as single cell models, based on the premise that activity is determined by two processes: a fast one described by an amplifying variable (positive feedback), and a slower one, represented by recovery variable (negative feedback), as is the case for excitable membranes. The analog of an action potential in a population model would be a population burst. Sustained spiking in a single cell model would correspond to a sustained oscillation in firing rate in the population model. If a population displays sustained oscillations with a minimum rate close to the basal population rate then the oscillation can be regarded as periodic population bursting (e.g. locomotion networks).

Our model of population activity (in an neural population which we will call An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e001.jpg) is motivated by the above analogy. The parameters of the model can be directly related to macroscopic biophysical aspects of a tissue of choice. The model is capable of reproducing tonic firing, and fast, nonlinear transitions from low, to high, back to low firing that resemble the excitability (spikes) in single cells. These tonic-phasic-tonic transitions are the population bursting described above, and their periodicity can be changed by varying different parameters. These mechanisms are captured by a two-dimensional dynamical system (i.e., two differential equations). In particular, drawing analogies from the phase plane analysis of single cell models, our model explains the qualitative transitions in terms of what can be regarded as population excitability.

For a simple visual illustration: Figure 1 shows the correspondence between the membrane potential dynamics of one (typical) neuron in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e002.jpg (represented by the one compartment, single cell minimal Av-Ron bursting model [29], described in Appendix S1) and the dynamics of the population An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e003.jpg as a whole. The transitions between different firing regimes in a single cell are shown in Figure 1A, both as a function of time (left) and in phase space (right). The corresponding time-dependent and phase space representations of the population firing, as simulated by our system, are shown in 1B.

Figure 1
Comparison of single cell transitions between firing and quiescence, and its analog in the population firing-rate based model (6)–(7).

In the rest of this article, and for the purposes of illustration, the parameters of the model are tuned to mimic the activity commonly observed in populations of midbrain dopaminergic neurons (MDNs); the values and ranges used are summarized in Table 1. The transitions between different population firing regimes will be described in terms of the population activity of MDNs.

Table 1
Parameters for the model (6)–(7).

Physiology of midbrain dopaminergic populations

In general terms, in vivo rodent and primate MDNs from the ventral tegmental area (VTA, A10) and susbtantia nigra compacta (SNc, A9) display basal tonic firing rates of up to 20 Hz in vivo [36][40]. Subsets of MDNs also burst at rates of up to 200 Hz in response to novel or partially unpredicted stimuli [41], [42]. Importantly, bursting behavior in the VTA is displayed both at the single cell and population levels [38], [43].

Midbrain dopaminergic cells receive typical synaptic inputs from external sources, which can be net excitatory or inhibitory. Excitatory synaptic input through cholinergic and glutamatergic terminals is received from several sources including the subthalamic nucleus, the peduculo-pointine tegmentum [22] and others [44][46]. Inhibitory GABAergic synapses are activated by cells from within the VTA and SNc and also from basal ganglia nuclei [47] and other sources [48]. Although for more particular distinctions and analyses it would be valuable to differentiate these as separate sources of input, in our model, for the sake of simplicity and generality, we group them together as a net extrinsic input rate.

In addition to fast excitatory and inhibitory chemical synaptic input, there is gap-junctional coupling between dopaminergic cells [33], [38], [49]. Electrical coupling is a widely observed phenomenon potentially significant in the synchronization of neuronal populations in some cases [33], [50], and/or acting as a frequency filter in other situations [49]. For instance, electrical coupling may be responsible for generation and stabilization of burst firing in hippocampal networks [51], [52]. Electrical conductance between neurons in not limited to early brain development, as previously believed. Even though the high number of gap junctions in the immature brain declines rapidly during development [53], [54], it is now known that electrical communication exists even between mature nerve cells. In the midbrain in particular, electrical connectivity decreases from 96% to about 20% [49] in a period of weeks, but the latter degree of connectivity persists throughout adult life. While lacking definitive experimental evidence on the contribution of gap junctions to the regulation of burst firing in MDN, it is generally agreed that it induces similar firing between coupled cells, and thus constitutes a source of internal amplification, or recruitment within MDN populations (see Appendix S1). Other sources of internal amplification may include the inactivation of A-type channels [55], NMDA receptor based excitation [56] and neuromodulator input [57]) (by analogy with thalamocortical population bursting, synchronized mainly via inhibitory thalamic reticular nucleus input [58]). Irrespective of the mechanism of this internal amplification, one of our model's aims is to establish whether it is a necessary factor underlying burst activity in the midbrain.

Animal studies show that firing rate of dopaminergic cells is also modulated via D2-dopaminergic receptor activation [38] [59]. This modulation results in down regulation of the firing in MDNs and can be regarded as a source of autoregulation within an MDN population. Extracellular dopamine is typically present in intervals of up to 300 milliseconds after a burst [43], [60], which means that autoregulation of MDN firing should occur within that window of time in the absence of other influences. Studies have pointed out a strong tie between dopamine autoreceptor dynamics and expectation of reward in rodents [59], and similarly, between midbrain autoregulatory factors and novelty-related traits in humans [61]. Notably, however, midbrain (VTA and SNc) dopamine autoreceptors may be characteristic to the rodent brain, since they have not yet been found in the same regions in humans. This suggests that other dopamine mechanisms may perform in humans this down-modulatory function [62]).

Bursting in MDNs can be triggered through several pathways. For instance, glutamatergic and cholinergic inputs from the pedunculo-pointine tegmentum [63] and laterodorsal tegmentum [47] are known to produce bursting in the VTA. In addition, glutamatergic inputs to nucleus accumbens and other striatal targets increase the tonic firing rate and burst firing in the VTA [64].

These firing rate ranges and sources of synaptic input were taken into consideration to constrain our model. Analysis of the different firing regimes displayed by the model as function of parameters was then conducted. Our simulations explain mechanisms by which different transitions between tonic, bursting, sustained bursting, and high tonic firing occur in our model, in terms of quantities relevant to the MDN system.

Modeling methods

Construction of the model

As argued above, we consider a neuronal population An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e020.jpg whose activity as a whole can be captured by a representative firing rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e021.jpg. The rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e022.jpg can be thought of as a weighted sum of the firing rates from the cells in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e023.jpg, or alternatively, as the firing rate of a prototypical cell in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e024.jpg (see Figure 1 and Appendix S1). It is assumed that An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e025.jpg changes as a function of two factors: (1) the average history of the cell's firing over a short time interval in the immediate past, and (2) the synaptic input and other modulatory influences. More specifically, as explained in the previous section, it is assumed that the input is a function of factors including intrinsic excitation (amplification), intrinsic inhibition (dampening), extrinsic excitation, and extrinsic inhibition. In the case of the MDNs in the rodent VTA, the intrinsic excitation can be thought of as resulting from a combination of gap-junctional coupling and NMDA receptor activation [33], [38], [49], [55], [56]. Intrinsic dampening results from spike frequency adaptation and autoregulation by dopamine [38]. Extrinsic excitation would come from glutamatergic and cholinergic synaptic inputs. Extrinsic inhibition can result mainly from activation of GABAergic synapses from within [65] and outside the VTA [48].

The inputs to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e026.jpg are integrated by means of a response function An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e027.jpg

equation image
(1)

with An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e029.jpg representing the gain of the response function (sec). An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e030.jpg has an increasing sigmoidal shape, with values between 0 and 1. That is, the population response saturates for large enough inputs and tends to zero as the total input rates decrease.

The population response to an incoming input An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e031.jpg will trigger changes in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e032.jpg within some delay An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e033.jpg as follows:

equation image
(2)

with An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e035.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e036.jpg representing, a time-weighted average of the recent firing between An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e037.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e038.jpg and, respectively, the maximum firing rate of the cells in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e039.jpg. The first factor on the right hand side of Eq. (2) accounts for the history of firing, and the second for the response to new inputs. The maximal theoretical output that can be generated in the population in response to integrated input An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e040.jpg is An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e041.jpg, if the population has no history of recent coordinated firing. Since delays between reception of input and response are typically very short (of the order of a few milliseconds), An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e042.jpg can be assumed to be small. As a consequence, the discrete-time formulation in Eq. (2) can be approximated by a continuous-time equation as:

equation image
(3)

In agreement with the analogy between cellular excitability and population excitability discussed above, we introduce a dynamic variable slower than An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e044.jpg representing negative feedback that tends to decrease the firing response to input. We achieve this by letting the population input be a weighted sum An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e045.jpg, where An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e046.jpg is the intrinsic excitation, An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e047.jpg is the net extrinsic excitation term, and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e048.jpg is the intrinsic dampening, with An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e049.jpg acting as a “sliding control” that tends to decrease the firing rate for a given input. The dynamics of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e050.jpg can be written as:

equation image
(4)

with An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e052.jpg representing the time constant of the negative feedback, and with steady state function

equation image
(5)

The parameter An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e054.jpg in Eq. (5) determines the maximum rate of change of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e055.jpg at steady state with respect to the firing rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e056.jpg. The firing rate at which the steady state of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e057.jpg is 1/2 is An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e058.jpg. In single cell models, An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e059.jpg would be similar to the gating of potassium channels that allows repolarization of the membrane during an action potential.

The gate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e060.jpg is not the only parameter which may change over time. However, in the present study the other parameters are kept fixed, and the dynamics of the system

equation image
(6)
equation image
(7)

are studied for a variety of fixed values using bifurcation analysis.

System parameters

Before beginning to analyze the dynamics of the system under perturbation, we note that quantitative choices for parameter values and ranges are also based on experimental evidence from MDN recordings in freely moving rodents.

Integration of inputs

In principle, we allow firing in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e063.jpg to range between zero and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e064.jpg Hz. However, An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e065.jpg is just a theoretical absolute maximum, and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e066.jpg cells spends most of their time at much lower firing rates (in fact, we show in the next section that firing is constrained to remain below an asymptotic limit of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e067.jpg Hz, consistent with the maximum rates found experimentally [39]). According to known distributions of firing rates found in experimental recordings, we considered the “critical” firing rate (at which the cells are “at half capacity,” and most sensitive to input) to be An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e068.jpg Hz. The time-constant for the DA neurons, representing the entire time interval over which input can be integrated, was found to range up to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e069.jpg msec [36]. We estimate the typical response-time of the population (time spent by An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e070.jpg on integrating incoming input before responding with an action potential) to be significantly lower, and in our analysis we considered values from An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e071.jpg from An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e072.jpg msec, to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e073.jpg msec.

External input

The cumulative firing rate of external inputs to MDN can be very large (the midbrain can get net excitatory inputs of the order of 1–2 kHZ). This translates post-synaptically to a much lower rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e074.jpg, on the order of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e075.jpg Hz, which is the effective rate integrated by the MDN to create an EPSP.

Intrinsic amplification

An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e076.jpg can be considered to act as a “coupling index”, that is the fraction of the maximal electrical coupling (i.e., An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e077.jpg if all cells within An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e078.jpg are mutually coupled). Dye coupling studies in midbrain regions observed that a single DA neuron may be coupled with one to five similar cells [33]. Grace and Bunney showed that 2 to 5 neurons out of the 18 SNc cells under study were electronically coupled [38], and Vandecasteele et al. found coupling indices from An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e079.jpg to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e080.jpg at various development ranges [49]. Electrical studies found pair indices from An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e081.jpg to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e082.jpg, and various average coupling conductances during the development [49]. To incorporate all theoretical possibilities, in our analysis we considered all values of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e083.jpg; however, as will be seen, most plausible dynamics occur for An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e084.jpg, which is consistent with the experimental findings.

Internal dampening

Since dopamine seems to be typically present extracellularly up to 300 milliseconds after a burst [43], [60], we considered that the autoregulation of MDN firing by DA autoreceptor occurs within that time-frame (An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e085.jpg). At its maximum capacity, we considered the effect of DA receptor autoregulation to be equivalent to that of a direct inhibitory input of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e086.jpg = 160 Hz. The highest receptor sensitivity An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e087.jpg to extracellular dopamine, corresponding to “half-capacity,” is obtained at the firing rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e088.jpg. Since one of the aims of our model is to study the effects of receptor sensitivity on burst firing, we allowed An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e089.jpg to be a free parameter, varying widely between An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e090.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e091.jpg Hz.

Fine tuning of these parameters does not substantially affect the big picture, which remains qualitatively robust in a neighborhood of the ranges considered. Finally, let us note that, although not all the model parameters are independent (e.g., changing An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e092.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e093.jpg clearly has the same effect), we prefer to keep these quantities separate, in order to preserve their physiological meaning, and we vary only the ones that are suspected to relate to regulation of bursting activity. In the discussion, we further show that many of the dynamic properties of the model, constrained to these parameter values and ranges, agree qualitatively and quantitatively with population firing considerations drawn from electrode recordings.

Basic analysis and bifurcations

Finding the transitions between different qualitative states of a dynamical system generally starts with plotting its nullclines (the curves in the system's phase-plane where either An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e094.jpg or An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e095.jpg), then calculating the equilibria An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e096.jpg of the system (i.e., the points where the two derivatives are simultaneously zero). In fact – as the Av-Ron single compartment cell burster [29], and the original Wilson-Cowan model [34] – our model falls in the class of Fitzhugh-Nagumo systems, which pioneered geometrical analysis of phase-portraits in the context of neurocomputational models [35]. While the general properties of the FitzHugh-Nagumo system have been studied and discussed at length [15], we will hereby focus on the quantitative phase plane characteristics, and on the parameter transitions more particular to our model.

In our model, a first easy remark concerns the region of the plane which contains the relevant dynamics. Indeed, start by noticing that, for any input An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e097.jpg and for any choice of parameters, the response An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e098.jpg. Since An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e099.jpg, this implies that An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e100.jpg for An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e101.jpg and that An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e102.jpg for An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e103.jpg. Hence, if An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e104.jpg, then An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e105.jpg, so that An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e106.jpg when An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e107.jpg. So all trajectories end asymptotically in the open strip An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e108.jpg. In addition to this, one can also note that we always have An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e109.jpg, which in turn implies that An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e110.jpg for An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e111.jpg and that An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e112.jpg for An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e113.jpg, constraining the long term dynamics to the open strip An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e114.jpg. In conclusion: the open region An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e115.jpg traps the asymptotic phase plane dynamics of system, in the sense that all trajectories eventually end in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e116.jpg and that any attracting sets (hence any features relevant to the long term dynamics like equilibria and cycles) are also contained in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e117.jpg.

The high nonlinearity of the system makes exact computation of the equilibria very difficult; however, it is not hard to show that the the monotonicity of the nullclines implies existence of one up to three intersection points (true in general for any FitzHugh-Nagumo system of ODE). The stability of these points changes with their position in the phase plane, so that different set of parameters can deliver different combinations for the stability of equilibria.

While the An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e118.jpg nullcline has sigmoidal shape, so is globally increasing, the piecewise monotonicity of the nullcline An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e119.jpg depends on the choice of parameters. It is its temporary slope reversal which gives rise to the possibility of multiple equilibria, hysteresis, cycles. (Equivalent conditions on the parameters for this to occur are hard to calculate, but Wilson and Cowan, for example, compute a sufficient condition in their original paper [34].) Obtaining any exact analytical information on the position and geometry of the cycles, when they exist, is even more intractable (as one would generally expect for this degree of nonlinearity). The fact that both nullclines are contained in the region An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e120.jpg implies, of course, that any phase plane cycle would have to also be contained in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e121.jpg, which is a useful initial estimate. Figure 2 shows simulation based illustrations of a few phase planes achievable by of our system, with the respective nullclines, equilibria and cycles (which are allowed to coexist in the same phase plane), together with representative trajectories. More thorough illustrations and descriptions of possible phase plane dynamics can be found in Appendixes S2 and S3.

Figure 2
Different firing behaviors and phase plane configurations.

A bifurcation in the dynamics of a system is by definition a (parameter) state of the system where these dynamics exhibit a qualitative change. For example, for a 2-dimensional continuous time system like (6)–(7), such a qualitative transition could be a change in the stability of its invariant sets: equilibria or cycles. A bifurcation diagram is a graph that shows these invariant sets, and possibly their type and stability, as a function of a single parameter (for co-dimension 1) or of two parameters varied simultaneously (for co-dimension 2) [66][67]. The stability of equilibria can be monitored computationally by following the nature (real or complex) and the sign of the the real part of the local Jacobian's eigenvalues; bifurcations characterized by a local change in stability of cycles are generally harder to establish.

In our study, we use bifurcation diagrams of co-dimensions 1 and 2 to assess the qualitative differences in firing rates and transitions between such regimes – resulting from varying the physiological parameters in the model. For illustrating and understanding these bifurcations, we use numerical computations. (The computations of the steady states and phase plane cycles as well as the bifurcation diagrams were performed in Matcont (version 2.4) [68], a bifurcation finding software based on continuation algorithms.) These are discussed throughout the remainder of the paper.

The stability-changing local bifurcations of codimension one that occur in our system are of three types: subcritical Hopf, limit point (or saddle node bifurcation of equilibria) and fold (saddle node bifurcation of cycles). The terminology in the field is not always consistent between authors; throughout the paper, we remained faithful to the system of terms used by the authors of the Matcont software; for clarification, all terms used are defined in the following paragraph. Each type of transitions corresponds to a different “crash”, or “degeneracy” of the respective equilibria or cycles.

A subcritical Hopf bifurcation appears generally when an equilibrium with complex conjugate eigenvalues changes stability (as one parameter of choice is varied), so that its Jacobian's eigenvalues transition from having positive real part (equilibrium is an unstable spiral) to negative real part (stable spiral), through a pure imaginary stage. This transition also creates a repelling cycle around the stable spiral, whose radius increases locally with the distance to the bifurcation point (see Appendix S2). A limit point (or saddle node) bifurcation appears when two equilibria (one stable and one unstable) traveling in the phase-plane collide and disappear, through an intermediate stage of a common half-stable equilibrium (i.e, with one null eigenvalue of the Jacobian). Similarly, in a fold bifurcation, two nested cycles of opposite stabilities collide and disappear, through a bifurcation stage of a common half-stable cycle. For illustrations of successions of such transitions, see Figure 3, or Appendixes S2 and S3.

Figure 3
Examples of two successions of local bifurcations.

When varying simultaneously two parameters, we encounter two types of codimension 2 bifurcations [66]. The presence of the hysteresis in the phase-plane opens the possibility for cusp bifurcations; these occur where two branches of limit point bifurcation curve meet tangentially. For neighboring parameter values, the system has three equilibria which collide and disappear pairwise via limit point bifurcations. A Bogdanov-Takens codimension 2 bifurcation appears generically at the intersection of a limit-point curve, a Hopf curve and a homoclinic curve. For nearby parameter values, the system has two equilibria (exactly one of which is a saddle) which collide and disappear via the limit-point bifurcation. The nonsaddle equilibrium undergoes a subcritical Hopf bifurcation generating an unstable cycle. This cycle degenerates into an orbit homoclinic to the saddle and disappears via a saddle homoclinic bifurcation. Although interesting in the context of the system's dynamics, the parameter ranges seem to force the physiological dynamics to remain away from the Bogdanov-Takens bifurcations, which is why they will not be discussed here in any further detail.

Results

Tonic versus phasic firing

The main focus of this paper is on mechanisms that trigger and stop bursting, and how these mechanisms can be explained by changes in the model parameters. In other words, such mechanisms will be studied by characterizing the bifurcation structure of the system. The geometry of the bursts, their timing and their kinetics can be very different, depending on the parameter values. Different ways in which parameters tune different features of the bursts are discussed below.

Figure 2 illustrates four possible firing regimes achieved by the system (6)–(7) with four different sets of parameters. The time evolution of the firing rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e163.jpg (Figure 2 left column), is paired in each case with the corresponding trajectory in the An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e164.jpg phase-plane (Figure 2, right column). Depending on the parameter values, phase-plane trajectories can exhibit two qualitatively different types of long-term behavior, determined by the presence or absence of two types of attractors. One such attractor is a tonic firing rate that corresponds to a stable equilibrium in the model. The other attractor is an oscillation, which corresponds to periodic transitions between low and high population firing rates; this regime may be regarded as periodic bursting of the cells in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e165.jpg. Figure 2A shows an example in which the system has a global stable equilibrium at An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e166.jpg33.9 Hz. When initialized at any other value, the firing may undergo an upturn, but always returns to the steady tonic rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e167.jpg. For the parameter values in Figure 2A, repeated bursting is not possible. For the parameters in Figure 2B, any trajectory in the An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e168.jpg-plane starts cycling. The firing rate describes an oscillation between An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e169.jpg200 and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e170.jpg0 Hz. These oscillations can be thought of as a bursting regime with no tonic firing. However, tonic firing and bursting are not mutually exclusive in this model. Figures 2C and 2D illustrate the coexistence of tonic firing and bursting. The right panel of Figure 2C shows a stable equilibrium at An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e171.jpg and a large stable cycle. The basins of attraction of these two attractors are separated by an unstable cycle (not shown). For instance, two nearby initial states could be situated, respectively, inside and outside the region enclosed by the unstable limit cycle. The trajectory of the point inside the unstable limit cycle would go toward An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e172.jpg11 Hz. The trajectory of the point outside the unstable cycle goes toward the limit cycle. That is, the system can asymptotically stabilize to either a low tonic firing, or toward bursting. A similar situation is pictured in Figure 2D, except that the stable equilibrium corresponds to high tonic rate of about An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e173.jpg189 Hz and the coexisting oscillations have longer peaks; note that the bursts are longer and the inter-burst intervals are shorter.

Bifurcation diagrams were constructed to investigate the transitions from tonic to oscillatory firing and back. As an example, two bifurcation diagrams for the half-maximal rate of dampening, An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e174.jpg, are compared in Figure 3. The two panels in Figure 3 illustrate how the equilibria of the system change as An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e175.jpg increases between 0 to 200 Hz. In Figure 3A, the amplification parameter was set to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e176.jpg = 0.5. The stable and unstable rates are shown, respectively, in blue and red. In this example, a limit cycle oscillation around the unstable fixed point emerges as An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e177.jpg increases toward 30 Hz and disappears as An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e178.jpg approaches An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e179.jpg140 Hz. The firing rates in the oscillatory regime alternate between An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e180.jpg200 Hz (top dotted lines) and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e181.jpg0 Hz (lower dotted lines). The transition from low tonic firing into the bursting regime takes place through a thin bistability window. In this window, a stable equilibrium and the stable cycle coexist, suggesting that cells in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e182.jpg can end up in either state, based on their recent history of firing. This coexistence has been observed experimentally in midbrain electrode recordings (for example by Grace et al. [69], [70]). The minimum and maximum rates of the repelling cycle that separates the stable equilibrium and the limit cycle are shown in the bottom inset as a red curve. The bistability onsets with the formation of a cycle around the stable state (fold bifurcation), and disappears when the unstable cycle becomes infinitesimally small and renders the equilibrium unstable (subcritical Hopf bifurcation). In the same figure, as An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e183.jpg continues to increase, the bursting eventually transitions sharply into high tonic firing. This transition occurs at the end of another small bistability window (top inset) formed between a subcritical Hopf and a fold bifurcations (see also Appendix S2 for a more detailed illustration of the phase plane transitions). Note that, strictly speaking, Hopf points mark the entry into the purely cycling regime, with no stable steady-state. Since the bistability windows are very small, and since the Hopf points are easily identifiable on the graph of the steady state, for the remainder of this paper we will consider the Hopf points to be the mark for the onset of sustained (periodic) oscillations.

In Figure 3B, we show the analogous bifurcation diagram for a higher intrinsic amplification An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e184.jpg = 0.75. The transition from low-tonic firing to bursting is qualitatively the same as in Figure 3A, i.e., through a fold and a Hopf bifurcation, and a subsequent bistability window between the two. The transition from bursting to high tonic firing is, however, more mathematically complex, and involves the appearance and disappearance of two additional fixed points, through two limit-point bifurcations (see also Appendix S3). Although of a slightly different nature than the one in Figure 3A, this transition produces the same final result, as oscillations cease for the larger rates of activation for intrinsic dampening.

Triggering sustained firing rate oscillations (tonic-bursting regime)

We study first how our population responds to extrinsic excitation with and without the intrinsic amplification, respectively, for An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e185.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e186.jpg (Figure 4A–B). Each curve corresponds to a different value of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e187.jpg. The steady state An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e188.jpg of the population firing rate increases, as the frequency of extrinsic excitatory input, An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e189.jpg, increases from 0 to 200 Hz.

Figure 4
Steady state firing rate as a function of extrinsic input An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e190.jpg.

As shown in Figure 4A, a unique, stable, and increasing steady state rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e211.jpg (as a function of level of excitation An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e212.jpg) occurs when there is no intrinsic dampening. In other words, the system has only one stable state, and no attracting cycles. That is, periodic bursting cannot be obtained if no amplification is present. The steady state rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e213.jpg increases as the dampening onset decreases (i.e., a larger An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e214.jpg rate places An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e215.jpg on a higher curve). That is, the model naturally predicts that the (asymptotic) population firing rate will be higher if dampening starts at higher rates. The values of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e216.jpg for which An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e217.jpg clearly departs from 0 depend on the half-activation of dampening, An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e218.jpg. If onset of dampening occurs at lower rates (red curve in the bottom of Figure 4), the input firing rate required to increase the population firing rate is much larger than if dampening occurs at higher rates.

To trigger oscillations, it is necessary for the system to have a nonzero intrinsic amplification (Figure 4B). However, although necessary, having An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e219.jpg is not sufficient to produce sustained bursting. As shown in Figure 4B, the onset of firing rate oscillations characteristic of periodic bursting also requires a large enough net balance between the extrinsic inputs and the onset of intrinsic dampening. In addition, bursting requires initial build-up to some low level of tonic firing to pass the stability threshold into the cycling regime. Such thresholds, located at the Hopf points, are shown as red stars in Figure 4B.

Modulating phasic firing

In general, the duration of high firing rates of the oscillatory regime can be modulated by intrinsic excitation and dampening (Figures 5 and and6B),6B), and extrinsic input (Figure 6A). First, consider the changes in burst duration as a function of the intrinsic amplification weight An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e220.jpg (Figure 5). The larger An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e221.jpg, the longer the bursts, with interburst intervals relatively constant. The transition from high/low firing rate oscillations into high tonic firing occurs for a sufficiently large increase in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e222.jpg (see also Figures 4B and and8B8B).

Figure 5
Increasing the duration of the high-firing phase of the oscillation with minimal changes to the low-firing phase.
Figure 6
Changing the duty cycle with minimal changes to the oscillation frequency.
Figure 8
Bifurcation diagrams illustrating the stability of the steady-state firing rate as a function of one parameter.

Other effects in the properties of sustained oscillations can be obtained by changing the external input An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e244.jpg or the intrinsic dampening activation rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e245.jpg. Figure 6 shows how the activation of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e246.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e247.jpg can be varied within relatively wide ranges (An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e248.jpg Hz for An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e249.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e250.jpg Hz for An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e251.jpg, respectively) so that the high firing rate portion of the cycle is longer and the quiescence intervals shorter, and without significantly changing the period of the oscillation. This effect is important, as it suggests a mechanism by which populations of bursting cells can regulate their duty cycles (burst duration/cycle duration) without altering the bursting frequency.

The frequency of the sustained oscillations in the model can be also controlled by varying the time constants An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e252.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e253.jpg. These parameters can be used to regulate the length of the burst and inter-burst intervals (Figure 7A), and the symmetry of the rising and falling phases (Figure 7B). One way to think about the effect of increasing An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e254.jpg is that the right hand side of equation (6) will increase for smaller values of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e255.jpg, yielding larger changes in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e256.jpg per unit time during an oscillation. On the other hand, increasing the time constant of the recovery variable, An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e257.jpg, decreases the time-dependent change in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e258.jpg, which results in slower decrease during the oscillation in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e259.jpg. As a consequence, the frequency of oscillations decreases, without significantly altering the duty cycle (i.e. the burst/inter-burst duration ratio).

Figure 7
Changing the oscillation frequency with minimal changes to the duty cycle.

Ending sustained oscillations

Similarly to the transition from equilibrium (tonic firing) into sustained firing rate oscillations, the converse transition (from oscillations to tonic firing) corresponds to a qualitative change in the dynamics of the system. To study these transitions, the attracting states of the system were calculated by varying only one parameter at a time. Two examples of these codimension 1 bifurcations are shown in Figure 8. Panel A illustrates the bifurcation diagrams for the rate of extrinsic input An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e274.jpg, for different values of the half activation of intrinsic dampening, (An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e275.jpg = 80, 100, 120, and 150 Hz, shown, respectively, in blue, purple, orange, and green). For instance, the lowest of the curves (green) in Figure 8A shows a monotonic increase in the steady state firing as a function of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e276.jpg, for a relatively low onset of intrinsic dampening (An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e277.jpg Hz). The transition (Hopf) points of the diagram are shown by the red stars. For contrast, the first curve from left to right (blue) was calculated assuming a higher onset of intrinsic dampening (An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e278.jpg Hz). Note that the curve is not monotonic, and two additional limit points (green stars) appear between the Hopf points. Recall that limit point (or saddle node) bifurcations are the states of the system where, under variation of one parameter, two equilibria collide and disappear, or, conversely, two equilibria emerge. For the blue curve in Figure 8A, corresponding to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e279.jpg = 150 Hz, a saddle and a repeller appear for large enough An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e280.jpg at the first limit point bifurcation (top green star), as new equilibria to a system already having an attracting equilibrium in the low tonic firing range. Under further increase of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e281.jpg, this low tonic equilibrium changes stability through a Hopf bifurcation (bottom red star), then steadily approaches the newly created saddle point, and eventually collides with it and vanishes at the second limit point bifurcation (bottom green star). The surviving (unstable) equilibrium will later undergo another Hopf bifurcation (top red star) and become attracting, causing the seizure-like, high tonic firing observable for very high values of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e282.jpg. The additional two equilibria are unstable, so they do not constitute attracting states for the firing in the system. However, they play a very important role in the transitions necessary to lead from the low tonic to the high tonic firing. This succession of transitions is mathematically more complex than the fold and Hopf pair observed along the lower (monotonic) curves in Figure 8A.

Panel B in Figure 8 shows the bifurcation diagram for the amplification weight An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e305.jpg, using different values of the dampening half activation (An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e306.jpg = 40, 80, and 120 Hz, respectively, in orange, green, and purple). Observe that the steady state firing rate An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e307.jpg is noticeably higher for higher values of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e308.jpg, and also, limit points are present for a larger range of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e309.jpg.

The cessation of bursting doesn't necessarily occur symmetrically with the bursting onset, and may involve a different succession of bifurcations (compared to the equilibrium-to-cycle transition). In principle, there are multiple and qualitatively different mechanisms of transition from bursting back into tonic firing, each corresponding to different parameter regimes. As has been observed in different experimental settings, the system can transition from bursting into either low or high tonic firing.

From oscillations into high tonic firing

Our first phase-plane and bifurcation plots, staring with Figure 2, suggested that cycling is possible for specific parameter interval. Figure 3 showed that there is an interval of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e310.jpg for which the system exhibits cycling, when all other parameters are held fixed. Along each such curve, there are two Hopf points – the approximate marks for the system entering and leaving the bursting mode. As noted earlier, increasing either type of excitation (i.e., increasing An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e311.jpg in Figure 5A or increasing An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e312.jpg in Figure 5B), or decreasing the intrinsic dampening (i.e., increasing An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e313.jpg in Figure 3) leads the system down a path from low tonic firing through an oscillation window, and into a regime of high tonic firing. However, it is not unreasonable to assume that the predicted seizure-like plateau after the end of oscillations (of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e314.jpg Hz) would be impossible to sustain in real neural populations, eventually leading to cellular death [71], [72] (see Discussion Section).

Transitions into low tonic firing

Sustained firing rate oscillations can be efficiently stopped by shifting the onset of dampening towards lower firing rates. One way in which neuron populations may experience this transition could be a decrease in the activation of dampening, shortly after entering the bursting mode. In MDNs, such a decrease could be triggered by a sudden increase in the local concentration of dopamine due to recruitment of cells into bursting mode. Shifting the dampening onset to lower firing rates can result into a regime where bursting is no longer possible (even though the extrinsic inputs that originally produced this bursting may still be present). The effects of this change can be visualized in Figures 8A–B as “pushing the system down” to a lower curve (corresponding to lower An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e315.jpg). If this drop is significant enough, the system's position on the new curve will be to the left of the Hopf point which marks the onset of bursting, in the range where the equilibrium is stable. Hence the system is forced to return to low tonic firing. It is possible that this forcing inhibitory term may be in effect for as long as necessary to prevent further bursting. For example, lowering the excitation An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e316.jpg may allow the return of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e317.jpg to the pre-oscillatory range, without any need for additional protection against bursting. Conversely, raising An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e318.jpg may require even higher inhibition (lowering the curve even more), if the return of bursting is not yet desirable in the system. This suggests a negative feedback modulation through intrinsic dampening, by a continuous readjustment of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e319.jpg according to the current state of the population firing (as further discussed in the last section).

Bifurcation diagrams in a two-dimensional parameter plane were used to study the global changes in the dynamics as two parameters change simultaneously and independently. In the figures in Appendix S4, we plotted the parameter planes An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e320.jpg, An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e321.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e322.jpg. In all three plots, we notice that the Hopf curve (i.e., the curve that contains all the Hopf points for the particular ranges of parameters) delimits a closed region. This is the region where “periodic bursting is possible” in the system. Recall that a Hopf point is the approximate mark for entering the cycling regime when one parameter is changed. Here, the Hopf curve is an approximate boundary for the 2-parameter region that permits periodic oscillations.

Some of these mechanisms can be more easily captured by allowing multiple system parameters to vary at the same time. Examples of the bifurcations displayed by the system as two parameters are varied simultaneously are shown in Appendix S4.

Discussion

The previous paragraphs describe a model that approximates the activity of a homogeneous population of neurons. The model exhibits transitions between bursting and tonic firing that resemble the typical behavior of several populations of cells in the nervous system. We proceed with some comments on features particular to the model, which we place within the context of MDN activity. We finish the discussion with more general remarks about potential applications of the model in studying the collective behavior of different populations of cells in the nervous system.

In systems without electrical coupling or other intrinsic amplification, solitary bursts might appear simply as a perturbation of the initial condition due to increase in external stimulus. For example, in thalamocortical neurons, population bursting during sleep or absence seizures is synchronized mainly via common inhibitory inputs from the reticular nucleus. However, since MDN are known to be partially electrically coupled (An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e323.jpg), the system has the potential for sustained bursting.

We believe that most of the dynamics observed in the MDN occur close to, and on either side of the pair of fold and Hopf bifurcations (in proximity of the low bistability window) This line of thinking is based on our general work hypothesis that the brain, in its attempt to solve the optimization problems faced by all biological systems, has to maintain its function in a dynamic range that provides sensitivity to a variety of different inputs [73]. Based on both theory and experiment, an entire line of research is working to support the hypothesis that the brain self-organizes as a complex system, functioning far from equilibrium, near a critical state [74], [75]. In the particular case of the MDN function, the proximity to the critical window maximized responsiveness in the sense that a slight increase in excitation may render singular or periodic bursts with rates up to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e324.jpg Hz, which then have to be (and can be) readily suppressed when necessary (e.g., when the triggering stimulus loses its novelty content).

Some experimental studies observed typical MDN burst rates of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e325.jpg Hz [42], [76]; in other studies, the heterogeneity of types, localization and behavior of midbrain cells transpired as wider distributions, ranging from (An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e326.jpg Hz) to very high (up to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e327.jpg Hz) [39]. Similarly, some spike counts quote An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e328.jpg spikes per burst [42], with relatively small interspike intervals of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e329.jpg msec [39]; however, more comprehensive descriptions have included higher quotes, and have additionally noticed that glutamate enhances bursting, increasing spike/burst counts to An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e330.jpg, while leaving inter-spike intervals relatively unchanged. Overall, accounts of burst duration have ranged widely from An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e331.jpg msec [39], [42], and inter-burst intervals from An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e332.jpg msec [42], values which are generally accommodated by the range of rhythms encompassed by our model in the neighborhood of the low bistability window.

Intrinsic amplification is necessary, but not sufficient in producing periodic bursting. As for solitary bursts, the transition into oscillations has to be produced by an additional factor, such as increase in the excitatory synaptic input (e.g., the pedunculo-pontine tegmentum – PPTg – gates bursting in MDNs [46]), or a decrease in the dampening (e.g., a reduction in the GABAergic input to VTA neurons has been correlated to bursting activity within the VTA [65]).

The behavior of the system (6)–(7) can be used to test the hypothesis that MDN firing states are independently regulated by two distinct afferent pathways, whose interactions control tonic and phasic activity in MDNs. On one hand, the mesolimbic dopamine system is strongly influenced by the hippocampus [22], [64]. In brief, infusions of NMDA into the ventral subiculum (vSub) increases the number of spontaneously active DA neurons (population activity), while having no effect on firing rate or average bursting activity. In contrast, NMDA activation of the pedunculopontine tegmentum resulted in a significant increase in DA neuron burst firing without significantly affecting population activity. Simultaneous excitation of the vSub and PPTg induced significant increases in population activity and burst firing in MDNs [48] [47]. Quite clearly, each of these inputs is not only necessary for regulating normal midbrain function (e.g., current injection triggers MDN bursting only if an NMDA-receptor agonist is added to the standard saline [77]), but also has a very particular effect on population firing. Regarding amplification, excitation and dampening as factors gating the bursting, our model illustrates the importance of an appropriate modulatory balance to the midbrain, for tuning the numbers of tonic versus bursting MDNs. A more detailed study of the influence of extrinsic inputs to populations of MDNs could be performed by extending the model to include more than one area, by coupling different populations like the one described by (6)–(7). Such coupling could incorporate specific inhibitory/excitatory pathways that are well known in these circuits, but whose functional influence on firing in MDN are still not well understood. The specific interactions between cells participating in feedback loops from/to the midbrain are only speculative at this point. In sum, a potential general attribute of the model presented here is the ability to verify or suggest theoretical alternatives to existing hypotheses about the network dynamics in such loops.

The transition into oscillations in the system (6)–(7) cannot occur if the initial firing rate is low. As shown in Figure 4B, the firing activity has to first build up to a level of tonic firing of An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e333.jpgHz, for bursting to occur, which agrees well with experimental evidence [47]. Once the sustained oscillations are triggered in the model, the duration of the high-frequency phases and those of the intervals between them can be tuned by a variety of factors. In the oscillatory regime, high (low) firing rate episodes get shorter/longer when the system is brought close to any point of transition into low (high) tonic firing. In contrast, high (low) firing rate episodes and get longer when the system is closer to a transition into high (low) tonic firing (e.g., seizure-like plateau).

The oscillations in the system (6)–(7) can be stopped by two different mechanisms corresponding to a forced transition into either low tonic or high tonic firing. As noted earlier, the transition between low-tonic to oscillations resemble the busting experimentally observed in MDN populations. The accumulation of dopamine in the vicinity of bursting cells has been shown to affect MDN bursting [38]. The model also predicts a result partially confirmed by experiment: that directly lowering the excitatory input (e.g., PPT input) or increasing inhibition (e.g, GABAergic tone from basal ganglia nuclei) may also stop oscillations [46]. Some limitations and possible generalizations of the model (6)–(7) are discussed next.

The functional details and parameter ranges used to illustrate the behaviour of the model are specifically applicable to populations of MDNs. Therefore, the model (6)–(7) is applicable to at least midbrain-like physiology and connectivity. However, given a region of interest in the nervous system, equations (6)–(7) could be set up to model different families of cells by choosing different parameter sets, thus enabling the study of larger systems. As noted before, some of the observed qualitative features could be applicable to networks that govern population bursting in other brain regions (e.g., the thalamus [78], or the subiculum [79]).

The relationship between the population activity in An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e334.jpg and its inputs was modeled by means of sigmoidal functions. This approach was chosen to capture the saturation effect that occurs in cells as their excitatory and inhibitory inputs are integrated [80]. The qualitative use of such functions in our context is justified by models and experiments that go back to Cowan, Boltzmann, and Hill [34], [81], [82]. Nevertheless, it is worth noticing that the parameters An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e335.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0012695.e336.jpg are only partially justified from a macroscopic perspective, at least for the MDN populations we discuss here.

The transitions between equilibrium, oscillations and back have been studied here from a static perspective, in the sense that the changes in parameters were imposed externally rather than being regulated through an intrinsic model feedback. By analogy with the single cell models by Av-Ron et al. [29], future work will address this issue by introducing particular bifurcation parameters as time-dependent system variables that may trigger the start and stop of oscillations.

General considerations

The focus of this paper has been placed on studying qualitative changes in the system dynamics that result from perturbations of its parameters, particularly from variations of one parameter at a time, without trying to achieve an overall picture of how these changes combine to determine the global behavior of the system. While the more parameters are considered simultaneously, the harder it is to visually illustrate the results, such a study would be very interesting and relevant to understanding the underlying mechanisms that produce these dynamics. One of the immediate goals of future work is to build upon our existing model of bursting in MDN as the center-piece of an entire network as proposed by the Lodge-Grace experiments [47], with network feedback loops captured by new equations. The theoretical framework used here potentially enables the possibility to produce testable hypotheses about mechanisms underlying the normal and pathological activities of networks involving the hippocampal formation, basal ganglia and midbrain monoaminergic nuclei.

The local and even global phenomena observed in the model (6)–(7) have been pointed out in single cell models [29], [30], [83], [84]. The similarity between the dynamics of population firing described here, and the membrane potential dynamics in single cells, is more substantial than phase plane analogies for occasional parameter values and is currently under rigorous examination. One possible direction to follow in this regard could be to use such phase-portrait analysis of dynamics to contextualize and understand the factors involved in regulating monoaminergic modulation systems, and study the ways in which these modulatory networks affect other networks [56]. There have been some attempts to study the effects of dopaminergic modulation and dysfunction using the rat as an animal model, [47], [70], and also in the clinical setting [85][87]. However, crucial information is yet missing for establishing a sustainable link between the two perspectives. A working translational model such as (6)–(7) may be the tool that would optimally combine clinical and basic research results.

Supporting Information

Appendix S1

(0.18 MB PDF)

Appendix S2

(0.22 MB PDF)

Appendix S3

(0.27 MB PDF)

Appendix S4

(0.15 MB PDF)

Acknowledgments

The author wishes to gratefully thank Marco Herrera for the invaluable consultations, as well as Erin McKiernan and Tom Hazy for all the useful suggestions and stimulating discussions.

Footnotes

Competing Interests: The author has declared that no competing interests exist.

Funding: The author has no support or funding to report.

References

1. Zhang X, Cui S, Wallace A, Hannesson D, Schmued L, et al. Relations between brain pathology and temporal lobe epilepsy. Journal of Neuroscience. 2002;22:6052–6061. [PubMed]
2. Friedman A, Friedman Y, Dremencov E, Yadid G. VTA dopamine neuron bursting is altered in an animal model of depression and corrected by desipramine. Journal of Molecular Neuroscience. 2008;34:201–209. [PubMed]
3. McCrimmon D, Zuperku E, Hayashi F, Dogas Z, Hinrichsen C, et al. Modulation of the synaptic drive to respiratory premotor and motor neurons. Respiration physiology. 1997;110:161–176. [PubMed]
4. Bailey E, Jones C, Reeder J, Fuller D, Fregosi R. Effect of pulmonary stretch receptor feedback and CO2 on upper airway and respiratory pump muscle activity in the rat. The Journal of Physiology. 2001;532:525–534. [PubMed]
5. Sqalli-Houssaini Y, Cazalets J, Clarac F. Oscillatory properties of the central pattern generator for locomotion in neonatal rats. Journal of neurophysiology. 1993;70:803–813. [PubMed]
6. Tierney A, Harris-Warrick R. Physiological role of the transient potassium current in the pyloric circuit of the lobster stomatogastric ganglion. Journal of neurophysiology. 1992;67:599–609. [PubMed]
7. Christensen T, Hildebrand J. Male-specific, sex pheromone-selective projection neurons in the antennal lobes of the mothManduca sexta. Journal of Comparative Physiology A: Neuroethology, Sensory, Neural, and Behavioral Physiology. 1987;160:553–569. [PubMed]
8. Wilson R, Turner G, Laurent G. Transformation of olfactory representations in the Drosophila antennal lobe. Science. 2004;303:366–370. [PubMed]
9. Bazhenov M, Stopfer M, Sejnowski T, Laurent G. Fast odor learning improves reliability of odor responses in the locust antennal lobe. Neuron. 2005;46:483–492. [PMC free article] [PubMed]
10. Krofczik S, Menzel R, Nawrot M. Rapid odor processing in the honeybee antennal lobe network. Frontiers in Computational Neuroscience. 2008;2:9. [PMC free article] [PubMed]
11. Heinbockel T, Kloppenburg P, Hildebrand J. Pheromone-evoked potentials and oscillations in the antennal lobes of the sphinx moth Manduca sexta. Journal of Comparative Physiology A: Neuroethology, Sensory, Neural, and Behavioral Physiology. 1998;182:703–714. [PubMed]
12. Joerges J, Küttner A, Galizia C, Menzel R. Representations of odours and odour mixtures visualized in the honeybee brain. Nature. 1997;387:285–288.
13. Buzsaki G. Two-stage model of memory trace formation: A role for. Neuroscience. 1989;31:551–570. [PubMed]
14. Wilson M, McNaughton B. Reactivation of hippocampal ensemble memories during sleep. Science. 1994;265:676–679. [PubMed]
15. Medvedev G, Kopell N. Synchronization and transient dynamics in the chains of electrically coupled FitzHugh-Nagumo oscillators. SIAM Journal on Applied Mathematics. 2001;61:1762–1801.
16. Wong K, Wang X. A recurrent network mechanism of time integration in perceptual decisions. Journal of Neuroscience. 2006;26:1314–1328. [PubMed]
17. Brown E, Gao J, Holmes P, Bogacz R, Gilzenrat M, et al. Simple neural networks that optimize decisions. International Journal of Bifurcation Chaos in Applied Sciences and Engineering. 2005;15:803–826.
18. Roxin A, Ledberg A. Neurobiological models of two-choice decision making can be reduced to a one-dimensional nonlinear diffusion equation. PLoS Computational Biology. 2008;4:e1000046. [PMC free article] [PubMed]
19. Buchanan J. Neural network simulations of coupled locomotor oscillators in the lamprey spinal cord. Biological Cybernetics. 1992;66:367–374. [PubMed]
20. Williams T. Phase coupling by synaptic spread in chains of coupled neuronal oscillators. Science. 1992;258:662–665. [PubMed]
21. Várkonyi P, Kiemel T, Hoffman K, Cohen A, Holmes P. On the derivation and tuning of phase oscillator models for lamprey central pattern generators. Journal of Computational Neuroscience. 2008;25:245–261. [PubMed]
22. Lisman J, Grace A. The hippocampal-VTA loop: Controlling the entry of information into long-term memory. Neuron. 2005;46:703–713. [PubMed]
23. Knight B. Dynamics of encoding in a population of neurons. Journal of General Physiology. 1972;59:734–766. [PMC free article] [PubMed]
24. Knight B. The relationship between the firing rate of a single neuron and the level of activity in a population of neurons: Experimental evidence for resonant enhancement in the population response. Journal of General Physiology. 1972;59:767–778. [PMC free article] [PubMed]
25. Knight B, Omurtag A, Sirovich L. The approach of a neuron population firing rate to a new equilibrium: an exact theoretical result. Neural Computation. 2000;12:1045–1055. [PubMed]
26. Nykamp D, Tranchina D. A population density approach that facilitates large-scale modeling of neural networks: Analysis and an application to orientation tuning. Journal of Computational Neuroscience. 2000;8:19–50. [PubMed]
27. Câteau H, Reyes A. Relation between single neuron and population spiking statistics and effects on network activity. Physical review letters. 2006;96:58101–58105. [PubMed]
28. Toyoizumi T, Rad K, Paninski L. Mean-field approximations for coupled populations of generalized linear model spiking neurons with Markov refractoriness. Neural computation. 2009;21:1203–1243. [PubMed]
29. Av-Ron E, Parnas H, Segel L. A basic biophysical model for bursting neurons. Biological cybernetics. 1993;69:87–95. [PubMed]
30. Izhikevich E. Dynamical systems in neuroscience: The geometry of excitability and bursting. The MIT press. 2006:267–280.
31. Rinzel J, Ermentrout G. Analysis of neural excitability and oscillations. Methods in neuronal modeling. 1998;2:251–292.
32. Traub R, Miles R, Wong R. Model of the origin of rhythmic population oscillations in the hippocampal slice. Science. 1989;243:1319–1325. [PubMed]
33. Komendantov A, Canavier C. Electrical coupling between model midbrain dopamine neurons: effects on firing pattern and synchrony. Journal of neurophysiology. 2002;87:1526–1541. [PubMed]
34. Wilson H, Cowan J. Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal. 1972;12:1–24. [PubMed]
35. Fitzhugh R. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal. 1961;1:445–466. [PubMed]
36. Grace A, Bunney B. Intracellular and extracellular electrophysiology of nigral dopaminergic neurons. I: Identification and characterization. Neuroscience. 1983;10:301–315. [PubMed]
37. Grace A, Bunney B. Intracellular and extracellular electrophysiology of nigral dopaminergic neurons. II: Action potential generating mechanisms and morphological correlates. Neuroscience. 1983;10:317–331. [PubMed]
38. Grace A, Bunney B. Intracellular and extracellular electrophysiology of nigral dopaminergic neurons. III: Evidence for electrotonic coupling. Neuroscience. 1983;10:333–348. [PubMed]
39. Kiyatkin E, Rebec G. Heterogeneity of ventral tegmental area neurons: single-unit recording and iontophoresis in awake, unrestrained rats. Neuroscience. 1998;85:1285–1309. [PubMed]
40. Schultz W, Apicella P, Ljungberg T. Responses of monkey dopamine neurons to reward and conditioned stimuli during successive steps of learning a delayed response task. Journal of Neuroscience. 1993;13:900–913. [PubMed]
41. Schultz W. Predictive reward signal of dopamine neurons. Journal of neurophysiology. 1998;80:1–27. [PubMed]
42. Hyland B, Reynolds J, Hay J, Perk C, Miller R. Firing modes of midbrain dopamine cells in the freely moving rat. Neuroscience. 2002;114:475–492. [PubMed]
43. Schultz W. Getting formal with dopamine and reward. Neuron. 2002;36:241–263. [PubMed]
44. Omelchenko N, Sesack S. Laterodorsal tegmental projections to identified cell populations in the rat ventral tegmental area. The Journal of Comparative Neurology. 2005;483:217–235. [PubMed]
45. Omelchenko N, Sesack S. Glutamate synaptic inputs to ventral tegmental area neurons in the rat derive primarily from subcortical sources. Neuroscience. 2007;146:1259–1274. [PMC free article] [PubMed]
46. Pan W, Hyland B. Pedunculopontine tegmental nucleus controls conditioned responses of midbrain dopamine neurons in behaving rats. Journal of Neuroscience. 2005;25:4725–4732. [PubMed]
47. Lodge D, Grace A. The laterodorsal tegmentum is essential for burst firing of ventral tegmental area dopamine neurons. Proceedings of the National Academy of Sciences. 2006;103:5167–5172. [PubMed]
48. Floresco S, West A, Ash B, Moore H, Grace A. Afferent modulation of dopamine neuron firing differentially regulates tonic and phasic dopamine transmission. Nature neuroscience. 2003;6:968–973. [PubMed]
49. Vandecasteele M, Glowinski J, Venance L. Electrical synapses between dopaminergic neurons of the substantia nigra pars compacta. Journal of Neuroscience. 2005;25:291–298. [PubMed]
50. Perez Velazquez J, Carlen P. Gap junctions, synchrony and seizures. Trends in Neurosciences. 2000;23:68–74. [PubMed]
51. Perez-Velazquez J, Valiante T, Carlen P. Modulation of gap junctional mechanisms during calcium-free induced field burst activity: a possible role for electrotonic coupling in epileptogenesis. Journal of Neuroscience. 1994;14:4308–4317. [PubMed]
52. Skinner F, Zhang L, Velazquez J, Carlen P. Bursting in inhibitory interneuronal networks: a role for gap-junctional coupling. Journal of neurophysiology. 1999;81:1274–1283. [PubMed]
53. Peinado A, Yuste R, Katz L. Extensive dye coupling between rat neocortical neurons during the period of circuit formation. Neuron. 1993;10:103–114. [PubMed]
54. Rozental R, Morales M, Mehler M, Urban M, Kremer M, et al. Changes in the properties of gap junctions during neuronal differentiation of hippocampal progenitor cells. Journal of Neuroscience. 1998;18:1753–1762. [PubMed]
55. Hoffman D, Magee J, Colbert C, Johnston D. K channel regulation of signal propagation in dendrites of hippocampal pyramidal neurons. Nature. 1997;387:869–875. [PubMed]
56. Zweifel L, Parker J, Lobb C, Rainwater A, Wall V, et al. Disruption of NMDAR-dependent burst firing by dopamine neurons provides selective assessment of phasic dopamine-dependent behavior. Proceedings of the National Academy of Sciences. 2009;106:7281–7288. [PubMed]
57. Pulver S, Griffith L. Spike integration and cellular memory in a rhythmic network from Na+/K+ pump current dynamics. Nature Neuroscience. 2009;13:53–59. [PMC free article] [PubMed]
58. Destexhe A, Sejnowski T. Synchronized oscillations in thalamic networks: insights from modeling studies. Thalamus. 1997;2:331–371.
59. Heinz A, Grace A, Beck A. The intricacies of dopamine neuron modulation. Biological Psychiatry. 2009;65:101–102. [PMC free article] [PubMed]
60. Gonon F. Prolonged and extrasynaptic excitatory action of dopamine mediated by D1 receptors in the rat striatum in vivo. Journal of Neuroscience. 1997;17:5972–5978. [PubMed]
61. Zald D, Cowan R, Riccardi P, Baldwin R, Ansari M, et al. Midbrain dopamine receptor availability is inversely associated with novelty-seeking traits in humans. Journal of Neuroscience. 2008;28:14372–14378. [PMC free article] [PubMed]
62. Meador-Woodruff J, Damask S, Watson S. Differential expression of autoreceptors in the ascending dopamine systems of the human brain. Proceedings of the National Academy of Sciences of the United States of America. 1994;91:8297–8301. [PubMed]
63. Lokwan S, Overton P, Berry M, Clark D. Stimulation of the pedunculopontine tegmental nucleus in the rat produces burst firing in A9 dopaminergic neurons. Neuroscience. 1999;92:245–254. [PubMed]
64. Floresco S, Todd C, Grace A. Glutamatergic afferents from the hippocampus to the nucleus accumbens regulate activity of ventral tegmental area dopamine neurons. Journal of Neuroscience. 2001;21:4915–4922. [PubMed]
65. Steffensen S, Svingos A, Pickel V, Henriksen S. Electrophysiological characterization of GABAergic neurons in the ventral tegmental area. Journal of neuroscience. 1998;18:8003–8015. [PubMed]
66. Guckenheimer J, Holmes P. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer; 1990.
67. Kuznetsov Y. Elements of applied bifurcation theory, volume 112 of Applied Mathematical Sciences. Springer-Verlag; 1995. pp. 295–403.
68. Dhooge A, Govaerts W, Kuznetsov Y. MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Transactions on Mathematical Software (TOMS) 2003;29:141–164.
69. Lodge D, Grace A. The hippocampus modulates dopamine neuron responsivity by regulating the intensity of phasic neuron activation. Neuropsychopharmacology. 2005;31:1356–1361. [PubMed]
70. Lodge D, Grace A. Aberrant hippocampal activity underlies the dopamine dysregulation in an animal model of schizophrenia. Journal of Neuroscience. 2007;27:11424–11430. [PubMed]
71. Fergestad T, Olson L, Patel K, Miller R, Palladino M, et al. Neuropathology in Drosophila mutants with increased seizure susceptibility. Genetics. 2008;178:947–956. [PubMed]
72. Meller R, Schindler C, Chu X, Xiong Z, Cameron J, et al. Seizure-like activity leads to the release of BAD from 14-3-3 protein and cell death in hippocampal neurons in vitro. Cell Death & Differentiation. 2003;10:539–547. [PubMed]
73. Kauffman S. The origins of order: Self organization and selection in evolution. USA: Oxford University Press; 1993.
74. Kitzbichler M, Smith M, Christensen S, Bullmore E. Broadband criticality of human brain network synchronization. PLoS computational biology. 2009;5:e1000314. [PMC free article] [PubMed]
75. Plenz D, Chialvo D. Scaling properties of neuronal avalanches are consistent with critical dynamics. Arxiv preprint arXiv. 2009:09125369.
76. Grace A, Bunney B. The control of firing pattern in nigral dopamine neurons: burst firing. Journal of neuroscience. 1984;4:2877–2890. [PubMed]
77. Johnson S, Seutin V, North R. Burst firing in dopamine neurons induced by N-methyl-D-aspartate: role of electrogenic sodium pump. Science. 1992;258:665–667. [PubMed]
78. He J, Hu B. Differential distribution of burst and single-spike responses in auditory thalamus. Journal of neurophysiology. 2002;88:2152–2156. [PubMed]
79. van Welie I, Remme M, van Hooft J, Wadman W. Different levels of Ih determine distinct temporal integration in bursting and regular-spiking neurons in rat subiculum. The Journal of Physiology. 2006;576:203–214. [PubMed]
80. Wilson C, Weyrick A, Terman D, Hallworth N, Bevan M. A model of reverse spike frequency adaptation and repetitive firing of subthalamic nucleus neurons. Journal of neurophysiology. 2004;91:1963–1980. [PubMed]
81. Brozović M, Abbott L, Andersen R. Mechanism of gain modulation at single neuron and network levels. Journal of Computational Neuroscience. 2008;25:158–168. [PubMed]
82. Marreiros A, Daunizeau J, Kiebel S, Friston K. Population dynamics: variance and the sigmoid activation function. Neuroimage. 2008;42:147–157. [PubMed]
83. Av-Ron E, Parnas H, Segel L. A minimal biophysical model for an excitable and oscillatory neuron. Biological cybernetics. 1991;65:487–500. [PubMed]
84. Terada K, Tanaka H, Yoshizawa S. Two-parameter bifurcations in the Hodgkin-Huxley equations for muscle fibers. Electronics and Communications in Japan(Part III Fundamental Electronic Science) 2000;83:86–94.
85. Tremblay L, Naranjo C, Graham S, Herrmann N, Mayberg H, et al. Functional neuroanatomical substrates of altered reward processing in major depressive disorder revealed by a dopaminergic probe. Archives of General Psychiatry. 2005;62:1228–1236. [PubMed]
86. Nestler E, Carlezon W. The mesolimbic dopamine reward circuit in depression. Biological psychiatry. 2006;59:1151–1159. [PubMed]
87. Silkis I. The role of dopamine-dependent negative feedback in the hippocampus-basal ganglia-thalamus-hippocampus loop in the extinction of responses. Neuroscience and Behavioral Physiology. 2008;38:399–405. [PubMed]

Articles from PLoS ONE are provided here courtesy of Public Library of Science