Network properties of a computational model of the ... - Semantic Scholar
Neural Networks (
)
–
Contents lists available at SciVerse ScienceDirect
Neural Networks journal homepage: www.elsevier.com/locate/neunet
2012 Special Issue
Network properties of a computational model of the dorsal raphe nucleus KongFatt Wong-Lin ∗ , Alok Joshi, Girijesh Prasad, T. Martin McGinnity Intelligent Systems Research Centre, University of Ulster, Magee Campus, Northland Road, BT48 7JL, Northern Ireland, UK
article
info
Keywords: Spiking neuronal network model Serotonin neurons Inhibitory fast-spiking non-serotonergic neurons Reward-based memory-guided decision task Theta rhythm
abstract Serotonin (5-HT) plays an important role in regulating mood, cognition and behaviour. The midbrain dorsal raphe nucleus (DRN) is one of the primary sources of 5-HT. Recent studies show that DRN neuronal activities can encode rewarding (e.g., appetitive) and unrewarding (e.g., aversive) behaviours. Experiments have also shown that DRN neurons can exhibit heterogeneous spiking behaviours. In this work, we build and study a basic spiking neuronal network model of the DRN constrained by neuronal properties observed in experiments. We use an efficient adaptive quadratic integrate-and-fire neuronal model to capture slow afterhyperpolarization current, occasional bursting behaviours in 5-HT neurons, and fast spiking activities in the non-5-HT inhibitory neurons. Provided that our noisy and heterogeneous spiking neuronal network model adopts a feedforward inhibitory network architecture, it is able to replicate the main features of DRN neuronal activities recorded in monkeys performing a reward-based memory-guided saccade task. The model exhibits theta band oscillation, especially among the non-5-HT inhibitory neurons during the rewarding outcome of a simulated trial, thus forming a model prediction. By varying the inhibitory synaptic strengths and the afferent inputs, we find that the network model can oscillate over a range of relatively low frequencies, allow co-existence of multiple stable frequencies, and spike synchrony can spread from within a local neural subgroup to global. Our model suggests plausible network architecture, provides interesting model predictions that can be experimentally tested, and offers a sufficiently realistic multi-scale model for 5-HT neuromodulation simulations. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction Serotonin (5-HT) is an important neurochemical that innervates throughout the brain, modulating neural activities, and thus regulating mood, cognition and behaviour (Giovanni, Matteo, & Esposito, 2008; Tseng & Atzori, 2007). The midbrain dorsal and median raphe nuclei are the sources of 5-HT. Within the dorsal raphe nucleus (DRN), about one- to two-thirds of the neurons contain 5-HT (Jacobs & Azmitia, 1992; Köhler & Steinbusch, 1982). Abnormal 5-HT activity levels have been implicated in devastating mental illnesses such as major depressive disorder, anxiety disorder and schizophrenia (Giovanni et al., 2008). Studies on psychotropic drugs, especially antidepressant drugs, have focused extensively on the 5-HT system, targeting its release, reuptake and various 5-HT receptors (Carr & Lucki, 2011; Cryan & Leonard, 2000; Giovanni et al., 2008). However, despite considerable research effort, an integrated understanding of how 5-HT and related drugs actually affect neuronal circuits in multifaceted ways remains unresolved (Pacher & Kecskemeti, 2004).
∗
Corresponding author. Tel.: +44 28 71375320; fax: +44 028 713 75435. E-mail addresses: [email protected] (K. Wong-Lin), [email protected] (A. Joshi), [email protected] (G. Prasad), [email protected] (T.M. McGinnity). 0893-6080/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.neunet.2012.02.009
Various electrophysiological, pharmacological, immunohistochemical and morphological studies of the raphe nuclei neurons have been conducted (Beck, Pan, Akanwa, & Kirby, 2004; Calizo et al., 2011; Giovanni et al., 2008; Hajós, Sharp, & Newberry, 1996; Kirby, Pernar, Valentino, & Beck, 2003; Kocsis, Varga, Dahan, & Sik, 2006; Li, Li, Kaneko, & Mizuno, 2001; Marinelli et al., 2004; Vandermaelen & Aghajanian, 1983). According to the classical and convenient identification, 5-HT neurons are presumed to have relatively slow and regular firing, broad action potentials, affinity to inhibitory feedback from 5-HT1A autoreceptors, or slow after-hyperpolarization (AHP) after a spike (Aghajanian & Vandermaelen, 1982; Hajós et al., 1996; Sprouse & Aghajanian, 1987; Vandermaelen & Aghajanian, 1983). However, more recent studies have shown that some of these suggested characteristics may not be highly reliable. Most electrophysiological properties are now found to be similar between 5-HT- and non-5-HT-containing neurons over various regions of the raphe nuclei (Allers & Sharp, 2003; Beck et al., 2004; Calizo et al., 2011; Kirby et al., 2003; Marinelli et al., 2004). Action potentials do not differ much between 5-HTand non-5-HT-containing neurons, and 5-HT1A autoreceptors can be found in both 5-HT- and non-5-HT-containing neurons (Calizo et al., 2011). Heterogeneity of neuronal activities in the DRN has also been observed in various experiments (Fig. 1), including behaving animals (Bromberg-Martin, Hikosaka, & Nakamura, 2010; Nakamura, Matsumoto, & Hikosaka, 2008; Ranade & Mainen,
2
K. Wong-Lin et al. / Neural Networks (
)
–
Fig. 1. Some heterogeneous neuronal properties found in experiments. (a–b) Decay time constants of afterhyperpolarization (AHP) in a 5-HT-containing (a) and a non-5HT-containing (b) neuron in the rat DRN (an in vitro intracellular study. (c–f) Variety of spiking behaviours of a presumably 5-HT neuron (c), a non-5-HT neuron (d), and a burst-firing 5-HT neuron (e), in the DRN of anaesthetized rat. Right panel of (e): asterisks denote spike doublets. Left panel: magnified to show spike doublets. (f) Spike singlets, doublets and triplets can be observed in a burst-firing 5-HT neuron. Scales for time and membrane potential magnitude are shown by the horizontal and vertical bars, respectively. Source: (a, b) adapted from Kirby et al. (2003) and (c–f) from Hajós et al. (1996).
2009). Recent advanced experimental tools have also shown that bursting DRN neurons, with similar electrophysiological properties as ‘‘classical’’ DRN neurons, are more reliably identified as 5-HT containing (Hajós et al., 2007, 1996; Kirby et al., 2003) (Fig. 1(e) and Fig. 1(f)). Overall, the mechanisms underlying these heterogeneous neuronal behaviours and their functions have yet to be fully understood. Previous computational models of the 5-HT systems have been useful in illuminating insights into either ‘‘low-level’’ single presynaptic terminal or ‘‘high-level’’ behavioural levels (Best, Nijhout, & Reed, 2010, 2011; Daw, Kakade, & Dayan, 2002; Dayan & Huys, 2009; Doya, 2002; Stoltenberg & Nag, 2010). However, neither accounts for the various neuronal properties of the 5-HT neurons nor attempts to understand at the intermediate circuit level. To further understand this complex and heterogeneous 5-HT system, it is pertinent to build multi-scale computational models with sufficient biological realism to provide more quantitative, systematic, predictive and integrated accounts of various neural and cognitive/behavioural effects (e.g. Wong-Lin, Prasad, & McGinnity, 2011). In this work, we make use of a spiking neuronal model to efficiently account for the heterogeneous DRN neuronal spiking behaviours. Then we combine the neurons to build a plausible network, comparing its activity to that recorded in animals performing a memory-guided decision task, and infer possible network architecture. Finally, we explore how the network model responds under various conditions, and discuss the implications. 2. Methods 2.1. Basic adaptive quadratic integrate-and-fire neuronal model To model the DRN neurons, we choose the quadratic integrateand-fire (QIF) neuronal model with a recovery variable as a compromise between computational efficiency and rich spiking dynamics properties (Izhikevich, 2003). The model, also known as the Izhikevich model, can be described by the coupled dynamical
equations of the membrane potential, V , and some recovery variable, U: dV dt dU dt
= 0.04V 2 + 5V + 140 − U + I
(1)
= a(bV − U )
(2)
with the after-spike resetting condition: if V ≤ Vpeak (in mV), then V is reset to a lower value c, and U is increased to U + d. Here Vpeak is the peak of an action potential, and a, b, c and d are constants. In our simulations, we found the parameter d does not affect our results much, and we fix it at 2. We retain the same value of parameter b at 0.2 as in Izhikevich (2003). We can emulate the slow AHP current by allowing slow decay dynamics of the recovery variable. This is achieved by setting c at higher values (∼ − 57 mV) and a = 0.005 such that recovery operates in the timescale of ∼200 ms (Hajós et al., 1996). As shown in Izhikevich (2003) and Wong-Lin et al. (2011), higher c values can result in bursting behaviours (two or more spikes per burst within less than 10 ms). This is consistent with experiments that had shown slow regular-spiking 5-HT neurons and bursting 5-HT neurons tend to have very similar basic electrophysiological properties. Inhibitory non-5-HT neurons are modelled in a way similar to the fast-spiking GABAergic interneurons in the cortex with parameters b = 0.28, and c = −65 mV (Izhikevich, 2003). The recovery variable U represents the activation of some inactivating currents (e.g. potassium ionic currents). For simplicity, we define I to be the afferent input current per membrane capacitance C (∼39 pF for both 5-HT neurons and non-5-HT neurons in Marinelli et al., 2004); hence the unit is in A/F, which we will henceforth not explicitly mention for the sake of brevity. The neuronal spiking behaviours are to some extent, similar to that of Hodgkin–Huxley type models as discussed in (Izhikevich, 2003, 2007). The term 0.04V2 + 5V + 140 is chosen so that V has a resting potential between −65 mV and −75 mV as observed in experiments (e.g. Hajós et al., 1996; Kirby et al., 2003), and more importantly to allow ease of extension to model other brain regions
K. Wong-Lin et al. / Neural Networks (
in future works (Izhikevich, 2003). We chose Vpeak = 0 mV to allow the spike amplitudes to be about 60–70 mV as in experiments (Hajós et al., 1996). There is experimental evidence that shows inhibitory synapses exert a stronger tonic influence than glutamatergic synapses in the DRN (Tao & Auerbach, 2003), which justify our focusing on the inhibitory neurons. Non-5-HT neurons can have almost all their electrophysiological characteristics similar to those of the 5-HT neurons, except a higher mean spontaneous firing rate, and a very fast AHP (Fig. 1(b)). Our model satisfies these non-5-HT neuronal criteria by adopting the neuronal model parameters a, b, c and d to be 0.02, 0.25, −67 mV and 2, respectively, similar to that of fastspiking inhibitory neurons in the cortex (Izhikevich, 2003). Thus, these neurons will have more sensitive response than the 5-HT neurons with respect to the input currents. We will not consider the smaller portion of non-5-HT neurons which contain glutamate, peptide transmitters, nitric oxide and other neurotransmitters (Allers & Sharp, 2003; Jacobs & Azmitia, 1992; Köhler & Steinbusch, 1982). 2.2. Synapses Inhibitory synapses: Inhibitory synaptic currents to all neurons are simulated with instantaneous dynamics (i.e. current-based) for simplicity: Ik = − gk
δ t − tj
(3)
j
where gk (>0) is the synaptic strength or weight on neuron k, and the sum is over Dirac-delta-like pre-synaptic spiking activities from all the inhibitory neurons. We can approximate such synapses when the latter has relatively fast dynamics (e.g. GABAA -mediated). Experimental recordings in rats’ DRN found that the decay time constant of inhibitory (GABAergic) synaptic gating variable is about 3.5 ms (4 ms at 32 °C in Lemos et al. (2006), and its rise time constant is very rapid at about 1 ms or less (Lemos et al., 2006). 2.3. Network We implement 300 neurons in our network, and adopt an all-toall connectivity architecture. To be consistent with experimental findings, we allow two-thirds (200 neurons) to be 5-HT-containing neurons while the rest (100 neurons) are non-5-HT-containing inhibitory neurons responsible for inhibitory synapses (e.g. Allers & Sharp, 2003). Although we have used a relatively small network (as compared to the actual number of ∼165 000 5-HT neurons in the human DRN (Baker et al., 1991)), the all-to-all connectivity allows scalability to a larger network (by appropriately adjusting the synaptic coupling strengths), and thus we expect our general results to hold for larger network size. We implement inhibitory synapses among the non-5-HT neurons, and from the non-5-HT neurons to 5-HT neurons. To prevent uncontrollable inhibitory firing, we have inhibitory synapses among the inhibitory neurons. Unless specified, the synaptic weight (gi ) from an inhibitory neuron to another inhibitory neuron, gi−i , is set to be at 2, while the synaptic weight from an inhibitory neuron to 5-HT neuron, gi−5HT , is 5. These values are used to allow the 5-HT (non-5-HT) neuronal population to attain a spontaneous or baseline firing rate of ∼0.5 Hz for 5-HT neurons and ∼5 Hz for the non-5-HT inhibitory neurons, consistent with experiments. Throughout any simulated trial, unless specified, there is a constant bias current to each neuron (3.57 for a 5-HT neuron, and 0.63 for a non-5-HT inhibitory neuron). Note that in this work, we
)
–
3
have adjusted the constant bias currents to be higher than that in Wong-Lin et al. (2011) to create more reasonable spontaneous firing rate levels (e.g. see Fig. 6). As compared to the single neuronal study, each neuron in the model receives a smaller amount of external Gaussian noise input in the afferent input I (standard deviation of 0.2 instead of 2) to compensate for the overall increase due to recurrent connections and heterogeneity within the network. Heterogeneity in neuronal excitability is generated before a simulated trial, by having the parameter c for both the 5-HT and non-5-HT neurons follow a Gaussian distribution with mean of −57 mV and −65 mV, respectively, and each a standard deviation of 2 mV. This ensures a subpopulation of cells to exhibit bursting behaviours while the rest produce regular single spikes. Furthermore, initial values of the membrane potential are randomized with Gaussian distribution with a mean of −65 mV and standard deviation of 5 mV. 2.4. Behavioural task in the experiments In the experiments of Bromberg-Martin et al. (2010) and Nakamura et al. (2008), subjects (monkeys) are trained to make fast eye movements (saccades) towards a visual target presented on a screen that may or may not be rewarded (see below). The subjects’ behavioural performance (reaction time and accuracy) and single neuronal spiking activity in the DRN are simultaneously recorded. Using the classical approach (based on spike duration, spiking regularity, and low baseline firing rate), putative 5-HT neurons in the DRN were identified and recorded during the tasks. Note that bursting neurons and fast-spiking non-5-HT neurons were not identified and recorded in these experiments. In the beginning of every experimental trial, the subjects were instructed to fixate their gaze on a visual fixation point on a screen (for 1000–1500 ms) prior to the appearance of a visual target either to the left or to the right of the fixation point (chosen pseudorandomly). In the memory-guided saccade (MGS) version of the task, a target flashed for a brief 100 ms duration before a delay period of 800 ms duration (with only the fixation point) appeared. The subject was required to maintain fixation of its gaze and remember the location of the target. At the end of this delay period the fixation point disappeared (Go cue) and the subject was required to make the appropriate saccade. A drop of liquid reward was delivered for the correct direction of the saccade made; no reward for an error. The reward was given after a 100 ms delay after the saccade and delivery duration was ∼500 ms. In a block of 20–32 trials, either the left or right choice target is rewarded. Without warning, the reward schedule will be switched to the other target location after a block of trials. This trial sequence is summarized in Fig. 2. The mean saccadic reaction times (defined as the time from Go cue to saccade made) were found to be ∼200 ms for the rewarding trials and slightly longer ∼280 ms for unrewarding trials (averaging over 2 subjects). 2.5. Behavioural task implemented in the model We assume for simplicity that one half of the 5-HT neurons (100 neurons) are used to encode reward and the other half (100 neurons) encode no reward. We implement the task stimuli and reward in our network model simply in the form of timevarying afferent inputs into our neuronal network. In addition to a constant bias current to all neurons, the total input current I, also includes a sum of input currents due to the presented stimuli: fixation point, visual targets, (un)rewarding outcome, and synaptic interactions among neurons. Depending on the trial type (rewarding or unrewarding), we insert an appropriate mean reaction time (mRT) when calculating the onset for rewarding or
4
K. Wong-Lin et al. / Neural Networks (
)
–
Fig. 2. Experimental set up. A typical trial in a memory-guided saccade task. Source: Adapted from Bromberg-Martin et al. (2010) with permission.
non-rewarding outcomes, in addition to a 100 ms delay in the reward. Each rewarding or unrewarding outcome lasts for 500 ms duration. Afferent excitatory inputs to the reward-encoding 5-HT neurons include that due to fixation point (0.1), target stimulus (0.1), and rewarding outcome (0.3); no input for non-rewarding outcome. Non-reward-encoding 5-HT neurons receive the same afferent visual stimulus inputs as reward-encoding 5-HT neurons except with a weaker input (0.1) due to the fixation point stimulus. However, non-reward encoding 5-HT neurons do not receive reward inputs, but instead receive inputs with unrewarding outcome (0.1). Afferent excitatory inputs to the non-5-HT inhibitory neurons turn out to be necessary to suppress neuronal activities at certain epochs within a trial. For simplicity, we assume non-5-HT inhibitory neurons receive the same input as that for reward-encoding neurons in rewarding trials, except that in nonrewarding outcomes they receive an additional input of 0.1, as for non-reward-encoding neurons. 2.6. Numerical simulations and analysis Numerical simulations: To mimic noise in the neuronal model, we include an additive noise Inoise to I with a random variable that follows a Gaussian distribution with a standard deviation of 2. This value was chosen to allow occasional bursts. Prior to every simulation, each neuron in the network has an initial random membrane potential such that its value is uniformly distributed from −65 mV to −15 mV. Time-averaged single-cell firing rates are averaged over 1.8 s. To numerically simulate the differential equations in the model, we use the forward Euler numerical scheme (Euler–Maruyama with noise) with a time step of 0.01 ms. Smaller time steps do not affect our results. Phase-plane analysis is carried out using the software XPPAUT. Population firing rate analysis: The average population firing rates are calculated by summing all the spikes in a single trial from a single sub-population or population of neurons divided by the total time taken in a trial. This is then repeated and averaged over several trials. Power spectral analysis: We emulate the local field potential by summing over all the membrane potentials of the sub-population of neurons. To calculate the power spectrum of oscillation for this summed membrane potentials, we use 105 sampling points per second; the total number of sampling points per second is calculated from dividing the total time per trial of 5500 ms by the time step of 0.01 ms and the total time of 5.5 s (5500/(0.01 ∗ 5.5) = 105 ). To reduce the size of the signal and to remove very high frequency components, we then downsampled the signal to 200 Hz (105 /500 = 200 Hz). According to Nyquist’s criterion, the maximum frequency is half of this value, i.e. 100 Hz. To select the frequency band, we fix a decision threshold at 55 dB.
3. Results 3.1. Single 5-HT neuronal model’s spiking patterns and mechanisms Before we study the neuronal network model, we shall first show that our neuronal model can demonstrate heterogeneity in spiking properties and is constrained by experiments. As mentioned earlier, most of the neuronal recordings in experiments performed on 5-HT containing neurons were sorted based on classically defined electrophysiological properties. In our model, we could conveniently capture such neuronal properties (e.g. Kirby et al., 2003; Li et al., 2001) by simply setting the reset membrane potential c at ∼− 60 mV or below (no bursting) while having some form of a spike frequency adaptation as shown in Fig. 3(a). The change in initial instantaneous firing rate was more gradual than its steady state over time (Fig. 3(b)) and applied current amplitude (Fig. 3(c)). In the absence of noise, the neuronal model exhibits a continuous change from no spiking to a finite spiking frequency of 2 Hz with a critical current (activation threshold) value where the change occurred very rapidly. Other than this brief initial rapid change, the slope of the model’s input–output function, or gain, is about 0.47 Hz/pA, which is comparable to the experimental values of neurons in several sites within the DRN (Crawford, Craige, & Beck, 2010). The addition of noise permits spiking at low firing rate and over a wider range of applied current (Fig. 3(c), dotted line). This is consistent with the observation that most putative 5-HT DRN neurons in cell cultures as well as in brain slice preparations do not generate spontaneous firing, unlike the in vivo intrinsically noisy condition where there is spontaneous synaptic inputs generating a range of spiking frequency observed in the range of ∼0.5–4 Hz (e.g. Allers & Sharp, 2003; Hajós et al., 1996; Sprouse & Aghajanian, 1987). As found in some experiments (e.g. Fig. 1(f)), our neuronal model can exhibit occasional bursting behaviour within a single noisy simulation even when its intrinsic neuronal properties in the noiseless condition does not permit bursting (spike doublets in Fig. 3(a), bottom). This is consistent with various experimental observations showing that a single neuron with slow regular spiking could exhibit occasional bursting of spike doublets and triplets especially in in vivo condition (e.g. Hajós et al., 2007, 1996; Rouchet et al., 2008; Vandermaelen & Aghajanian, 1983). Conceptually, this fluctuation-driven bursting behaviour can be understood by performing a phase-plane analysis similar to that in Wong-Lin et al. (2011) (see also Guckenheimer & Holmes, 1983; Izhikevich, 2007; Strogatz, 2001). In Fig. 3(d), we see that, given sufficient noise amplitude, an intrinsically non-bursting neuron can be perturbed out of the V-nullcline region in the presence of noise (dotted elliptical region), thus resulting in a quick second spike, i.e. creating a spike doublet (black trajectory). We would have expected this neuron to return to its refractory state after a single spike in the absence of this noise perturbation (grey dashed
K. Wong-Lin et al. / Neural Networks (
)
–
5
Fig. 3. Neuronal spiking properties and mechanisms. (a) Single neuronal model with c = −60 mV can show from spike singlets to doublets as the applied current (mean or noise) increases. Horizontal axis: time from current onset. A Gaussian noise with standard deviation 2.5 is used to generate noise-induced bursts e.g. spike doublets (labelled by asterisks). To be compared with (d). (b) Firing rate adaptation over time. Instantaneous firing rates calculated from the inverse of the interspike intervals over different time points. Initial instantaneous firing rates (calculate from the first 2 spikes) increase faster than the steady-state values as applied current increases. c = −60 mV. (c) Current-vs-frequency (input–output) curve under deterministic and noisy conditions. Noise can permit low (subcritical) firing rates. Dotted line: averaged over 30 trial simulations with Inoise = 0.2. (d) Fluctuation-driven bursts can be generated by noise even when the neuronal properties are not intrinsically bursting (with low reset c values). Phase plane of model superimposed with the same temporal evolution of membrane potential from bottom panel of panel (a). Grey (black) part of a single trajectory: periodic regular spike singlets (occasional noise-driven spike doublet). Although noise is more visible along the V-nullcline, it is added at every time step of the trajectory, which can occasionally drive the trajectory outside the V-nullcline resulting in a burst of spike doublet (elliptical region). Inset: zoomed in to the elliptical region; arrow denote point of escaping out of the V-nullcline.
trajectory). A similar argument for burst with more than two spikes can be made. 3.2. Neuronal network model during a reward-based memory-guided saccade task The above discussion on neuronal spiking mechanism would allow us to generate heterogeneous neuronal spiking patterns. We will use a network model to simulate DRN neuronal activity of non-human primates performing a simple perceptual decision task for both rewarding and unrewarding trials (Bromberg-Martin et al., 2010; Nakamura et al., 2008). This will then inform possible underlying neuronal circuit architecture, afferent inputs to the DRN circuit, and predict plausible responses from the inhibitory neurons. Before that, we shall first discuss the general results in the experiments of Bromberg-Martin et al. (2010) and Nakamura et al. (2008). 3.2.1. Main experimental results In Bromberg-Martin et al. (2010), three interesting types of DRN neurons were identified. One type of (putative 5-HT) DRN neurons elevated its activity during the reward epoch of the trial (Fig. 4, top, light grey), but not in an unrewarding trial (Fig. 4, top, dark grey). These neurons seemingly encode (from their elevated activities) rewards from target onset through the delay period, and especially towards the end of the trial from the Go cue through the rewarding epoch of the trial. Another group of (putative 5-HT) DRN neurons showed the reverse pattern: elevated activity during unrewarding trial, but suppressed activity during rewarding trial during the target-delay and outcome epochs. Fig. 4 shows that these two types of neurons showed firing rate activities in a reciprocal way.
A third group of (putative) non-5-HT neurons showed only phasic activation for any stimulus (fixation point or target) onset (not shown, but see Bromberg-Martin et al., 2010). Thus, their work shows that DRN neurons that are tonically excited during the task are also excited during the rewarding epoch of the trial, while DRN neurons that are tonically suppressed below baseline during the beginning of fixation onset, are excited during the unrewarding epoch of the trial. This latter group of negative or no reward encoding neurons could play an important role in encoding some form of ‘‘disappointment’’ to a lack of reward, which could be broadcasted to other parts of the brain. This could in turn modulate mood, cognition and subsequent behaviours. 3.2.2. Network model simulations of behavioural task We shall now attempt to mimic the stimulus inputs as in the experiments (compare Fig. 5(a) with Fig. 2). To simulate the neuronal activities found in the experiments with the simplest possible network architecture, we are constrained to implement an inhibitory feedforward network architecture (Fig. 5(b)). If we assume inhibition from inhibitory neurons to 5-HT neurons to be, on average, equal, we found that the fixation input current has to be much smaller for non-reward-encoding neurons than for rewardencoding neurons. The feedforward inhibition from the non-5-HT inhibitory neurons to the 5-HT neurons also has to be about 2.5 times that of inhibition among the inhibitory neurons. Satisfying these requirements allows the network model to replicate the essential neural activity patterns of the two types of DRN neurons found in the Bromberg-Martin et al. (2010) experiment described in 3.2.1 (compare Fig. 4 (top panel) with Fig. 5(c), and Fig. 4 (bottom panel) with Fig. 5(d)). Additionally, although our implemented input currents are step functions
6
K. Wong-Lin et al. / Neural Networks (
)
–
Fig. 4. Neuronal activities in the DRN recorded of behaving monkeys. Two types of neurons found in the study. Top (bottom): a sample (non-)reward-encoding neuronal spiking frequency averaged over trials. Light grey (dark grey): (un)rewarded trial. Double asterisks show significant changes in neuronal activities between the two trial types. Source: Adapted from Bromberg-Martin et al. (2010) with permission.
Fig. 5. Spiking neuronal network model. (a) Time course of afferent inputs to 5-HT and non-5-HT neurons in the model. Time course implemented to mimic the memoryguided saccade task. mRT: mean reaction time. (b) A feedforward network architecture. Inh: non-5-HT inhibitory neurons. + and −: excitatory and inhibitory synapses. 5-HT neurons separated into reward-encoding and non-reward-encoding types. (c, d) Single-trial spike raster diagram of 5-HT neurons (neuron #1–200) and non-5-HT inhibitory neurons (neuron #201–300) from our network model. Non-reward-encoding neurons (neuron #1–100); reward-encoding 5-HT neurons (neuron #101–200). Memory-guided saccade task was modelled for a rewarding (c) and an unrewarding (d) trial. Labels (i) to (iv) denote the epochs in a trial as in (a).
(Fig. 5(a)), the neurons are capable of firing phasically locked to stimulus onset followed by adaptation, thus resembling the spiking
activities in the experiments (Fig. 4). Although the target input is implemented briefly for only 100 ms, the non-reward-encoding
K. Wong-Lin et al. / Neural Networks (
Fig. 6. Spontaneous network activity with bursting behaviour. The averaged population firing rate for 5-HT (black, neuron #1–200) and non-5-HT inhibitory (grey, neuron #201–300) neurons are ∼0.5 Hz and ∼5 Hz, respectively. I5-HT = 3.57 and Ii = 0.63. Squares denote bursts (two or more spikes within ∼10 ms).
neuronal activities continue to be reduced as in the experiment. This is due to the lower fixation input current than that of rewardencoding neurons. During the Go cue (disappearance of fixation point and where subjects are required to make a saccadic choice), we find in the model that the rewarding signal has to begin from the Go signal onset till about 500 ms after the outcome onset. Otherwise, due to the absence of the fixation input, the neuronal spiking activities will drop to baseline (not observed in the recorded neurons) and not sustained over time (data not shown). This may suggest that the Go signal during a rewarding block (of trials) might actually be associated with the rewarding outcome on the reward-encoding neurons, perhaps through reinforcement (e.g. temporal difference) learning during the training phase of the subjects. Interestingly, the non-rewarding outcome input has to appear after a much longer delay than that of the rewarding outcome input. This means that non-rewarding outcome may not be directly attached to the immediate decision. As suggested in BrombergMartin et al. (2010), during this epoch, an afferent input may come from other brain regions, e.g. the lateral habenula, which seem to encode negative reward (Hikosaka, 2010). Lastly, we observed in our simulations that the non-5-HT inhibitory neurons can exhibit strong oscillations (∼12.5 Hz) especially during the rewarding outcome epoch (Fig. 5(c)). This distinct oscillatory behaviour may be due to the overall higher total input current into these neurons, thus creating more synchrony. 3.2.3. Spontaneous bursting in the network model When we remove all behavioural task stimuli, the 5-HT neuronal population in our network model can attain an averaged spontaneous firing rate of ∼0.5 Hz for the 5-HT neurons, comparable to the base values found in experiments (Hajós et al., 1996; Vandermaelen & Aghajanian, 1983). The non-5-HT inhibitory neuronal population typically can have a higher firing rate of ∼7 Hz, which is comparable to the ∼5 Hz in our model. In Fig. 6, we can clearly see this difference in spike frequencies between the 5-HT and non-5-HT neurons. We have allowed the model to exhibit occasional bursts over time by having noise (varying afferent input) and heterogeneity (varying neuronal parameter c) in the system. Because the bursting behaviour is an intrinsic property of the neurons, we can observe that a neuron that has produced a burst at an earlier time will have a higher chance of bursting again at a later time (Fig. 6, squares). 3.2.4. Effects of inhibitory synapses in the behavioural task Inhibitory synapses are known to have important effects on the DRN circuit (e.g. Tao & Auerbach, 2003). Thus they are also the target of some experimental drug studies. Studies have also shown the general importance of inhibitory neurons and synapses on oscillatory network behaviour (e.g. Whittington & Traub, 2003). Therefore, here, we shall investigate, in the context of our network
)
–
7
model, how the strengths or weights of the synapses can affect the network dynamics in this specific behavioural task paradigm. We will focus only on rewarding trials here as the non-rewarding trials show similar effects. We found that as the synaptic weights, gi−5-HT , increases from 5 to the more extreme value 50, 5-HT neurons are almost fully suppressed, except during the phasic bursts upon the onsets of stimuli, where the 5-HT neurons are synchronized within their sub-population. Furthermore, immediately after the rewarding epoch, the reward-encoding 5-HT neurons fire synchronously with the non-reward-encoding 5-HT neurons (Fig. 7(a)). This could be due to the sudden relaxation of inhibitory neuronal activity, having a huge simultaneous effect on all the 5-HT neurons. Thus, synchronous firing can spread from a local sub-population of 5-HT neurons to more global effects on all the 5-HT neurons within a single trial. When the synaptic weight, gi−i , is instead increased from 2 to 50, there is no dramatic change to the firing rates (compare Fig. 7(b) with Fig. 5(c)). However, at such high gi−i value there is some subtle increase in synchronous firing among the inhibitory neurons, resulting in slightly less suppression from the inhibitory neurons to the non-reward-encoding 5-HT neurons. 3.3. Synaptic strengths and afferent currents on population network firing rate and oscillation frequency 3.3.1. Population firing rates Having observed some possible network behavioural changes above, we shall now perform a more systematic study of the network effects due to changes in the afferent current amplitudes and synaptic weights. Unless specified to be varied, the parameter values gi−i , gi−5-HT , I5-HT , and Ii are fixed at 2, 5, 3.57 and 0.63, respectively, for studying changes in the population firing rates. We first vary the input current into the inhibitory neurons, Ii , and found that the gain or slope of the population current vs firing rate function systematically changes (Fig. 8(a)). We can understand this in a heuristic way. Suppose the inhibitory population firing rate at steady state can be described by ri = F (−gi−i ri + Ib ), where Ib is some activation threshold or bias current, and that F is approximately linear. Then effectively, ri = Ib ∗ gain where gain = 1/(1 + gi−i ). Thus increasing gi−i can decrease the overall gain of the inhibitory population. When gi−5-HT is varied, we do not observe similar change in the input–output gain of the 5-HT neuronal population (Fig. 8(b)), as changing the gi−5-HT is equivalent to changing the afferent inputs into the 5-HT neurons. Comparing the population input–output functions of the inhibitory neurons compared to the 5-HT neurons, we can see that the activation threshold (current) is lower while the gain is higher for the former than that of the latter. Thus, the inhibitory neurons in the model are more easily excited than the 5-HT neurons, consistent with experiments. Next, we gradually increase the synaptic weight gi−i . As expected, the overall inhibitory neuronal population firing rate generally drops as self-inhibition within the population is enhanced, and this is not affected by gi−5-HT since there is no feedback connection in this simplified model architecture (Fig. 8(c)). It is, however, interesting to note that the inhibitory population firing rate does not decrease too much as previously seen in Fig. 7(b). Increasing gi−5-HT produces a similar trend for the 5-HT neuronal population firing rate, but the change is now more dramatic (Fig. 8(d)). This change is attenuated with higher gi−i . This is due to the slightly lower firing rate of the inhibitory neurons (Fig. 8(c)), and thus there will be weaker inhibition from the inhibitory neurons to the 5-HT neurons. Note that the gradual saturation of the change in Fig. 8(d) as gi−i increases is directly related to the saturating change in Fig. 8(c).
8
K. Wong-Lin et al. / Neural Networks (
)
–
Fig. 7. Synchronous activities in the network model. Raster diagram with (a) gi−i = 2 and gi−5-HT = 50; (b) gi−i = 50 and gi−5-HT = 5. Label of time epoch as in Fig. 5(c).
Fig. 8. Change in 5-HT and non-5-HT neuronal population firing rates with various afferent currents (a–b) and synaptic weights (c–d).
3.3.2. Population oscillation frequency We shall now investigate how the 5-HT neuronal population oscillation frequency can change with different synaptic strengths. We first sum the membrane potentials of all the 200 5-HT neurons to mimic the local field potential of the 5-HT neuronal population. We then identify distinctive peak(s) in the power spectrum. To study network oscillations, we raise the afferent currents I5-HT , and Ii to 4.25 and 1.1, respectively, so that the 5-HT neuronal population can now oscillate about 4 Hz. We settle for a ∼4 Hz oscillation frequency (Fig. 9(a)) to constrain it to the typical range of ∼0.5–4 Hz (e.g. (Allers & Sharp, 2003; Hajós et al., 1996; Sprouse & Aghajanian, 1987). Also, there is some experimental evidence of theta band oscillation found in 5-HT neurons in the raphe nuclei (Kocsis et al., 2006; Kocsis & Vertes, 1992). The power spectrum in Fig. 9(b) clearly shows the distinct peak around this value, verifying our technique. When gi−5-HT synaptic strength is increased, multi-stable oscillation frequencies starts to appear (Fig. 9(c) and (d)). The
double peaks, one at about 3.9 Hz and the other at 12.5 Hz, appearing in the power spectrum in Fig. 9(d) confirm this. The 12.5 Hz higher frequency oscillation does not seem to be caused by single neuronal sub-threshold oscillation, but appear when we summed the membrane potentials together (Fig. 9(c) insets). As gi−5-HT increases to some critical value (∼20), the lower 3.9 Hz frequency decreases slightly to 3.1 Hz, while the higher frequency abruptly appears (Fig. 9(e)). If the synaptic strength gi−i is higher (20), the previous 12.5 Hz frequency is moved to an even higher value (∼37 Hz) without changing the lower frequency much. In Fig. 9(e), we can also see that the critical gi−5-HT is right-shifted. This phenomenon could be understood if we consider the gi−5-HT as a ‘‘gate’’ for oscillatory signals from the inhibitory population to influence the subthreshold membrane potentials of the 5-HT population. At certain critical gi−5-HT value, the gate is opened for information to flow from one population to another. This critical value has to be increased with higher gi−i as the lower inhibitory population activity is more suppressed, and a higher ‘‘gain’’ is needed to transmit the oscillatory information.
K. Wong-Lin et al. / Neural Networks (
)
–
9
Fig. 9. 5-HT population oscillation frequencies. (a,b) Control condition with about 4 Hz oscillation frequency in the summed membrane potentials (simulating local field potential). (c, d) Increased gi−5-HT synaptic strength results in multi-stable frequency. The higher frequencies do not come from single neuronal subthreshold oscillation (left inset in (c), single neuronal membrane potential timecourse) but appear when summed (right inset in (c), with 10 5-HT neurons). (d) Two peaks in the power spectrum (compared to one peak in (b)), one at 3.9 Hz and the other at 12.5 Hz. (b, d) A threshold of 55 dB is used to select the more distinctive oscillation frequencies. (e) Summary of the range of oscillation frequencies in the summed membrane potentials of 5-HT neurons with increasing gi−5-HT synaptic strength. With gi−i = 2 and gi−5-HT = 20, higher frequency appears at ∼12.5 Hz. With a higher gi−i value (e.g. of 20), higher oscillation frequencies can be obtained (e.g. ∼37 Hz, grey).
4. Discussion We have developed an efficient spiking neuronal network model of the DRN which consists of coupled 5-HT and non-5HT inhibitory neurons. At the neuronal level, we have shown, using computational simulations and dynamical systems analysis, that with variability in a slow recovery current, overall afferent input, and intrinsic neuronal noise, we can efficiently emulate heterogeneous spiking patterns (regular spiking, bursting or mixture of both). This is achieved even when the neurons have similar basic electrophysiological properties. These neuronal properties are actually consistent with experimental findings. With these constraints, we coupled the heterogeneous and noisy 5-HT-containing and non-5HT-containing inhibitory neurons to construct an important part of the DRN circuit network. The circuit can produce spontaneous bursting behaviour that is inherently a property of the neurons; neurons within a network that exhibit burst at some point will have a higher tendency to show bursting again at a later time. We constrain our network model architecture to a feedforward type so that the neuronal data from the Bromberg-Martin et al. (2010) experiment can be replicated. This suggests that inputs activating the inhibitory neurons may play an important role in providing fast inhibition to ‘‘competing’’ reward-encoding and non-reward-encoding 5-HT neurons. Input due to rewarding outcome (encoded by some reward-selective 5-HT neurons) seems to be activated during the Go cue rather than around the time of the actual reward, possibly suggesting some form of temporal difference learning is involved. Non-rewarding outcome (encoded by some non-reward-selective 5-HT neurons) appears after some significant delay, which suggests influence through different neural pathways, possibly via the lateral habenula. Increased synaptic strength of gi−5-HT not only suppresses 5-HT neuronal
spiking activity, but also enhances their tendency to synchronize, from local sub-population to globally. It can also decrease the gain of the inhibitory population input–output function. However, an increased self-inhibitory synaptic strength of gi−i does not seem to change the network as much, showing some form of robustness within the inhibitory population. We have seen how with the behavioural task stimuli that the network can show theta band oscillation especially around the rewarding outcome epoch of the trial. By increasing the constant bias input currents to the neurons, oscillations in the local field potentials of the network can appear consistently. We found that under certain parameter regimes, multiple oscillation frequencies can co-exist—multi-stability. Interestingly, under this regime of high gi−5-HT , the lower stable frequencies (∼2–4 Hz) seem to be related to the neuronal spike timings while the higher frequencies (∼12–40 Hz) of the 5-HT population’s local field potential is due to the transmission of the inhibitory population’s oscillatory signals. Theta oscillations of 5-HT and non-5-HT neurons in the raphe nuclei have been found in experiments (Kocsis et al., 2006; Kocsis & Vertes, 1992). This theta oscillation in the raphe nuclei can be phase locked to the theta oscillation of other brain regions, e.g. the hippocampal EEG (Kocsis et al., 2006), and a recent work has even suggested that serotonin can control hippocampal theta rhythm (Sörman, Wang, Hajós, & Kocsis, 2011). Thus theta oscillation may be a way to transfer information between the raphe nuclei and other brain regions like the hippocampus, and so it may be important to know whether coupling with other brain regions is necessary to generate theta oscillation in the raphe nuclei. Our model suggests that the DRN is capable of generating its own intrinsic theta oscillations. Slower oscillations e.g. ∼1–2 Hz have been observed in some experiments (Allers & Sharp, 2003; Kocsis et al., 2006; Wang & Aghajanian, 1982). Our simple network model is able to show
10
K. Wong-Lin et al. / Neural Networks (
oscillation frequency of about 2 Hz when the synaptic strength gi−5-HT is sufficiently high (Fig. 9(e)). Our model predicts that around this regime, we may also expect to observe high frequency oscillations in the local field potential of the 5-HT neurons. This could be verified in experiments. However, we have been unable to find parameter regimes where our network model can oscillate at a frequency lower than 2 Hz. Thus very slow oscillations may possibly be due to the presence of some slow ionic currents (e.g. mediated by NMDA, GABAB , or calcium) not implemented in our current model. Thus, a possible direction for future computational work may be to implement more specific ionic channels especially those with slower dynamics, e.g. (Kocsis et al., 2006; Penington, Kelly, & Fox, 1991; Rouchet et al., 2008). To conclude, we have successfully constructed an efficient neuronal circuit model of the DRN that may allow convenient investigation of serotonergic modulation in large-scale brain circuits. Moreover, the same adaptive spiking neuronal model can simulate the spiking properties of cortical and subcortical neurons with minimal parameter changes, even with incomplete experimental data (e.g. specific ionic channel properties). Incorporating genetic and intracellular dynamical processes (Best et al., 2010; Stoltenberg & Nag, 2010) into such models would result in a highly unified and multi-scale computational model. Acknowledgements The authors would like to thank Shahjahan Shahid for advice on power spectral analysis of noisy oscillatory signals. This study was initially supported under the CNRT award by the Northern Ireland Department for Employment and Learning through its ‘‘Strengthening the All-Island Research Base’’ initiative, and more recently via the Centre of Excellence in Intelligent Systems award, funded by InvestNI and the Integrated Development Fund, through its local facilitator, ILEX. References Aghajanian, G. K., & Vandermaelen, C. P. (1982). Intracellular identification of central noradrenergic and serotonergic neurons by a new double labeling procedure. The Journal of Neuroscience, 2, 1786–1792. Allers, K. A., & Sharp, T. (2003). Neurochemical and anatomical identification of fastand slow-firing neurones in the rat dorsal raphe nucleus using juxtacellular labelling methods in vivo. Neuroscience, 122, 193–204. Baker, K. G., Halliday, G. M., Hornung, J. P., Geffen, L. B., Cotton, R. G., & Törk, I. (1991). Distribution, morphology and number of monoaminesynthesizing and substance P-containing neurons in the human dorsal raphe nucleus. Neuroscience, 42(3), 757–775. Beck, S. G., Pan, Y.-Z., Akanwa, A. C., & Kirby, L. G. (2004). Median and dorsal raphe neurons are not electrophysiologically identical. Journal of Neurophysiology, 91, 994–1005. Best, J., Nijhout, H. F., & Reed, M. (2010). Serotonin synthesis, release and reuptake in terminals: a mathematical model. Theoretical Biology and Medical Modelling, 7, 34. doi:10.1186/1742-4682-7-34. Best, J., Nijhout, H. F., & Reed, M. (2011). Bursts and the efficacy of selective serotonin reuptake inhibitors. Pharmacopsychiatry, 44, S76–S83. doi:10.1055/s0031-1273697. Bromberg-Martin, E. S., Hikosaka, O., & Nakamura, K. (2010). Coding of task reward value in the dorsal raphe nucleus. The Journal of Neuroscience, 30, 6262–6272. Calizo, L. H., Ma, X., Pan, Y., Lemos, J., Craige, C., Heemstra, L. A., et al. (2011). Raphe serotonin neurons are not homogenous: electrophysiological, morphological and neurochemical evidence. Neuropharmacology, 61, 524–543. Carr, V. G., & Lucki, I. (2011). The role of serotonin receptor subtypes in treating depression: a review of animal studies. Psychopharmacology (Berl), 213, 265–287. Crawford, L. K., Craige, C. P., & Beck, S. G. (2010). Increased intrinsic excitability of lateral wing serotonin neurons of the dorsal raphe: a mechanism for selective activation in stress circuits. Journal of Neurophysiology, 103, 2652–2663. Cryan, J. F., & Leonard, B. E. (2000). 5-HT1A and beyond: the role of serotonin and its receptors in depression and the antidepressant response. Human Psychopharmacology, 15, 113–135. Daw, N. D., Kakade, D., & Dayan, P. (2002). Opponent interactions between serotonin and dopamine. Neural Networks, 15, 603–616.
)
–
Dayan, P., & Huys, Q. J. (2009). Serotonin in affective control. Annual Review of Neuroscience, 32, 95–126. Doya, K. (2002). Metalearning and neuromodulation. Neural Networks, 15, 495–506. Giovanni, G. D., Matteo, V. D., & Esposito, E. (2008). Progr. brain res.: Vol. 172. Serotonin-dopamine interaction: experimental evidence and therapeutic relevance. Elsevier Science. Guckenheimer, J., & Holmes, P. (1983). Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. New York: Springer-Verlag. Hajós, M., Allers, K. A., Jennings, K., Sharp, T., Charette, G., Sík, A., et al. (2007). Neurochemical identification of stereotypic burst-firing neurons in the rat dorsal raphe nucleus using juxtacellular labelling methods. European Journal of Neuroscience, 25, 119–126. Hajós, M., Sharp, T., & Newberry, N. R. (1996). Intracellular recordings from burstfiring presumed serotonergic neurones in the rat dorsal raphe nucleus in vivo. Brain Research, 737, 308–312. Hikosaka, O. (2010). The habenula: from stress evasion to value-based decisionmaking. Nature Reviews Neuroscience, 11, 503–513. Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Transactions on Neural Networks, 14, 1569–1572. Izhikevich, E. M. (2007). Dynamical systems in neuroscience. Cambridge, MA: MIT Press. Jacobs, B. L., & Azmitia, E. C. (1992). Structure and function of the brain serotonin system. Physiological Reviews, 72, 165–229. Kirby, L. G., Pernar, L., Valentino, R. J., & Beck, S. G. (2003). Distinguishing characteristics of serotonin and non-serotonin-containing cells in the dorsal raphe nucleus: electrophysiological and immunohistochemical studies. Neuroscience, 116, 669–683. Kocsis, B., Varga, V., Dahan, L., & Sik, A. (2006). Serotonergic neuron diversity: identification of raphe neurons with discharges time-locked to the hippocampal theta rhythm. Proceedings of the National Academy of Sciences USA, 103, 1059–1064. Kocsis, B., & Vertes, R. P. (1992). Dorsal raphe neurons: synchronous discharge with the theta rhythm of the hippocampus in the freely behaving rat. Journal of Neurophysiology, 68, 1463–1467. Köhler, C., & Steinbusch, H. (1982). Identification of serotonin and non-serotonincontaining neurons of the mid-brain raphe projecting to the entorhinal area and the hippocampal formation. A combined immunohistochemical and fluorescent retrograde tracing study in the rat brain. Neuroscience, 7, 951–975. Lemos, J. C., Pan, Y.-Z., Ma, X., Lamy, C., Akanwa, A. C., & Beck, S. G. (2006). Selective 5-HT1B receptor inhibition of glutamatergic and GABAergic synaptic activity in the rat dorsal and median raphe. European Journal of Neuroscience, 24, 3415–3430. Li, Y.-Q., Li, H., Kaneko, T., & Mizuno, N. (2001). Morphological features and electrophysiological properties of serotonergic and non-serotonergic projection neurons in the dorsal raphe nucleus. An intracellular recording and labeling study in rat brain slices. Brain Research, 900, 110–118. Marinelli, S., Schnell, S. A., Hack, S. P., Christie, M. J., Wessendorf, M. W., & Vaughan, C. W. (2004). Serotonergic and nonserotonergic dorsal raphe neurons are pharmacologically and electrophysiologically heterogeneous. Journal of Neurophysiology, 92, 3532–3537. Nakamura, K., Matsumoto, M., & Hikosaka, O. (2008). Reward-dependent modulation of neuronal activity in the primate dorsal raphe nucleus. The Journal of Neuroscience, 28, 5331–5343. Pacher, P., & Kecskemeti, V. (2004). Trends in the development of new antidepressants. Is there a light at the end of the tunnel? Current Medicinal Chemistry, 11, 925–943. Penington, N. J., Kelly, J. S., & Fox, A. P. (1991). A study of the mechanism of Ca2+ current inhibition produced by serotonin in rat dorsal raphe neurons. The Journal of Neuroscience, 11, 3594–3609. Ranade, S. P., & Mainen, Z. F. (2009). Transient firing of dorsal raphe neurons encodes diverse and specific sensory, motor, and reward events. Journal of Neurophysiology, 102, 3026–3037. Rouchet, N., Waroux, O., Lamy, C., Massotte, L., Scuvée-Moreau, J., Liégeois, J. F., et al. (2008). SK channel blockade promotes burst firing in dorsal raphe serotonergic neurons. European Journal of Neuroscience, 28, 1108–1115. Sörman, E., Wang, D., Hajós, M., & Kocsis, B. (2011). Control of hippocampal theta rhythm by serotonin: role of 5-HT2c receptors. Neuropharmacology, 61, 489–494. Sprouse, J. S., & Aghajanian, G. K. (1987). Electrophysiological responses of serotoninergic dorsal raphe neurons to 5-HT1A and 5-HT1B agonists. Synapse, 1, 582–588. Stoltenberg, S. F., & Nag, P. (2010). Description and validation of a dynamical systems model of presynaptic serotonin function: genetic variation, brain activation and impulsivity. Behavior Genetics, 40, 262–279. Strogatz, S. H. (2001). Nonlinear dynamics and chaos: with applications to physics, biology, chemistry and engineering. Cambridge, MA: Perseus Books Group. Tao, R., & Auerbach, S. B. (2003). Influence of inhibitory and excitatory inputs on serotonin efflux differs in the dorsal raphe and median raphe nuclei. Brain Research, 961, 109–120. Tseng, K.-Y., & Atzori, M. (2007). Monoaminergic modulation of cortical excitability. Springer.
K. Wong-Lin et al. / Neural Networks ( Vandermaelen, C. P., & Aghajanian, G. K. (1983). Electrophysiological and pharmacological characterization of serotonergic dorsal raphe neurons recorded extracellularly and intracellularly in rat brain slices. Brain Research, 289, 109–119. Wang, R. Y., & Aghajanian, G. K. (1982). Correlative firing patterns of serotonergic neurons in rat dorsal raphe nucleus. The Journal of Neuroscience, 2, 11–16.
)
–
11
Whittington, M. A., & Traub, R. D. (2003). Interneuron diversity series: inhibitory interneurons and network oscillations in vitro. Trends in Neurosciences, 26, 676–682. Wong-Lin, K. F., Prasad, G., & McGinnity, T. M. (2011). A spiking neuronal network model of the dorsal raphe nucleus. In Proceedings of 2011 international joint conference on neural networks. San Jose, CA, USA.
)
–
Contents lists available at SciVerse ScienceDirect
Neural Networks journal homepage: www.elsevier.com/locate/neunet
2012 Special Issue
Network properties of a computational model of the dorsal raphe nucleus KongFatt Wong-Lin ∗ , Alok Joshi, Girijesh Prasad, T. Martin McGinnity Intelligent Systems Research Centre, University of Ulster, Magee Campus, Northland Road, BT48 7JL, Northern Ireland, UK
article
info
Keywords: Spiking neuronal network model Serotonin neurons Inhibitory fast-spiking non-serotonergic neurons Reward-based memory-guided decision task Theta rhythm
abstract Serotonin (5-HT) plays an important role in regulating mood, cognition and behaviour. The midbrain dorsal raphe nucleus (DRN) is one of the primary sources of 5-HT. Recent studies show that DRN neuronal activities can encode rewarding (e.g., appetitive) and unrewarding (e.g., aversive) behaviours. Experiments have also shown that DRN neurons can exhibit heterogeneous spiking behaviours. In this work, we build and study a basic spiking neuronal network model of the DRN constrained by neuronal properties observed in experiments. We use an efficient adaptive quadratic integrate-and-fire neuronal model to capture slow afterhyperpolarization current, occasional bursting behaviours in 5-HT neurons, and fast spiking activities in the non-5-HT inhibitory neurons. Provided that our noisy and heterogeneous spiking neuronal network model adopts a feedforward inhibitory network architecture, it is able to replicate the main features of DRN neuronal activities recorded in monkeys performing a reward-based memory-guided saccade task. The model exhibits theta band oscillation, especially among the non-5-HT inhibitory neurons during the rewarding outcome of a simulated trial, thus forming a model prediction. By varying the inhibitory synaptic strengths and the afferent inputs, we find that the network model can oscillate over a range of relatively low frequencies, allow co-existence of multiple stable frequencies, and spike synchrony can spread from within a local neural subgroup to global. Our model suggests plausible network architecture, provides interesting model predictions that can be experimentally tested, and offers a sufficiently realistic multi-scale model for 5-HT neuromodulation simulations. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction Serotonin (5-HT) is an important neurochemical that innervates throughout the brain, modulating neural activities, and thus regulating mood, cognition and behaviour (Giovanni, Matteo, & Esposito, 2008; Tseng & Atzori, 2007). The midbrain dorsal and median raphe nuclei are the sources of 5-HT. Within the dorsal raphe nucleus (DRN), about one- to two-thirds of the neurons contain 5-HT (Jacobs & Azmitia, 1992; Köhler & Steinbusch, 1982). Abnormal 5-HT activity levels have been implicated in devastating mental illnesses such as major depressive disorder, anxiety disorder and schizophrenia (Giovanni et al., 2008). Studies on psychotropic drugs, especially antidepressant drugs, have focused extensively on the 5-HT system, targeting its release, reuptake and various 5-HT receptors (Carr & Lucki, 2011; Cryan & Leonard, 2000; Giovanni et al., 2008). However, despite considerable research effort, an integrated understanding of how 5-HT and related drugs actually affect neuronal circuits in multifaceted ways remains unresolved (Pacher & Kecskemeti, 2004).
∗
Corresponding author. Tel.: +44 28 71375320; fax: +44 028 713 75435. E-mail addresses: [email protected] (K. Wong-Lin), [email protected] (A. Joshi), [email protected] (G. Prasad), [email protected] (T.M. McGinnity). 0893-6080/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.neunet.2012.02.009
Various electrophysiological, pharmacological, immunohistochemical and morphological studies of the raphe nuclei neurons have been conducted (Beck, Pan, Akanwa, & Kirby, 2004; Calizo et al., 2011; Giovanni et al., 2008; Hajós, Sharp, & Newberry, 1996; Kirby, Pernar, Valentino, & Beck, 2003; Kocsis, Varga, Dahan, & Sik, 2006; Li, Li, Kaneko, & Mizuno, 2001; Marinelli et al., 2004; Vandermaelen & Aghajanian, 1983). According to the classical and convenient identification, 5-HT neurons are presumed to have relatively slow and regular firing, broad action potentials, affinity to inhibitory feedback from 5-HT1A autoreceptors, or slow after-hyperpolarization (AHP) after a spike (Aghajanian & Vandermaelen, 1982; Hajós et al., 1996; Sprouse & Aghajanian, 1987; Vandermaelen & Aghajanian, 1983). However, more recent studies have shown that some of these suggested characteristics may not be highly reliable. Most electrophysiological properties are now found to be similar between 5-HT- and non-5-HT-containing neurons over various regions of the raphe nuclei (Allers & Sharp, 2003; Beck et al., 2004; Calizo et al., 2011; Kirby et al., 2003; Marinelli et al., 2004). Action potentials do not differ much between 5-HTand non-5-HT-containing neurons, and 5-HT1A autoreceptors can be found in both 5-HT- and non-5-HT-containing neurons (Calizo et al., 2011). Heterogeneity of neuronal activities in the DRN has also been observed in various experiments (Fig. 1), including behaving animals (Bromberg-Martin, Hikosaka, & Nakamura, 2010; Nakamura, Matsumoto, & Hikosaka, 2008; Ranade & Mainen,
2
K. Wong-Lin et al. / Neural Networks (
)
–
Fig. 1. Some heterogeneous neuronal properties found in experiments. (a–b) Decay time constants of afterhyperpolarization (AHP) in a 5-HT-containing (a) and a non-5HT-containing (b) neuron in the rat DRN (an in vitro intracellular study. (c–f) Variety of spiking behaviours of a presumably 5-HT neuron (c), a non-5-HT neuron (d), and a burst-firing 5-HT neuron (e), in the DRN of anaesthetized rat. Right panel of (e): asterisks denote spike doublets. Left panel: magnified to show spike doublets. (f) Spike singlets, doublets and triplets can be observed in a burst-firing 5-HT neuron. Scales for time and membrane potential magnitude are shown by the horizontal and vertical bars, respectively. Source: (a, b) adapted from Kirby et al. (2003) and (c–f) from Hajós et al. (1996).
2009). Recent advanced experimental tools have also shown that bursting DRN neurons, with similar electrophysiological properties as ‘‘classical’’ DRN neurons, are more reliably identified as 5-HT containing (Hajós et al., 2007, 1996; Kirby et al., 2003) (Fig. 1(e) and Fig. 1(f)). Overall, the mechanisms underlying these heterogeneous neuronal behaviours and their functions have yet to be fully understood. Previous computational models of the 5-HT systems have been useful in illuminating insights into either ‘‘low-level’’ single presynaptic terminal or ‘‘high-level’’ behavioural levels (Best, Nijhout, & Reed, 2010, 2011; Daw, Kakade, & Dayan, 2002; Dayan & Huys, 2009; Doya, 2002; Stoltenberg & Nag, 2010). However, neither accounts for the various neuronal properties of the 5-HT neurons nor attempts to understand at the intermediate circuit level. To further understand this complex and heterogeneous 5-HT system, it is pertinent to build multi-scale computational models with sufficient biological realism to provide more quantitative, systematic, predictive and integrated accounts of various neural and cognitive/behavioural effects (e.g. Wong-Lin, Prasad, & McGinnity, 2011). In this work, we make use of a spiking neuronal model to efficiently account for the heterogeneous DRN neuronal spiking behaviours. Then we combine the neurons to build a plausible network, comparing its activity to that recorded in animals performing a memory-guided decision task, and infer possible network architecture. Finally, we explore how the network model responds under various conditions, and discuss the implications. 2. Methods 2.1. Basic adaptive quadratic integrate-and-fire neuronal model To model the DRN neurons, we choose the quadratic integrateand-fire (QIF) neuronal model with a recovery variable as a compromise between computational efficiency and rich spiking dynamics properties (Izhikevich, 2003). The model, also known as the Izhikevich model, can be described by the coupled dynamical
equations of the membrane potential, V , and some recovery variable, U: dV dt dU dt
= 0.04V 2 + 5V + 140 − U + I
(1)
= a(bV − U )
(2)
with the after-spike resetting condition: if V ≤ Vpeak (in mV), then V is reset to a lower value c, and U is increased to U + d. Here Vpeak is the peak of an action potential, and a, b, c and d are constants. In our simulations, we found the parameter d does not affect our results much, and we fix it at 2. We retain the same value of parameter b at 0.2 as in Izhikevich (2003). We can emulate the slow AHP current by allowing slow decay dynamics of the recovery variable. This is achieved by setting c at higher values (∼ − 57 mV) and a = 0.005 such that recovery operates in the timescale of ∼200 ms (Hajós et al., 1996). As shown in Izhikevich (2003) and Wong-Lin et al. (2011), higher c values can result in bursting behaviours (two or more spikes per burst within less than 10 ms). This is consistent with experiments that had shown slow regular-spiking 5-HT neurons and bursting 5-HT neurons tend to have very similar basic electrophysiological properties. Inhibitory non-5-HT neurons are modelled in a way similar to the fast-spiking GABAergic interneurons in the cortex with parameters b = 0.28, and c = −65 mV (Izhikevich, 2003). The recovery variable U represents the activation of some inactivating currents (e.g. potassium ionic currents). For simplicity, we define I to be the afferent input current per membrane capacitance C (∼39 pF for both 5-HT neurons and non-5-HT neurons in Marinelli et al., 2004); hence the unit is in A/F, which we will henceforth not explicitly mention for the sake of brevity. The neuronal spiking behaviours are to some extent, similar to that of Hodgkin–Huxley type models as discussed in (Izhikevich, 2003, 2007). The term 0.04V2 + 5V + 140 is chosen so that V has a resting potential between −65 mV and −75 mV as observed in experiments (e.g. Hajós et al., 1996; Kirby et al., 2003), and more importantly to allow ease of extension to model other brain regions
K. Wong-Lin et al. / Neural Networks (
in future works (Izhikevich, 2003). We chose Vpeak = 0 mV to allow the spike amplitudes to be about 60–70 mV as in experiments (Hajós et al., 1996). There is experimental evidence that shows inhibitory synapses exert a stronger tonic influence than glutamatergic synapses in the DRN (Tao & Auerbach, 2003), which justify our focusing on the inhibitory neurons. Non-5-HT neurons can have almost all their electrophysiological characteristics similar to those of the 5-HT neurons, except a higher mean spontaneous firing rate, and a very fast AHP (Fig. 1(b)). Our model satisfies these non-5-HT neuronal criteria by adopting the neuronal model parameters a, b, c and d to be 0.02, 0.25, −67 mV and 2, respectively, similar to that of fastspiking inhibitory neurons in the cortex (Izhikevich, 2003). Thus, these neurons will have more sensitive response than the 5-HT neurons with respect to the input currents. We will not consider the smaller portion of non-5-HT neurons which contain glutamate, peptide transmitters, nitric oxide and other neurotransmitters (Allers & Sharp, 2003; Jacobs & Azmitia, 1992; Köhler & Steinbusch, 1982). 2.2. Synapses Inhibitory synapses: Inhibitory synaptic currents to all neurons are simulated with instantaneous dynamics (i.e. current-based) for simplicity: Ik = − gk
δ t − tj
(3)
j
where gk (>0) is the synaptic strength or weight on neuron k, and the sum is over Dirac-delta-like pre-synaptic spiking activities from all the inhibitory neurons. We can approximate such synapses when the latter has relatively fast dynamics (e.g. GABAA -mediated). Experimental recordings in rats’ DRN found that the decay time constant of inhibitory (GABAergic) synaptic gating variable is about 3.5 ms (4 ms at 32 °C in Lemos et al. (2006), and its rise time constant is very rapid at about 1 ms or less (Lemos et al., 2006). 2.3. Network We implement 300 neurons in our network, and adopt an all-toall connectivity architecture. To be consistent with experimental findings, we allow two-thirds (200 neurons) to be 5-HT-containing neurons while the rest (100 neurons) are non-5-HT-containing inhibitory neurons responsible for inhibitory synapses (e.g. Allers & Sharp, 2003). Although we have used a relatively small network (as compared to the actual number of ∼165 000 5-HT neurons in the human DRN (Baker et al., 1991)), the all-to-all connectivity allows scalability to a larger network (by appropriately adjusting the synaptic coupling strengths), and thus we expect our general results to hold for larger network size. We implement inhibitory synapses among the non-5-HT neurons, and from the non-5-HT neurons to 5-HT neurons. To prevent uncontrollable inhibitory firing, we have inhibitory synapses among the inhibitory neurons. Unless specified, the synaptic weight (gi ) from an inhibitory neuron to another inhibitory neuron, gi−i , is set to be at 2, while the synaptic weight from an inhibitory neuron to 5-HT neuron, gi−5HT , is 5. These values are used to allow the 5-HT (non-5-HT) neuronal population to attain a spontaneous or baseline firing rate of ∼0.5 Hz for 5-HT neurons and ∼5 Hz for the non-5-HT inhibitory neurons, consistent with experiments. Throughout any simulated trial, unless specified, there is a constant bias current to each neuron (3.57 for a 5-HT neuron, and 0.63 for a non-5-HT inhibitory neuron). Note that in this work, we
)
–
3
have adjusted the constant bias currents to be higher than that in Wong-Lin et al. (2011) to create more reasonable spontaneous firing rate levels (e.g. see Fig. 6). As compared to the single neuronal study, each neuron in the model receives a smaller amount of external Gaussian noise input in the afferent input I (standard deviation of 0.2 instead of 2) to compensate for the overall increase due to recurrent connections and heterogeneity within the network. Heterogeneity in neuronal excitability is generated before a simulated trial, by having the parameter c for both the 5-HT and non-5-HT neurons follow a Gaussian distribution with mean of −57 mV and −65 mV, respectively, and each a standard deviation of 2 mV. This ensures a subpopulation of cells to exhibit bursting behaviours while the rest produce regular single spikes. Furthermore, initial values of the membrane potential are randomized with Gaussian distribution with a mean of −65 mV and standard deviation of 5 mV. 2.4. Behavioural task in the experiments In the experiments of Bromberg-Martin et al. (2010) and Nakamura et al. (2008), subjects (monkeys) are trained to make fast eye movements (saccades) towards a visual target presented on a screen that may or may not be rewarded (see below). The subjects’ behavioural performance (reaction time and accuracy) and single neuronal spiking activity in the DRN are simultaneously recorded. Using the classical approach (based on spike duration, spiking regularity, and low baseline firing rate), putative 5-HT neurons in the DRN were identified and recorded during the tasks. Note that bursting neurons and fast-spiking non-5-HT neurons were not identified and recorded in these experiments. In the beginning of every experimental trial, the subjects were instructed to fixate their gaze on a visual fixation point on a screen (for 1000–1500 ms) prior to the appearance of a visual target either to the left or to the right of the fixation point (chosen pseudorandomly). In the memory-guided saccade (MGS) version of the task, a target flashed for a brief 100 ms duration before a delay period of 800 ms duration (with only the fixation point) appeared. The subject was required to maintain fixation of its gaze and remember the location of the target. At the end of this delay period the fixation point disappeared (Go cue) and the subject was required to make the appropriate saccade. A drop of liquid reward was delivered for the correct direction of the saccade made; no reward for an error. The reward was given after a 100 ms delay after the saccade and delivery duration was ∼500 ms. In a block of 20–32 trials, either the left or right choice target is rewarded. Without warning, the reward schedule will be switched to the other target location after a block of trials. This trial sequence is summarized in Fig. 2. The mean saccadic reaction times (defined as the time from Go cue to saccade made) were found to be ∼200 ms for the rewarding trials and slightly longer ∼280 ms for unrewarding trials (averaging over 2 subjects). 2.5. Behavioural task implemented in the model We assume for simplicity that one half of the 5-HT neurons (100 neurons) are used to encode reward and the other half (100 neurons) encode no reward. We implement the task stimuli and reward in our network model simply in the form of timevarying afferent inputs into our neuronal network. In addition to a constant bias current to all neurons, the total input current I, also includes a sum of input currents due to the presented stimuli: fixation point, visual targets, (un)rewarding outcome, and synaptic interactions among neurons. Depending on the trial type (rewarding or unrewarding), we insert an appropriate mean reaction time (mRT) when calculating the onset for rewarding or
4
K. Wong-Lin et al. / Neural Networks (
)
–
Fig. 2. Experimental set up. A typical trial in a memory-guided saccade task. Source: Adapted from Bromberg-Martin et al. (2010) with permission.
non-rewarding outcomes, in addition to a 100 ms delay in the reward. Each rewarding or unrewarding outcome lasts for 500 ms duration. Afferent excitatory inputs to the reward-encoding 5-HT neurons include that due to fixation point (0.1), target stimulus (0.1), and rewarding outcome (0.3); no input for non-rewarding outcome. Non-reward-encoding 5-HT neurons receive the same afferent visual stimulus inputs as reward-encoding 5-HT neurons except with a weaker input (0.1) due to the fixation point stimulus. However, non-reward encoding 5-HT neurons do not receive reward inputs, but instead receive inputs with unrewarding outcome (0.1). Afferent excitatory inputs to the non-5-HT inhibitory neurons turn out to be necessary to suppress neuronal activities at certain epochs within a trial. For simplicity, we assume non-5-HT inhibitory neurons receive the same input as that for reward-encoding neurons in rewarding trials, except that in nonrewarding outcomes they receive an additional input of 0.1, as for non-reward-encoding neurons. 2.6. Numerical simulations and analysis Numerical simulations: To mimic noise in the neuronal model, we include an additive noise Inoise to I with a random variable that follows a Gaussian distribution with a standard deviation of 2. This value was chosen to allow occasional bursts. Prior to every simulation, each neuron in the network has an initial random membrane potential such that its value is uniformly distributed from −65 mV to −15 mV. Time-averaged single-cell firing rates are averaged over 1.8 s. To numerically simulate the differential equations in the model, we use the forward Euler numerical scheme (Euler–Maruyama with noise) with a time step of 0.01 ms. Smaller time steps do not affect our results. Phase-plane analysis is carried out using the software XPPAUT. Population firing rate analysis: The average population firing rates are calculated by summing all the spikes in a single trial from a single sub-population or population of neurons divided by the total time taken in a trial. This is then repeated and averaged over several trials. Power spectral analysis: We emulate the local field potential by summing over all the membrane potentials of the sub-population of neurons. To calculate the power spectrum of oscillation for this summed membrane potentials, we use 105 sampling points per second; the total number of sampling points per second is calculated from dividing the total time per trial of 5500 ms by the time step of 0.01 ms and the total time of 5.5 s (5500/(0.01 ∗ 5.5) = 105 ). To reduce the size of the signal and to remove very high frequency components, we then downsampled the signal to 200 Hz (105 /500 = 200 Hz). According to Nyquist’s criterion, the maximum frequency is half of this value, i.e. 100 Hz. To select the frequency band, we fix a decision threshold at 55 dB.
3. Results 3.1. Single 5-HT neuronal model’s spiking patterns and mechanisms Before we study the neuronal network model, we shall first show that our neuronal model can demonstrate heterogeneity in spiking properties and is constrained by experiments. As mentioned earlier, most of the neuronal recordings in experiments performed on 5-HT containing neurons were sorted based on classically defined electrophysiological properties. In our model, we could conveniently capture such neuronal properties (e.g. Kirby et al., 2003; Li et al., 2001) by simply setting the reset membrane potential c at ∼− 60 mV or below (no bursting) while having some form of a spike frequency adaptation as shown in Fig. 3(a). The change in initial instantaneous firing rate was more gradual than its steady state over time (Fig. 3(b)) and applied current amplitude (Fig. 3(c)). In the absence of noise, the neuronal model exhibits a continuous change from no spiking to a finite spiking frequency of 2 Hz with a critical current (activation threshold) value where the change occurred very rapidly. Other than this brief initial rapid change, the slope of the model’s input–output function, or gain, is about 0.47 Hz/pA, which is comparable to the experimental values of neurons in several sites within the DRN (Crawford, Craige, & Beck, 2010). The addition of noise permits spiking at low firing rate and over a wider range of applied current (Fig. 3(c), dotted line). This is consistent with the observation that most putative 5-HT DRN neurons in cell cultures as well as in brain slice preparations do not generate spontaneous firing, unlike the in vivo intrinsically noisy condition where there is spontaneous synaptic inputs generating a range of spiking frequency observed in the range of ∼0.5–4 Hz (e.g. Allers & Sharp, 2003; Hajós et al., 1996; Sprouse & Aghajanian, 1987). As found in some experiments (e.g. Fig. 1(f)), our neuronal model can exhibit occasional bursting behaviour within a single noisy simulation even when its intrinsic neuronal properties in the noiseless condition does not permit bursting (spike doublets in Fig. 3(a), bottom). This is consistent with various experimental observations showing that a single neuron with slow regular spiking could exhibit occasional bursting of spike doublets and triplets especially in in vivo condition (e.g. Hajós et al., 2007, 1996; Rouchet et al., 2008; Vandermaelen & Aghajanian, 1983). Conceptually, this fluctuation-driven bursting behaviour can be understood by performing a phase-plane analysis similar to that in Wong-Lin et al. (2011) (see also Guckenheimer & Holmes, 1983; Izhikevich, 2007; Strogatz, 2001). In Fig. 3(d), we see that, given sufficient noise amplitude, an intrinsically non-bursting neuron can be perturbed out of the V-nullcline region in the presence of noise (dotted elliptical region), thus resulting in a quick second spike, i.e. creating a spike doublet (black trajectory). We would have expected this neuron to return to its refractory state after a single spike in the absence of this noise perturbation (grey dashed
K. Wong-Lin et al. / Neural Networks (
)
–
5
Fig. 3. Neuronal spiking properties and mechanisms. (a) Single neuronal model with c = −60 mV can show from spike singlets to doublets as the applied current (mean or noise) increases. Horizontal axis: time from current onset. A Gaussian noise with standard deviation 2.5 is used to generate noise-induced bursts e.g. spike doublets (labelled by asterisks). To be compared with (d). (b) Firing rate adaptation over time. Instantaneous firing rates calculated from the inverse of the interspike intervals over different time points. Initial instantaneous firing rates (calculate from the first 2 spikes) increase faster than the steady-state values as applied current increases. c = −60 mV. (c) Current-vs-frequency (input–output) curve under deterministic and noisy conditions. Noise can permit low (subcritical) firing rates. Dotted line: averaged over 30 trial simulations with Inoise = 0.2. (d) Fluctuation-driven bursts can be generated by noise even when the neuronal properties are not intrinsically bursting (with low reset c values). Phase plane of model superimposed with the same temporal evolution of membrane potential from bottom panel of panel (a). Grey (black) part of a single trajectory: periodic regular spike singlets (occasional noise-driven spike doublet). Although noise is more visible along the V-nullcline, it is added at every time step of the trajectory, which can occasionally drive the trajectory outside the V-nullcline resulting in a burst of spike doublet (elliptical region). Inset: zoomed in to the elliptical region; arrow denote point of escaping out of the V-nullcline.
trajectory). A similar argument for burst with more than two spikes can be made. 3.2. Neuronal network model during a reward-based memory-guided saccade task The above discussion on neuronal spiking mechanism would allow us to generate heterogeneous neuronal spiking patterns. We will use a network model to simulate DRN neuronal activity of non-human primates performing a simple perceptual decision task for both rewarding and unrewarding trials (Bromberg-Martin et al., 2010; Nakamura et al., 2008). This will then inform possible underlying neuronal circuit architecture, afferent inputs to the DRN circuit, and predict plausible responses from the inhibitory neurons. Before that, we shall first discuss the general results in the experiments of Bromberg-Martin et al. (2010) and Nakamura et al. (2008). 3.2.1. Main experimental results In Bromberg-Martin et al. (2010), three interesting types of DRN neurons were identified. One type of (putative 5-HT) DRN neurons elevated its activity during the reward epoch of the trial (Fig. 4, top, light grey), but not in an unrewarding trial (Fig. 4, top, dark grey). These neurons seemingly encode (from their elevated activities) rewards from target onset through the delay period, and especially towards the end of the trial from the Go cue through the rewarding epoch of the trial. Another group of (putative 5-HT) DRN neurons showed the reverse pattern: elevated activity during unrewarding trial, but suppressed activity during rewarding trial during the target-delay and outcome epochs. Fig. 4 shows that these two types of neurons showed firing rate activities in a reciprocal way.
A third group of (putative) non-5-HT neurons showed only phasic activation for any stimulus (fixation point or target) onset (not shown, but see Bromberg-Martin et al., 2010). Thus, their work shows that DRN neurons that are tonically excited during the task are also excited during the rewarding epoch of the trial, while DRN neurons that are tonically suppressed below baseline during the beginning of fixation onset, are excited during the unrewarding epoch of the trial. This latter group of negative or no reward encoding neurons could play an important role in encoding some form of ‘‘disappointment’’ to a lack of reward, which could be broadcasted to other parts of the brain. This could in turn modulate mood, cognition and subsequent behaviours. 3.2.2. Network model simulations of behavioural task We shall now attempt to mimic the stimulus inputs as in the experiments (compare Fig. 5(a) with Fig. 2). To simulate the neuronal activities found in the experiments with the simplest possible network architecture, we are constrained to implement an inhibitory feedforward network architecture (Fig. 5(b)). If we assume inhibition from inhibitory neurons to 5-HT neurons to be, on average, equal, we found that the fixation input current has to be much smaller for non-reward-encoding neurons than for rewardencoding neurons. The feedforward inhibition from the non-5-HT inhibitory neurons to the 5-HT neurons also has to be about 2.5 times that of inhibition among the inhibitory neurons. Satisfying these requirements allows the network model to replicate the essential neural activity patterns of the two types of DRN neurons found in the Bromberg-Martin et al. (2010) experiment described in 3.2.1 (compare Fig. 4 (top panel) with Fig. 5(c), and Fig. 4 (bottom panel) with Fig. 5(d)). Additionally, although our implemented input currents are step functions
6
K. Wong-Lin et al. / Neural Networks (
)
–
Fig. 4. Neuronal activities in the DRN recorded of behaving monkeys. Two types of neurons found in the study. Top (bottom): a sample (non-)reward-encoding neuronal spiking frequency averaged over trials. Light grey (dark grey): (un)rewarded trial. Double asterisks show significant changes in neuronal activities between the two trial types. Source: Adapted from Bromberg-Martin et al. (2010) with permission.
Fig. 5. Spiking neuronal network model. (a) Time course of afferent inputs to 5-HT and non-5-HT neurons in the model. Time course implemented to mimic the memoryguided saccade task. mRT: mean reaction time. (b) A feedforward network architecture. Inh: non-5-HT inhibitory neurons. + and −: excitatory and inhibitory synapses. 5-HT neurons separated into reward-encoding and non-reward-encoding types. (c, d) Single-trial spike raster diagram of 5-HT neurons (neuron #1–200) and non-5-HT inhibitory neurons (neuron #201–300) from our network model. Non-reward-encoding neurons (neuron #1–100); reward-encoding 5-HT neurons (neuron #101–200). Memory-guided saccade task was modelled for a rewarding (c) and an unrewarding (d) trial. Labels (i) to (iv) denote the epochs in a trial as in (a).
(Fig. 5(a)), the neurons are capable of firing phasically locked to stimulus onset followed by adaptation, thus resembling the spiking
activities in the experiments (Fig. 4). Although the target input is implemented briefly for only 100 ms, the non-reward-encoding
K. Wong-Lin et al. / Neural Networks (
Fig. 6. Spontaneous network activity with bursting behaviour. The averaged population firing rate for 5-HT (black, neuron #1–200) and non-5-HT inhibitory (grey, neuron #201–300) neurons are ∼0.5 Hz and ∼5 Hz, respectively. I5-HT = 3.57 and Ii = 0.63. Squares denote bursts (two or more spikes within ∼10 ms).
neuronal activities continue to be reduced as in the experiment. This is due to the lower fixation input current than that of rewardencoding neurons. During the Go cue (disappearance of fixation point and where subjects are required to make a saccadic choice), we find in the model that the rewarding signal has to begin from the Go signal onset till about 500 ms after the outcome onset. Otherwise, due to the absence of the fixation input, the neuronal spiking activities will drop to baseline (not observed in the recorded neurons) and not sustained over time (data not shown). This may suggest that the Go signal during a rewarding block (of trials) might actually be associated with the rewarding outcome on the reward-encoding neurons, perhaps through reinforcement (e.g. temporal difference) learning during the training phase of the subjects. Interestingly, the non-rewarding outcome input has to appear after a much longer delay than that of the rewarding outcome input. This means that non-rewarding outcome may not be directly attached to the immediate decision. As suggested in BrombergMartin et al. (2010), during this epoch, an afferent input may come from other brain regions, e.g. the lateral habenula, which seem to encode negative reward (Hikosaka, 2010). Lastly, we observed in our simulations that the non-5-HT inhibitory neurons can exhibit strong oscillations (∼12.5 Hz) especially during the rewarding outcome epoch (Fig. 5(c)). This distinct oscillatory behaviour may be due to the overall higher total input current into these neurons, thus creating more synchrony. 3.2.3. Spontaneous bursting in the network model When we remove all behavioural task stimuli, the 5-HT neuronal population in our network model can attain an averaged spontaneous firing rate of ∼0.5 Hz for the 5-HT neurons, comparable to the base values found in experiments (Hajós et al., 1996; Vandermaelen & Aghajanian, 1983). The non-5-HT inhibitory neuronal population typically can have a higher firing rate of ∼7 Hz, which is comparable to the ∼5 Hz in our model. In Fig. 6, we can clearly see this difference in spike frequencies between the 5-HT and non-5-HT neurons. We have allowed the model to exhibit occasional bursts over time by having noise (varying afferent input) and heterogeneity (varying neuronal parameter c) in the system. Because the bursting behaviour is an intrinsic property of the neurons, we can observe that a neuron that has produced a burst at an earlier time will have a higher chance of bursting again at a later time (Fig. 6, squares). 3.2.4. Effects of inhibitory synapses in the behavioural task Inhibitory synapses are known to have important effects on the DRN circuit (e.g. Tao & Auerbach, 2003). Thus they are also the target of some experimental drug studies. Studies have also shown the general importance of inhibitory neurons and synapses on oscillatory network behaviour (e.g. Whittington & Traub, 2003). Therefore, here, we shall investigate, in the context of our network
)
–
7
model, how the strengths or weights of the synapses can affect the network dynamics in this specific behavioural task paradigm. We will focus only on rewarding trials here as the non-rewarding trials show similar effects. We found that as the synaptic weights, gi−5-HT , increases from 5 to the more extreme value 50, 5-HT neurons are almost fully suppressed, except during the phasic bursts upon the onsets of stimuli, where the 5-HT neurons are synchronized within their sub-population. Furthermore, immediately after the rewarding epoch, the reward-encoding 5-HT neurons fire synchronously with the non-reward-encoding 5-HT neurons (Fig. 7(a)). This could be due to the sudden relaxation of inhibitory neuronal activity, having a huge simultaneous effect on all the 5-HT neurons. Thus, synchronous firing can spread from a local sub-population of 5-HT neurons to more global effects on all the 5-HT neurons within a single trial. When the synaptic weight, gi−i , is instead increased from 2 to 50, there is no dramatic change to the firing rates (compare Fig. 7(b) with Fig. 5(c)). However, at such high gi−i value there is some subtle increase in synchronous firing among the inhibitory neurons, resulting in slightly less suppression from the inhibitory neurons to the non-reward-encoding 5-HT neurons. 3.3. Synaptic strengths and afferent currents on population network firing rate and oscillation frequency 3.3.1. Population firing rates Having observed some possible network behavioural changes above, we shall now perform a more systematic study of the network effects due to changes in the afferent current amplitudes and synaptic weights. Unless specified to be varied, the parameter values gi−i , gi−5-HT , I5-HT , and Ii are fixed at 2, 5, 3.57 and 0.63, respectively, for studying changes in the population firing rates. We first vary the input current into the inhibitory neurons, Ii , and found that the gain or slope of the population current vs firing rate function systematically changes (Fig. 8(a)). We can understand this in a heuristic way. Suppose the inhibitory population firing rate at steady state can be described by ri = F (−gi−i ri + Ib ), where Ib is some activation threshold or bias current, and that F is approximately linear. Then effectively, ri = Ib ∗ gain where gain = 1/(1 + gi−i ). Thus increasing gi−i can decrease the overall gain of the inhibitory population. When gi−5-HT is varied, we do not observe similar change in the input–output gain of the 5-HT neuronal population (Fig. 8(b)), as changing the gi−5-HT is equivalent to changing the afferent inputs into the 5-HT neurons. Comparing the population input–output functions of the inhibitory neurons compared to the 5-HT neurons, we can see that the activation threshold (current) is lower while the gain is higher for the former than that of the latter. Thus, the inhibitory neurons in the model are more easily excited than the 5-HT neurons, consistent with experiments. Next, we gradually increase the synaptic weight gi−i . As expected, the overall inhibitory neuronal population firing rate generally drops as self-inhibition within the population is enhanced, and this is not affected by gi−5-HT since there is no feedback connection in this simplified model architecture (Fig. 8(c)). It is, however, interesting to note that the inhibitory population firing rate does not decrease too much as previously seen in Fig. 7(b). Increasing gi−5-HT produces a similar trend for the 5-HT neuronal population firing rate, but the change is now more dramatic (Fig. 8(d)). This change is attenuated with higher gi−i . This is due to the slightly lower firing rate of the inhibitory neurons (Fig. 8(c)), and thus there will be weaker inhibition from the inhibitory neurons to the 5-HT neurons. Note that the gradual saturation of the change in Fig. 8(d) as gi−i increases is directly related to the saturating change in Fig. 8(c).
8
K. Wong-Lin et al. / Neural Networks (
)
–
Fig. 7. Synchronous activities in the network model. Raster diagram with (a) gi−i = 2 and gi−5-HT = 50; (b) gi−i = 50 and gi−5-HT = 5. Label of time epoch as in Fig. 5(c).
Fig. 8. Change in 5-HT and non-5-HT neuronal population firing rates with various afferent currents (a–b) and synaptic weights (c–d).
3.3.2. Population oscillation frequency We shall now investigate how the 5-HT neuronal population oscillation frequency can change with different synaptic strengths. We first sum the membrane potentials of all the 200 5-HT neurons to mimic the local field potential of the 5-HT neuronal population. We then identify distinctive peak(s) in the power spectrum. To study network oscillations, we raise the afferent currents I5-HT , and Ii to 4.25 and 1.1, respectively, so that the 5-HT neuronal population can now oscillate about 4 Hz. We settle for a ∼4 Hz oscillation frequency (Fig. 9(a)) to constrain it to the typical range of ∼0.5–4 Hz (e.g. (Allers & Sharp, 2003; Hajós et al., 1996; Sprouse & Aghajanian, 1987). Also, there is some experimental evidence of theta band oscillation found in 5-HT neurons in the raphe nuclei (Kocsis et al., 2006; Kocsis & Vertes, 1992). The power spectrum in Fig. 9(b) clearly shows the distinct peak around this value, verifying our technique. When gi−5-HT synaptic strength is increased, multi-stable oscillation frequencies starts to appear (Fig. 9(c) and (d)). The
double peaks, one at about 3.9 Hz and the other at 12.5 Hz, appearing in the power spectrum in Fig. 9(d) confirm this. The 12.5 Hz higher frequency oscillation does not seem to be caused by single neuronal sub-threshold oscillation, but appear when we summed the membrane potentials together (Fig. 9(c) insets). As gi−5-HT increases to some critical value (∼20), the lower 3.9 Hz frequency decreases slightly to 3.1 Hz, while the higher frequency abruptly appears (Fig. 9(e)). If the synaptic strength gi−i is higher (20), the previous 12.5 Hz frequency is moved to an even higher value (∼37 Hz) without changing the lower frequency much. In Fig. 9(e), we can also see that the critical gi−5-HT is right-shifted. This phenomenon could be understood if we consider the gi−5-HT as a ‘‘gate’’ for oscillatory signals from the inhibitory population to influence the subthreshold membrane potentials of the 5-HT population. At certain critical gi−5-HT value, the gate is opened for information to flow from one population to another. This critical value has to be increased with higher gi−i as the lower inhibitory population activity is more suppressed, and a higher ‘‘gain’’ is needed to transmit the oscillatory information.
K. Wong-Lin et al. / Neural Networks (
)
–
9
Fig. 9. 5-HT population oscillation frequencies. (a,b) Control condition with about 4 Hz oscillation frequency in the summed membrane potentials (simulating local field potential). (c, d) Increased gi−5-HT synaptic strength results in multi-stable frequency. The higher frequencies do not come from single neuronal subthreshold oscillation (left inset in (c), single neuronal membrane potential timecourse) but appear when summed (right inset in (c), with 10 5-HT neurons). (d) Two peaks in the power spectrum (compared to one peak in (b)), one at 3.9 Hz and the other at 12.5 Hz. (b, d) A threshold of 55 dB is used to select the more distinctive oscillation frequencies. (e) Summary of the range of oscillation frequencies in the summed membrane potentials of 5-HT neurons with increasing gi−5-HT synaptic strength. With gi−i = 2 and gi−5-HT = 20, higher frequency appears at ∼12.5 Hz. With a higher gi−i value (e.g. of 20), higher oscillation frequencies can be obtained (e.g. ∼37 Hz, grey).
4. Discussion We have developed an efficient spiking neuronal network model of the DRN which consists of coupled 5-HT and non-5HT inhibitory neurons. At the neuronal level, we have shown, using computational simulations and dynamical systems analysis, that with variability in a slow recovery current, overall afferent input, and intrinsic neuronal noise, we can efficiently emulate heterogeneous spiking patterns (regular spiking, bursting or mixture of both). This is achieved even when the neurons have similar basic electrophysiological properties. These neuronal properties are actually consistent with experimental findings. With these constraints, we coupled the heterogeneous and noisy 5-HT-containing and non-5HT-containing inhibitory neurons to construct an important part of the DRN circuit network. The circuit can produce spontaneous bursting behaviour that is inherently a property of the neurons; neurons within a network that exhibit burst at some point will have a higher tendency to show bursting again at a later time. We constrain our network model architecture to a feedforward type so that the neuronal data from the Bromberg-Martin et al. (2010) experiment can be replicated. This suggests that inputs activating the inhibitory neurons may play an important role in providing fast inhibition to ‘‘competing’’ reward-encoding and non-reward-encoding 5-HT neurons. Input due to rewarding outcome (encoded by some reward-selective 5-HT neurons) seems to be activated during the Go cue rather than around the time of the actual reward, possibly suggesting some form of temporal difference learning is involved. Non-rewarding outcome (encoded by some non-reward-selective 5-HT neurons) appears after some significant delay, which suggests influence through different neural pathways, possibly via the lateral habenula. Increased synaptic strength of gi−5-HT not only suppresses 5-HT neuronal
spiking activity, but also enhances their tendency to synchronize, from local sub-population to globally. It can also decrease the gain of the inhibitory population input–output function. However, an increased self-inhibitory synaptic strength of gi−i does not seem to change the network as much, showing some form of robustness within the inhibitory population. We have seen how with the behavioural task stimuli that the network can show theta band oscillation especially around the rewarding outcome epoch of the trial. By increasing the constant bias input currents to the neurons, oscillations in the local field potentials of the network can appear consistently. We found that under certain parameter regimes, multiple oscillation frequencies can co-exist—multi-stability. Interestingly, under this regime of high gi−5-HT , the lower stable frequencies (∼2–4 Hz) seem to be related to the neuronal spike timings while the higher frequencies (∼12–40 Hz) of the 5-HT population’s local field potential is due to the transmission of the inhibitory population’s oscillatory signals. Theta oscillations of 5-HT and non-5-HT neurons in the raphe nuclei have been found in experiments (Kocsis et al., 2006; Kocsis & Vertes, 1992). This theta oscillation in the raphe nuclei can be phase locked to the theta oscillation of other brain regions, e.g. the hippocampal EEG (Kocsis et al., 2006), and a recent work has even suggested that serotonin can control hippocampal theta rhythm (Sörman, Wang, Hajós, & Kocsis, 2011). Thus theta oscillation may be a way to transfer information between the raphe nuclei and other brain regions like the hippocampus, and so it may be important to know whether coupling with other brain regions is necessary to generate theta oscillation in the raphe nuclei. Our model suggests that the DRN is capable of generating its own intrinsic theta oscillations. Slower oscillations e.g. ∼1–2 Hz have been observed in some experiments (Allers & Sharp, 2003; Kocsis et al., 2006; Wang & Aghajanian, 1982). Our simple network model is able to show
10
K. Wong-Lin et al. / Neural Networks (
oscillation frequency of about 2 Hz when the synaptic strength gi−5-HT is sufficiently high (Fig. 9(e)). Our model predicts that around this regime, we may also expect to observe high frequency oscillations in the local field potential of the 5-HT neurons. This could be verified in experiments. However, we have been unable to find parameter regimes where our network model can oscillate at a frequency lower than 2 Hz. Thus very slow oscillations may possibly be due to the presence of some slow ionic currents (e.g. mediated by NMDA, GABAB , or calcium) not implemented in our current model. Thus, a possible direction for future computational work may be to implement more specific ionic channels especially those with slower dynamics, e.g. (Kocsis et al., 2006; Penington, Kelly, & Fox, 1991; Rouchet et al., 2008). To conclude, we have successfully constructed an efficient neuronal circuit model of the DRN that may allow convenient investigation of serotonergic modulation in large-scale brain circuits. Moreover, the same adaptive spiking neuronal model can simulate the spiking properties of cortical and subcortical neurons with minimal parameter changes, even with incomplete experimental data (e.g. specific ionic channel properties). Incorporating genetic and intracellular dynamical processes (Best et al., 2010; Stoltenberg & Nag, 2010) into such models would result in a highly unified and multi-scale computational model. Acknowledgements The authors would like to thank Shahjahan Shahid for advice on power spectral analysis of noisy oscillatory signals. This study was initially supported under the CNRT award by the Northern Ireland Department for Employment and Learning through its ‘‘Strengthening the All-Island Research Base’’ initiative, and more recently via the Centre of Excellence in Intelligent Systems award, funded by InvestNI and the Integrated Development Fund, through its local facilitator, ILEX. References Aghajanian, G. K., & Vandermaelen, C. P. (1982). Intracellular identification of central noradrenergic and serotonergic neurons by a new double labeling procedure. The Journal of Neuroscience, 2, 1786–1792. Allers, K. A., & Sharp, T. (2003). Neurochemical and anatomical identification of fastand slow-firing neurones in the rat dorsal raphe nucleus using juxtacellular labelling methods in vivo. Neuroscience, 122, 193–204. Baker, K. G., Halliday, G. M., Hornung, J. P., Geffen, L. B., Cotton, R. G., & Törk, I. (1991). Distribution, morphology and number of monoaminesynthesizing and substance P-containing neurons in the human dorsal raphe nucleus. Neuroscience, 42(3), 757–775. Beck, S. G., Pan, Y.-Z., Akanwa, A. C., & Kirby, L. G. (2004). Median and dorsal raphe neurons are not electrophysiologically identical. Journal of Neurophysiology, 91, 994–1005. Best, J., Nijhout, H. F., & Reed, M. (2010). Serotonin synthesis, release and reuptake in terminals: a mathematical model. Theoretical Biology and Medical Modelling, 7, 34. doi:10.1186/1742-4682-7-34. Best, J., Nijhout, H. F., & Reed, M. (2011). Bursts and the efficacy of selective serotonin reuptake inhibitors. Pharmacopsychiatry, 44, S76–S83. doi:10.1055/s0031-1273697. Bromberg-Martin, E. S., Hikosaka, O., & Nakamura, K. (2010). Coding of task reward value in the dorsal raphe nucleus. The Journal of Neuroscience, 30, 6262–6272. Calizo, L. H., Ma, X., Pan, Y., Lemos, J., Craige, C., Heemstra, L. A., et al. (2011). Raphe serotonin neurons are not homogenous: electrophysiological, morphological and neurochemical evidence. Neuropharmacology, 61, 524–543. Carr, V. G., & Lucki, I. (2011). The role of serotonin receptor subtypes in treating depression: a review of animal studies. Psychopharmacology (Berl), 213, 265–287. Crawford, L. K., Craige, C. P., & Beck, S. G. (2010). Increased intrinsic excitability of lateral wing serotonin neurons of the dorsal raphe: a mechanism for selective activation in stress circuits. Journal of Neurophysiology, 103, 2652–2663. Cryan, J. F., & Leonard, B. E. (2000). 5-HT1A and beyond: the role of serotonin and its receptors in depression and the antidepressant response. Human Psychopharmacology, 15, 113–135. Daw, N. D., Kakade, D., & Dayan, P. (2002). Opponent interactions between serotonin and dopamine. Neural Networks, 15, 603–616.
)
–
Dayan, P., & Huys, Q. J. (2009). Serotonin in affective control. Annual Review of Neuroscience, 32, 95–126. Doya, K. (2002). Metalearning and neuromodulation. Neural Networks, 15, 495–506. Giovanni, G. D., Matteo, V. D., & Esposito, E. (2008). Progr. brain res.: Vol. 172. Serotonin-dopamine interaction: experimental evidence and therapeutic relevance. Elsevier Science. Guckenheimer, J., & Holmes, P. (1983). Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. New York: Springer-Verlag. Hajós, M., Allers, K. A., Jennings, K., Sharp, T., Charette, G., Sík, A., et al. (2007). Neurochemical identification of stereotypic burst-firing neurons in the rat dorsal raphe nucleus using juxtacellular labelling methods. European Journal of Neuroscience, 25, 119–126. Hajós, M., Sharp, T., & Newberry, N. R. (1996). Intracellular recordings from burstfiring presumed serotonergic neurones in the rat dorsal raphe nucleus in vivo. Brain Research, 737, 308–312. Hikosaka, O. (2010). The habenula: from stress evasion to value-based decisionmaking. Nature Reviews Neuroscience, 11, 503–513. Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Transactions on Neural Networks, 14, 1569–1572. Izhikevich, E. M. (2007). Dynamical systems in neuroscience. Cambridge, MA: MIT Press. Jacobs, B. L., & Azmitia, E. C. (1992). Structure and function of the brain serotonin system. Physiological Reviews, 72, 165–229. Kirby, L. G., Pernar, L., Valentino, R. J., & Beck, S. G. (2003). Distinguishing characteristics of serotonin and non-serotonin-containing cells in the dorsal raphe nucleus: electrophysiological and immunohistochemical studies. Neuroscience, 116, 669–683. Kocsis, B., Varga, V., Dahan, L., & Sik, A. (2006). Serotonergic neuron diversity: identification of raphe neurons with discharges time-locked to the hippocampal theta rhythm. Proceedings of the National Academy of Sciences USA, 103, 1059–1064. Kocsis, B., & Vertes, R. P. (1992). Dorsal raphe neurons: synchronous discharge with the theta rhythm of the hippocampus in the freely behaving rat. Journal of Neurophysiology, 68, 1463–1467. Köhler, C., & Steinbusch, H. (1982). Identification of serotonin and non-serotonincontaining neurons of the mid-brain raphe projecting to the entorhinal area and the hippocampal formation. A combined immunohistochemical and fluorescent retrograde tracing study in the rat brain. Neuroscience, 7, 951–975. Lemos, J. C., Pan, Y.-Z., Ma, X., Lamy, C., Akanwa, A. C., & Beck, S. G. (2006). Selective 5-HT1B receptor inhibition of glutamatergic and GABAergic synaptic activity in the rat dorsal and median raphe. European Journal of Neuroscience, 24, 3415–3430. Li, Y.-Q., Li, H., Kaneko, T., & Mizuno, N. (2001). Morphological features and electrophysiological properties of serotonergic and non-serotonergic projection neurons in the dorsal raphe nucleus. An intracellular recording and labeling study in rat brain slices. Brain Research, 900, 110–118. Marinelli, S., Schnell, S. A., Hack, S. P., Christie, M. J., Wessendorf, M. W., & Vaughan, C. W. (2004). Serotonergic and nonserotonergic dorsal raphe neurons are pharmacologically and electrophysiologically heterogeneous. Journal of Neurophysiology, 92, 3532–3537. Nakamura, K., Matsumoto, M., & Hikosaka, O. (2008). Reward-dependent modulation of neuronal activity in the primate dorsal raphe nucleus. The Journal of Neuroscience, 28, 5331–5343. Pacher, P., & Kecskemeti, V. (2004). Trends in the development of new antidepressants. Is there a light at the end of the tunnel? Current Medicinal Chemistry, 11, 925–943. Penington, N. J., Kelly, J. S., & Fox, A. P. (1991). A study of the mechanism of Ca2+ current inhibition produced by serotonin in rat dorsal raphe neurons. The Journal of Neuroscience, 11, 3594–3609. Ranade, S. P., & Mainen, Z. F. (2009). Transient firing of dorsal raphe neurons encodes diverse and specific sensory, motor, and reward events. Journal of Neurophysiology, 102, 3026–3037. Rouchet, N., Waroux, O., Lamy, C., Massotte, L., Scuvée-Moreau, J., Liégeois, J. F., et al. (2008). SK channel blockade promotes burst firing in dorsal raphe serotonergic neurons. European Journal of Neuroscience, 28, 1108–1115. Sörman, E., Wang, D., Hajós, M., & Kocsis, B. (2011). Control of hippocampal theta rhythm by serotonin: role of 5-HT2c receptors. Neuropharmacology, 61, 489–494. Sprouse, J. S., & Aghajanian, G. K. (1987). Electrophysiological responses of serotoninergic dorsal raphe neurons to 5-HT1A and 5-HT1B agonists. Synapse, 1, 582–588. Stoltenberg, S. F., & Nag, P. (2010). Description and validation of a dynamical systems model of presynaptic serotonin function: genetic variation, brain activation and impulsivity. Behavior Genetics, 40, 262–279. Strogatz, S. H. (2001). Nonlinear dynamics and chaos: with applications to physics, biology, chemistry and engineering. Cambridge, MA: Perseus Books Group. Tao, R., & Auerbach, S. B. (2003). Influence of inhibitory and excitatory inputs on serotonin efflux differs in the dorsal raphe and median raphe nuclei. Brain Research, 961, 109–120. Tseng, K.-Y., & Atzori, M. (2007). Monoaminergic modulation of cortical excitability. Springer.
K. Wong-Lin et al. / Neural Networks ( Vandermaelen, C. P., & Aghajanian, G. K. (1983). Electrophysiological and pharmacological characterization of serotonergic dorsal raphe neurons recorded extracellularly and intracellularly in rat brain slices. Brain Research, 289, 109–119. Wang, R. Y., & Aghajanian, G. K. (1982). Correlative firing patterns of serotonergic neurons in rat dorsal raphe nucleus. The Journal of Neuroscience, 2, 11–16.
)
–
11
Whittington, M. A., & Traub, R. D. (2003). Interneuron diversity series: inhibitory interneurons and network oscillations in vitro. Trends in Neurosciences, 26, 676–682. Wong-Lin, K. F., Prasad, G., & McGinnity, T. M. (2011). A spiking neuronal network model of the dorsal raphe nucleus. In Proceedings of 2011 international joint conference on neural networks. San Jose, CA, USA.
