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

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 shortlasting 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][16][17][18], central pattern generator [19][20][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][24][25][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 X ) 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 twodimensional 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 X (represented by the one compartment, single cell minimal Av-Ron bursting model [29], described in Appendix S1) and the dynamics of the population X 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.
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.

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][37][38][39][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][45][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 X whose activity as a whole can be captured by a representative firing rate F . The rate F can be thought of as a weighted sum of the firing rates from the cells in X , or Figure 1. Comparison of single cell transitions between firing and quiescence, and its analog in the population firing-rate based model (6)- (7). A. We used the single cell biophysical bursting model of AvRon et al. [29] to illustrate the time evolution of the membrane potential V (t) (left) as well as the coupled phase-plane dynamics of the potential V and the slower recovery variable W (right). The equations in the reference and the parameters used for the simulation are given in Appendix S1. B. The population model (6)-(7) illustrates the time evolution of the firing rate F (t) (left) as well as the trajectories in the (b, F ) plane (right). It is notable that the single cell model employs three equations to trigger the bursts (i.e., transitions between quiescence and oscillations of the membrane potential), while our model captures the phenomenon simply as a high-low oscillation of the population firing rate, by a system of two equations. The parameters used in conjunction with our model to match the firing rate in B with the firing rate in A are also given in Appendix S1. doi:10.1371/journal.pone.0012695.g001 alternatively, as the firing rate of a prototypical cell in X (see Figure 1 and Appendix S1). It is assumed that F 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 X are integrated by means of a response function S with k S representing the gain of the response function (sec). S 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 y(t) will trigger changes in F within some delay t F as follows: withF F (t)*F (t) and F max representing, a time-weighted average of the recent firing between t and tzt F and, respectively, the maximum firing rate of the cells in X . 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 y is F~F max : S(y), 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), t F 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: In agreement with the analogy between cellular excitability and population excitability discussed above, we introduce a dynamic variable slower than F representing negative feedback that tends to decrease the firing response to input. We achieve this by letting the population input be a weighted sum y~aF {b max bzP, where aF is the intrinsic excitation, P is the net extrinsic excitation term, and b max b is the intrinsic dampening, with b [ ½0, 1 acting as a ''sliding control'' that tends to decrease the firing rate for a given input. The dynamics of b can be written as: with t b representing the time constant of the negative feedback, and with steady state function The parameter k b in Eq. (5) determines the maximum rate of change of b at steady state with respect to the firing rate F . The firing rate at which the steady state of b is 1/2 is F b . In single cell models, b would be similar to the gating of potassium channels that allows repolarization of the membrane during an action potential.
The gate b 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 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 X to range between zero and F max~4 00 Hz. However, F max is just a theoretical absolute maximum, and X 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 F max =2~200 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 y S~8 0 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 12 msec [36]. We estimate the typical response-time of the population (time Table 1. Parameters for the model (6)-(7). spent by X on integrating incoming input before responding with an action potential) to be significantly lower, and in our analysis we considered values from t F from 1=800 sec~1:25 msec, to 1=200 sec~5 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 P, on the order of 100{200 Hz, which is the effective rate integrated by the MDN to create an EPSP.
Intrinsic amplification. a [ ½0, 1 can be considered to act as a ''coupling index'', that is the fraction of the maximal electrical coupling (i.e., a~1 if all cells within X 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 17% to 40% at various development ranges [49]. Electrical studies found pair indices from 20% to 96%, and various average coupling conductances during the development [49]. To incorporate all theoretical possibilities, in our analysis we considered all values of a [ ½0, 1; however, as will be seen, most plausible dynamics occur for a [ ½0:1, 0:4, 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 (t b *0:03 sec). At its maximum capacity, we considered the effect of DA receptor autoregulation to be equivalent to that of a direct inhibitory input of b max = 160 Hz. The highest receptor sensitivity k b to extracellular dopamine, corresponding to ''half-capacity,'' is obtained at the firing rate F~F b . Since one of the aims of our model is to study the effects of receptor sensitivity on burst firing, we allowed F b to be a free parameter, varying widely between F b~0 and F b~2 00 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 y S and P 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 dF dt~0 or db dt~0 ), then calculating the equilibria (b Ã ,F Ã ) 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 y and for any choice of parameters, the response S(y) [ (0,1). Since dF =dt~{F z(F max {F )S(y), this implies that dF =dtv0 for F §F max and that dF =dtw0 for F ƒ0. Hence, if 0vF vF max , then {F vdF =dtvF max {2F , so that dF =dtv0 when F wF max =2. So all trajectories end asymptotically in the open strip 0vFvF max =2~200. In addition to this, one can also note that we always have 0vb ? (F )v1, which in turn implies that db=dtv0 for b §1 and that db=dtw0 for bƒ0, constraining the long term dynamics to the open strip 0vbv1. In conclusion: the open region R~(0, 1)|(0, 200) traps the asymptotic phase plane dynamics of system, in the sense that all trajectories eventually end in R and that any attracting sets (hence any features relevant to the long term dynamics like equilibria and cycles) are also contained in R.
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 db=dt~0 nullcline has sigmoidal shape, so is globally increasing, the piecewise monotonicity of the nullcline dF =dt~0 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 R implies, of course, that any phase plane cycle would have to also be contained in R, 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.
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 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 phaseplane 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.
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 limitpoint 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.

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 F (t) (Figure 2 left column), is paired in each case with the corresponding trajectory in the (b,F ) phaseplane ( 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 X . Figure 2A shows an example in which the system has a global stable equilibrium at F ? *33.9 Hz. When initialized at any other value, the firing may undergo an upturn, but always returns to the steady tonic rate F ? . For the parameter values in Figure 2A, repeated bursting is not possible. For the parameters in Figure 2B, any trajectory in the (b, F)-plane starts cycling. The firing rate describes an oscillation between *200 and *0 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 (b ? , F ? ) 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  , F ), describing a continuous curve as F b increases from 0 Hz to 200 Hz. The stability of the equilibria is color coded as follows: blue for stable spiral, red for unstable spiral. The two changes in stability are thus noticeable along the curve (red stars) and are caused by two subcritical Hopf bifurcations in the dynamics (paired with fold bifurcations, which are represented by the purple vertical lines). Cycles emerge at these bifurcations. The evolution of the large stable cycles as F b increases is represented by the blue dotted curve (described by its highest and lowest F values). Since the bistability windows between the fold and the Hopf bifurcations are very small, we inserted them as two magnified frames; in each of these frames, the fold bifurcation is visible as a purple vertical line, and tops and bottoms of the repelling cycle as F b increases are drawn as a red curve. In panel A, we fixed a = 0.5 and P = 120 Hz. B. The steady state curve shows two fold bifurcations (purple vertical lines), two Hopf bifurcations (red stars), and two limit point -or saddle node -bifurcations (light green stars). In this case a~0:75, P~100 Hz, and F b increasing from 0 Hz to 200 Hz. Additional explanations of the transitions in the dynamics and illustrations of the phase planes for these transitions are presented in Appendix S3. Common fixed parameters for A and B: k b = 0.025 sec, k S = 0.2 sec, y S = 80 Hz, b max = 160 Hz, t F = 2.5 msec, t b = 33 msec. doi:10.1371/journal.pone.0012695.g003 outside the region enclosed by the unstable limit cycle. The trajectory of the point inside the unstable limit cycle would go toward F ? *11 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 F ? *189 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, F b , are compared in Figure 3. The two panels in Figure 3 illustrate how the equilibria of the system change as F b increases between 0 to 200 Hz. In Figure 3A, the amplification parameter was set to a = 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 F b increases toward 30 Hz and disappears as F b approaches *140 Hz. The firing rates in the oscillatory regime alternate between *200 Hz (top dotted lines) and *0 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 X 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 F b 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 a = 0.75. The transition from lowtonic 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 (tonicbursting regime) We study first how our population responds to extrinsic excitation with and without the intrinsic amplification, respectively, for a~0 and aw0 ( Figure 4A-B). Each curve corresponds to a different value of F b . The steady state F ? of the population firing rate increases, as the frequency of extrinsic excitatory input, P, increases from 0 to 200 Hz.
As shown in Figure 4A, a unique, stable, and increasing steady state rate F ? (as a function of level of excitation P) 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 F ? increases as the dampening onset decreases (i.e., a larger F b rate places F ? 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 P for which F clearly departs from 0 depend on the halfactivation of dampening, F b . 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 aw0 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 6B), and extrinsic input ( Figure 6A). First, consider the changes in burst duration as a function of the intrinsic amplification weight a ( Figure 5). The larger a, 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 a (see also Figures 4B and 8B).
Other effects in the properties of sustained oscillations can be obtained by changing the external input P or the intrinsic dampening activation rate F b . Figure 6 shows how the activation of P and F b can be varied within relatively wide ranges (100{120 Hz for P and 80{140 Hz for F b , 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 t F and t b . 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 t F is that the right hand side of equation (6) will increase for smaller values of t F , yielding larger changes in F per unit time during an oscillation. On the other hand, increasing the time constant of the recovery variable, t b , decreases the timedependent change in b, which results in slower decrease during the oscillation in F . As a consequence, the frequency of oscillations decreases, without significantly altering the duty cycle (i.e. the burst/inter-burst duration ratio).

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 P, for different values of the half activation of intrinsic dampening, (F b = 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 P, for a relatively low onset of intrinsic dampening (F b~8 0 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 (F b~1 50 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 F b = 150 Hz, a saddle and a repeller appear for large enough P 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 P, 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 P. 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 There are no Hopf points on any of these curves, hence periodic bursting can't be triggered for a = 0. B. In the presence of sufficient intrinsic amplification, an extrinsic excitation appropriately balancing the dampening will trigger periodic bursting. Along these curves, Hopf points are marked with red stars, and limit points are marked with light green stars. The steady states are unstable (spirals or saddles) between the Hopf points, and stable otherwise (stability not shown). To maintain clarity of the diagram, and since their approximate placement is clear, the fold bifurcations and the corresponding cycles are not shown. Parameters: F b = 100 Hz, and increasing levels of intrinsic amplification. From lower to higher curve: a~0 (green curve, identical with the green curve in A), a = 0.1 (purple curve), a = 0.25 (cyan curve), a~0:5 (orange curve) and a~0:75 (blue curve). Fixed parameters common to A and B: k b = 0.025 sec, k S = 0.2 sec, y S = 80 Hz, P = 100 Hz, b max = 160 Hz, t F = 2.5 msec, t b = 33 msec. doi:10.1371/journal.pone.0012695.g004 Figure 5. Increasing the duration of the high-firing phase of the oscillation with minimal changes to the low-firing phase. The high firing rate portion of the oscillation (burst) can be increased by increasing the weight of the intrinsic amplification a. The figure shows the dependence on a of: the burst intervals (red curve), the inter-burst intervals (green curve) and their sum, i.e., the period of the oscillation (blue curve). As the weight of intrinsic amplification increases, the length of the bursts is longer, but the relaxation intervals between bursts remain approximately unchanged. Fixed parameter values: k S = 0.2, k b = 0.025, y S = 80 Hz, F b = 60 Hz, P = 120 Hz, b max = 160 Hz, t F = 2.5 msec, t b = 33 msec. doi:10.1371/journal.pone.0012695.g005 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 a, using different values of the dampening half activation (F b = 40, 80, and 120 Hz, respectively, in orange, green, and purple). Observe that the steady state firing rate F ? is noticeably higher for higher values of F b , and also, limit points are present for a larger range of F b .
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 F b 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 P in Figure 5A or increasing a in Figure 5B), or decreasing the intrinsic dampening (i.e., increasing F b 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 *200 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 F b ). 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 P may allow the return of F b to the preoscillatory range, without any need for additional protection against bursting. Conversely, raising P 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 F b 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 (P, a), (P, F b ) and (a, F b ). 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 (aw0), 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 F~200 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 20{30 Hz [42,76]; in other studies, the heterogeneity of types, localization and behavior of midbrain cells transpired as wider distributions, ranging from (20{30 Hz) to very high (up to 150{200 Hz) [39]. Similarly, some spike counts quote 2{5 spikes per burst [42], with relatively small interspike intervals of 6{12 msec [39]; however, more comprehensive descriptions have included higher quotes, and have additionally noticed that glutamate enhances bursting, increasing spike/burst counts to 8{10, while leaving inter-spike intervals relatively unchanged. Overall, accounts of burst duration have ranged widely from 20{120 msec [39,42], and inter-burst intervals from 50{250 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 pedunculopontine 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 2{10Hz, 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 X 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 k S and k b 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][86][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.