Temporal derivative computation in the dorsal raphe network revealed by an experimentally driven augmented integrate-and-fire modeling framework

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

955 views

Version of Record published: February 27, 2023 (This version)
Accepted Manuscript published: January 19, 2023 (Go to version)
Accepted: December 19, 2022
Received: August 10, 2021
Preprint posted: June 27, 2021 (Go to version)

1. Of interest
Systems level identification of a matrisome-associated macrophage polarisation state in multi-organ fibrosis

John F Ouyang, Kunal Mishra ... Jacques Behmoaras

Research Article Oct 3, 2023
Further reading

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Appendix 1
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

By means of an expansive innervation, the serotonin (5-HT) neurons of the dorsal raphe nucleus (DRN) are positioned to enact coordinated modulation of circuits distributed across the entire brain in order to adaptively regulate behavior. Yet the network computations that emerge from the excitability and connectivity features of the DRN are still poorly understood. To gain insight into these computations, we began by carrying out a detailed electrophysiological characterization of genetically identified mouse 5-HT and somatostatin (SOM) neurons. We next developed a single-neuron modeling framework that combines the realism of Hodgkin-Huxley models with the simplicity and predictive power of generalized integrate-and-fire models. We found that feedforward inhibition of 5-HT neurons by heterogeneous SOM neurons implemented divisive inhibition, while endocannabinoid-mediated modulation of excitatory drive to the DRN increased the gain of 5-HT output. Our most striking finding was that the output of the DRN encodes a mixture of the intensity and temporal derivative of its input, and that the temporal derivative component dominates this mixture precisely when the input is increasing rapidly. This network computation primarily emerged from prominent adaptation mechanisms found in 5-HT neurons, including a previously undescribed dynamic threshold. By applying a bottom-up neural network modeling approach, our results suggest that the DRN is particularly apt to encode input changes over short timescales, reflecting one of the salient emerging computations that dominate its output to regulate behavior.

Editor's evaluation

To characterize physiological properties of dorsal raphe serotonin neurons, the authors applied the approach called an augmented generalized integrate-and-fire [aGIF] model, which incorporates a relatively small number of salient biophysical properties of a specific neuron type, and whose parameters are optimized based on voltage dynamics obtained experimentally. The results showed that after-hyperpolarization and A-type potassium currents, in combination with heterogeneous feedforward inhibition from local GABA neurons, give rise to a derivative-like input-output relationship in serotonin neurons.

https://doi.org/10.7554/eLife.72951.sa0

Introduction

The forebrain-projecting serotonin (5-HT) neurons of the dorsal raphe nucleus (DRN) play a key role in regulating behavior in dynamic environments, but the precise nature of this role is still not well understood (Young et al., 1985; Delgado, 1994; Warden et al., 2012; Dayan and Huys, 2015). DRN serotonin neurons have been proposed to modulate a wide range of cognitive processes, such as encouraging patience for future rewards (Miyazaki et al., 2014; Fonseca et al., 2015), signaling the beneficialness of current actions or states (Luo et al., 2016), complementing reinforcement signals of dopamine (Daw et al., 2002; Maier and Watkins, 2005; Nakamura et al., 2008; Ranade and Mainen, 2009; Tops et al., 2009; Cools et al., 2011; Li et al., 2016), and, partially as a corollary, regulating both learning (Soubrié, 1986; Deakin, 1991; Daw et al., 2002; Dayan and Huys, 2009; Matias et al., 2017; Grossman et al., 2022) and mood (Savitz et al., 2009; Fava and Kendler, 2000; Donaldson et al., 2013; Cipriani et al., 2018). While the remarkable diversity of roles attributed to this single neurotransmitter has historically been perplexing, recent findings are beginning to provide insight (see Okaty et al., 2019 for review). For example, the unsuspected organization of 5-HT neurons into anatomical sub-modules that differentially regulate behavior (Abrams et al., 2004; Lowry et al., 2005; Commons, 2015; Muzerelle et al., 2016; Ren et al., 2018; Paquelet et al., 2022), or the observation that 5-HT neurons can encode distinct salient features of the environment over different timescales (Trulson and Jacobs, 1979; Schweimer and Ungless, 2010; Ranade and Mainen, 2009; Cohen et al., 2015; Zhong et al., 2017) is a compelling mechanism that may contribute to the multiplicity of 5-HT’s actions. These anatomical and dynamical perspectives on 5-HT diversity need not be mutually exclusive. A clearer understanding of the biophysical mechanisms that contribute to the coding features of raphe neurons over multiple timescales has the potential to substantially increase our understanding of how 5-HT regulates behavior.

The spiking statistics of 5-HT neurons necessarily shape and constrain their computational role. For instance, the slow firing rate (~5 Hz) of 5-HT neurons, in large part attributable to a large after-hyperpolarization potential (AHP) (Aghajanian and Vandermaelen, 1982; Vandermaelen and Aghajanian, 1983), may appear to preclude signaling on faster timescales. However, fast signaling despite slow firing can arise naturally in ensemble-rate codes (Knight, 1972; Gerstner, 2000). Consistent with this idea, the in vivo population activity of 5-HT neurons has been observed to track impending rewards over second to sub-second timescales (Zhong et al., 2017), and the trial-averaged ensemble rates of individual 5-HT neurons can track environmental changes over the millisecond timescale (Ranade and Mainen, 2009; Cohen et al., 2015). In addition, the fact that 5-HT receptor subtypes can regulate the excitability of target neurons over different timescales, including ionotropic 5-HT3 receptors with millisecond gating kinetics (Béïque et al., 2004; Béïque et al., 2007; Andrade, 2011; Varga et al., 2009), at the very least suggests that the 5-HT system is capable of fast information transmission, an observation mirrored by the fast dynamics of neurons which project to the DRN (Amo et al., 2014; Matsumoto and Hikosaka, 2007). If fast and slow signaling by the DRN are manifest, it is less clear which cellular mechanisms regulate the interplay between these timescales, nor which input features are represented on which timescales.

Computational modeling is a standard approach to link levels of description and is thus well suited to delineate how network-level function emerges from excitability features identified at the single-cell level. In spite of their conceptual utility, the most detailed single cell models, including those of DRN neurons (Tuckwell and Penington, 2014; Wong-Lin et al., 2011), do not lend themselves with ease to bottom-up modeling efforts because of the substantial technical difficulty of obtaining sufficiently accurate values for a large number of interacting model parameters (Prinz et al., 2004; Gerstner and Naud, 2009). Mathematically simpler generalized integrate-and-fire (GIF) models provide a strong foundation for network modeling because their small number of parameters can be estimated with a high degree of precision (Mensi et al., 2012; Pozzorini et al., 2013; Teeter et al., 2018). This precision comes at a price, however the process of distilling the effects of many biophysical mechanisms into a small number of model parameters makes it difficult to study a specific mechanism (e.g. a subthreshold ion channel) in isolation. A hybrid approach based on a reductionistic GIF model augmented with a limited set of biophysical mechanisms could leverage the precision of GIFs while allowing the ability to link specific biophysical mechanisms with higher-order network function.

In this study, we developed and validated for DRN neurons a hybrid modeling approach that lies between reductionist GIF and biophysical Hodgkin-Huxley-type models to capture excitability features of individual neurons for accurate simulations of population dynamics and, by extension, network computation inference. To this end, we carried out cellular electrophysiological recordings from genetically identified DRN 5-HT and SOM neurons to (1) extract and validate, from sets of noisy inputs, parameters for the automatic development of accurate GIF models and (2) experimentally define complementary biophysical mechanisms to be grafted onto the GIF models to iteratively improve their prediction accuracy (augmented GIFs). This approach recapitulated and extended past findings on DRN neurons by showing that the best-performing models of 5-HT neurons featured slow membrane time constants, an A-type potassium current, and strong adaptation mechanisms. Network simulations of optimized GIF models of both 5-HT and GABAergic SOM neurons organized in a feed-forward inhibitory circuit revealed that 5-HT neuron populations context-dependently encode a mixture of the intensity and temporal derivative of their inputs. Our overall approach further allowed us to trace back specific features of these population responses (e.g. gain) to defined excitability features of DRN neurons.

Results

Salient electrophysiological features of DRN neurons

Our main goal was to develop an experimentally grounded model of the DRN to better understand its computational properties. As a first step toward this goal, we carried out experiments to constrain a set of single-neuron models of the two main cell types found in the DRN: 5-HT and SOM GABA neurons. We performed whole-cell electrophysiological recordings from genetically identified 5-HT (Figure 1A1; SERT-Cre::Rosa-TdTomato mice) and SOM (Figure 1A2; Table 1; SOM-Cre::Rosa-TdTomato mice) neurons in slices. In keeping with previous descriptions (e.g. Vandermaelen and Aghajanian, 1983; Calizo et al., 2011), in the majority of the 5-HT neurons recorded in our dataset, current steps induced strongly adapting action potential firing accompanied by large AHPs, and a characteristic kink in the voltage trace leading up to the first spike (Figure 1B). Qualitatively distinct firing patterns of 5-HT neurons were, however, occasionally observed (Figure 1—figure supplement 1). Recordings from SOM neurons revealed spiking patterns that were more heterogeneous (Figure 1—figure supplement 2). Comparing the relationship between the injected currents and firing frequencies between these populations, we found that SOM neurons were generally more sensitive to changes in input current (gain) than 5-HT neurons and responded to weaker inputs (Figure 1B). The gain showed greater variability in SOM neurons than in 5-HT neurons (Brown-Forsythe equality of variance test p=0.001 on N = 17 5-HT and N = 7 SOM neurons). In line with this observation, SOM neurons also consistently exhibited a wider range of firing frequencies for a given input (e.g. for a 50 pA input 5-HT neurons fired at 2.81 ± 2.22 Hz vs 8.16 ± 5.70 Hz; Brown-Forsythe test p=0.005 in $N = 17$ 5-HT neurons and $N = 14$ SOM cells). Together, these observations outlined three salient cellular-level features of DRN neurons, namely the strong AHP and voltage kink of 5-HT neurons as well as noticeable heterogeneous excitability of SOM neurons.

Figure 1 with 7 supplements see all

Download asset Open asset

Physiology of dorsal raphe nucleus (DRN) neurons.

(**A, B**) Morphology, current steps (A), and spike frequency vs. input ( $f / I$ ) curves (B) of genetically identified DRN neurons. Non-monotonically increasing $f / I$ curves (gray) and linear fits to monotonically increasing curves (red dashed) for $N = 17$ serotonin (5-HT) (B1) and $N = 14$ somatostatin (SOM) (B2) neurons. (B3) Rheobase of 33.8 ± 21.0 pA in 5-HT neurons vs. 11.3 ± 16.0 pA in SOM neurons. (B4) Gain of 52.2 ± 22.2 Hz/nA in 5-HT neurons vs. 87.2 ± 33.0 Hz/nA in $N = 7$ SOM neurons with monotonically increasing $f / I$ curves. (C) Leak-subtracted whole-cell currents evoked by a depolarizing step. Each trace is one cell; $N = 13$ 5-HT and $N = 11$ SOM cells. Traces without a transient outward current are shown in gray. (D) Proportion of neurons with a transient outward current by cell type. (E) Quantification of transient outward currents in each cell type. $N = 3$ SOM cells without a transient outward current were excluded from analysis, leaving $N = 13$ 5-HT and $N = 8$ SOM neurons. Annotations reflect Mann-Whitney U-tests. Non-parametric Brown-Forsythe equality of variance tests indicated significantly more variable time to peak (p=1.11e-4; E2) and inactivation time constant (p=1.97e-4; E3) in SOM cells.

Table 1

Membrane parameters of DRN neurons.

Parameters obtained from recordings from mPFC L5 pyramidal neurons used to fit GIF models as a point of comparison are also shown. Data are presented as mean ± SD. Distributions are shown in Figure 1—figure supplement 3.

Cell type	R (GΩ)	C (pF)	τ (ms)	N
5-HT	1.16±0.55	67.0±17.1	75.2±33.8	96
SOM	1.07±0.58	43.5±15.5	42.2±19.8	28
mPFC	0.188±0.130	160.6±48.2	27.4±13.2	25

DRN, dorsal raphe nucleus; PFC, prefrontal cortex; GIF, generalized integrate-and-fire model; 5-HT, serotonin; SOM, somatostatin; mPFC, medial prefrontal cortex; R, memrane resistance; C, membrane capacitance; τ, membrane time constant.

The characteristic kink in the voltage leading up to the first spike in 5-HT neurons in principle may be caused by near-threshold activation of voltage-gated potassium channels (VGKCs; Connor and Stevens, 1971; Connor et al., 1977; Drion et al., 2015). We therefore examined whole-cell currents evoked by voltage steps (from –90 mV to –20 mV) in both 5-HT and SOM neurons to look for evidence of such a VGKC. In 5-HT cells, these experiments revealed a large (peak amplitude 928 ± 249 pA, leak-subtracted), partly inactivating (steady-state amplitude 142 ± 45 pA, leak-subtracted) outward current (Figure 1C1) that was sensitive to K_v4-selective potassium channel blockers (Figure 1—figure supplement 4). This current activated rapidly (peak latency 7.46 ± 0.21 ms) and inactivated over tens of milliseconds (inactivation time constant $τ_{h} = 42.9 \pm 9.4$ ms; kinetics are similar at near-physiological temperature, see Figure 1—figure supplement 5). The gating and kinetic profile (Table 2, Figure 1—figure supplements 6 and 7) of the inactivating component of this conductance in 5-HT neurons are broadly similar to those expected of the A-type potassium currents (I_A) characterized in great detail in several other cell types (e.g. Storm, 1989). Because these parameters are sufficient to construct a model of this conductance (see below), we have not attempted to determine its molecular identity further. For the sake of simplicity, we refer to the inactivating component herein as I_A (in keeping with the previous literature; see Aghajanian, 1985; Tuckwell and Penington, 2014) and the steady-state component as $I_{K}$ . Thus, an I_A -like inactivating VGKC is a consistent feature of DRN 5-HT neurons.

Table 2

5-HT I_A current gating parameters.

Gating curves shown in Figure 1—figure supplement 6B were fitted with the scaled Boltzmann function $g_{\infty} (V) / g_{\infty} (V_{r e f}) = x_{\infty} = A_{x} / (1 + \exp [- k_{x} (V - V_{x}^{*})])$ . Values are based on experiments from N=13 cells.

Gate	$V_{r e f} (m v)$	$A_{x}$	$k_{x} (m V^{- 1})$	$V_{x}^{*} (m V)$
$m_{\infty}$	–20	1.61	0.0985	–23.7
$h_{\infty}$	–80	1.03	–0.165	–59.2
$n_{\infty}$	–20	1.55	0.216	–24.3

5-HT, serotonin.

The same voltage-clamp protocol applied to SOM neurons, in contrast, triggered a mixture of outward and inward currents that varied widely from cell to cell (Figure 1D2). A significant proportion of SOM neurons did not express a transient outward current at all (27.3 %, Figure 1E), while the remaining cells had currents that were significantly smaller ( $p = 0.003$ ), activated more slowly ( $p = 0.003$ ), and exhibited much more heterogeneous kinetic profiles than those found in 5-HT neurons (Figure 1E3). Together, these results show that the expression of this subthreshold voltage-gated current is substantially more variable in SOM neurons than in 5-HT neurons, in line with the distinctive heterogeneity of excitability features observed in this DRN cell type (Figure 1C2, Figure 1—figure supplement 2).

I_A regulates initial firing rate via a control of spike time jitter

To develop an intuition for how I_A impacts the firing patterns of 5-HT populations, we first created a toy leaky integrate-and-fire (LIF) model that captured the effect of this conductance on single-cell voltage dynamics (see Methods). In keeping with previous studies, I_A introduced a kink in the subthreshold voltage leading up to spike threshold (Getting, 1983; Segal, 1985; McCormick, 1991; Figure 2A) and increased the latency to the first spike evoked by a square step stimulus, particularly when starting from a hyperpolarized voltage at which I_A is free from inactivation (Figure 2B and C). The effect of I_A on spike latency depends at least to some extent on its effective magnitude and inactivation kinetics (defined as the ratio of maximal A-type conductance to inverse membrane resistance and the ratio of the inactivation time constant of I_A to the membrane time constant; see Methods). When we set the corresponding parameters in our toy model to experimentally determined values from 5-HT neurons, we observed the same qualitative relationship between spike latency and initial voltage (Figure 2—figure supplement 1), further pointing toward a functional effect of I_A in this cell type. The predicted relationship between initial voltage and latency was experimentally recapitulated in whole-cell recordings from identified 5-HT neurons (Figure 2D–F). In particular, the onset of spiking was delayed by hyperpolarization (Figure 2D), and the magnitude of this effect was significantly reduced by the partial pharmacological block of I_A with 4-AP ( p=3.2e-6 for initial voltages between –90 mV and –70 mV; Figure 2E). Finally, we also observed an inflection point predicted by the model in the normalized initial-voltage/latency relationship (Figure 2F, compare with model prediction in Figure 2C). In summary, our toy model captured the expected effects of I_A in single 5-HT cells.

Figure 2 with 1 supplement see all

Download asset Open asset

I_A qualitatively alters the relationship between initial voltage and spike timing.

(**A–C**) A toy leaky integrate-and-fire (LIF) model neuron with I_A predicts a non-linear effect of voltage history on spike timing in a simple experiment. (**D–F**) Experiments in serotonin (5-HT) neurons fulfill predictions of the toy model. (F) Latency curves for $N = 8$ 5-HT neurons, normalized to the maximum latency for each cell. Each color is one cell. (**G–I**) I_A causes an increase in spike latency and jitter in the presence of noise. Models and input are the same as in A–C. Spike latency histograms for populations of 600 toy neurons are shown in I. The width of the histogram reflects jitter in the timing of the first spike, while the height of the histogram approximates the peak instantaneous firing rate. Note that as jitter increases, the height of the histogram decreases. The toy model with 100% I_A has an effective I_A conductance ${\bar{g}}_{A}^{'} / g_{l} = 10$ and effective inactivation time constant $τ_{h} / τ_{m e m} = 1.2$ .

Next, we used our experimentally validated toy model to understand how I_A impacts the spiking responses of whole neuronal populations. To do this, we simulated the effect of a shared step input to a population of 600 toy neurons each receiving independent background noise (corresponding to naturalistic fluctuations in synaptic inputs). Whereas subthreshold fluctuations yielded time-locked spikes without I_A (Figure 2G), they induced spiking with larger jitter across the simulated population when I_A was present (Figure 2H). This desynchronizing effect of I_A also decreased the peak population rate at the time corresponding to the mean latency (Figure 2I) since the peak rate corresponds to the coincidence rate from an ensemble of cells with similar properties. (The same effects were also observed in the toy models with parameters constrained to experimentally determined values; Figure 2—figure supplement 1.) Taken together, results from these toy models revealed a role of I_A in regulating the degree of synchronization of a population following sudden inputs, suggesting that I_A may regulate the gain of the DRN network to time-varying inputs. This intuition gleaned from this toy model is examined in more detail with optimized GIF models (see below).

Extensions to GIF models are required to capture the excitability of DRN neurons

We next sought to develop a model able to capture the essential biophysical features of DRN neurons and accurately predict their responses to naturalistic inputs. GIF models offer a flexible modeling framework well suited to this purpose because they can be trained to accurately reproduce the firing patterns of individual neurons using less than 5 min of electrophysiological data per neuron (Gerstner et al., 2014; Teeter et al., 2018; Paninski et al., 2005; Mensi et al., 2012; Pozzorini et al., 2015). In this framework, individual neurons are described in terms of three core components: (1) a passive membrane filter, $κ$ , which transforms input currents into a subthreshold membrane potential; (2) a stochastic spiking process, which transforms the subthreshold membrane potential into action potentials; and (3) two adaptation mechanisms, namely a spike-triggered current mediating the commonly observed AHP, $η$ , and change in firing threshold, $γ$ (Figure 3A1, see Methods). These components are described by parameters, the values of which are inferred from the electrophysiological data using a combination of least-squares multi-linear regression and gradient ascent of a likelihood function. The flexibility and data efficiency of this framework lend itself well to capturing the functional properties of single neurons and, by extension, heterogeneous neural populations.

Figure 3 with 3 supplements see all

Download asset Open asset

Single neuron models accurately predict the subthreshold voltage and spike times of dorsal raphe nucleus neurons.

(A) Components of candidate single neuron models: $λ$ intensity of stochastic spike-generating process; $η$ spike-triggered current (positive values indicate a hyperpolarizing current); $γ$ spike-triggered threshold movement. (B) Generalized integrate-and-fire (GIF) model extensions. (C) Representative experiment used to train and validate neuron models. Training set consists of repetitions of 60 s of frozen Ornstein-Uhlenbeck (OU) noise, and the validation set consists of repetitions of a different frozen OU noise stimulus lasting 10 s (only one repetition is shown). (D) Representative validation data and model predictions for each cell type. OU noise input current (top), recorded and predicted voltage traces (middle), and recorded and predicted spike times across all repetitions of the validation stimulus. Stimulus parameters were adjusted for each cell type, see Figure 3—figure supplement 1. (**E–G**) Quantification of model performance in terms of $R^{2}$ on the training subthreshold voltage derivative and on the validation spike-train similarity metric ${M_{d}}^{*}$ . GIF and inactivation GIF (iGIF) models have the same subthreshold performance because the subthreshold components of these models are identical (see Methods). Benchmarks are for models fitted to $N = 18$ serotonin (5-HT), $N = 14$ somatostatin (SOM), and $N = 7$ medial prefrontal cortex (mPFC) neurons. aGIF, augmented GIF.

Our results outlined in Figure 2 show that I_A regulates spike timing in 5-HT neurons because of its nonlinear subthreshold effects. Foreseeing that the presence of this prominent current may limit the accuracy of canonical GIF models—which are not designed to capture nonlinear subthreshold effects—we first augmented the canonical GIF model (aGIF; Figure 3A2) with a simplified Hodgkin-Huxley-type model of the subthreshold voltage-dependent currents we recorded in 5-HT neurons (see Methods). To assess whether incorporating additional biophysical details into the aGIF model might further improve its predictive performance, we turned to the previously described sodium channel-inactivation GIF model (iGIF; Figure 3A3), which extends the GIF model of Mensi et al., 2012 by adding a non-parametric voltage coupling function to the dynamic spike threshold (Mensi et al., 2016; see Methods). Although this GIF model extension was initially conceived specifically to capture the influence of subthreshold sodium channel inactivation on firing threshold (hence, its name), the non-parametric definition of the threshold coupling function gives it the capacity to account for a wide range of other subthreshold biophysical mechanisms which regulate spiking, notably including, but not limited to, I_A. Comparing the performance of the more parsimonious aGIF model to that of the iGIF model enabled us to assess whether accounting for additional mechanisms that regulate spiking beyond I_A might further improve our DRN neuron models.

To establish comparative GIF model benchmarks across cell types, we carried out whole-cell electrophysiological recordings not only from DRN 5-HT and SOM cells but also from canonical deep-layer pyramidal neurons of the medial prefrontal cortex (mPFC). For each recording, we applied two distinct instantiations of noisy in vivo-like inputs (see Supplemental methods, Figure 3—figure supplement 1), one of which was used to determine the model parameters while the other was reserved for post hoc evaluation of the models’ accuracy (i.e. ‘training’ data and ‘validation’ data, respectively; see Figure 3C). Accuracy was assessed by comparing models with recorded data across cell types in terms of Figure 3D: (1) subthreshold voltage changes on training data, R², and; (2) spike timing on validation data, ${M_{d}}^{*}$ (where ${M_{d}}^{*} = 1$ is the best possible performance and ${M_{d}}^{*} = 0$ is the chance level; see Methods).

The canonical GIF model predicted both the subthreshold dynamics and spike timing of mPFC pyramidal neurons with high accuracy (R² = 0.431 ± 0.249; ${M_{d}}^{*}$ = 0.783 ± 0.134; Figure 3E), consistent with previous reports on cortical pyramidal neurons (Mensi et al., 2012; Pozzorini et al., 2015; Mensi et al., 2016; Teeter et al., 2018). While our aGIF model slightly better predicted the voltage of mPFC neurons (R² = 0.544 ± 0.280, $p = 0.028$ , Figure 3E1), this did not translate into more accurate spike predictions ( ${M_{d}}^{*}$ = 0.743 ± 0.180, $p = 0.710$ , Figure 3E2), consistent with the observation that I_A is not a significant conductance recorded from the cell body of mPFC pyramidal neurons (Figure 3—figure supplement 2 and see Dong and White, 2003; Dong et al., 2005). On the basis of spike timing prediction, the canonical GIF model thus offered the most parsimonious account of the behavior of mPFC neurons.

With this point of comparison established, we next quantified the performance of each of our candidate GIF models (GIF, aGIF, and iGIF) in 5-HT neurons. As previously intuited, the canonical GIF model performed rather poorly in 5-HT neurons (Figure 3F), predicting <15% of the variance of the subthreshold voltage (R² = 0.128 ± 0.135) and achieving ${M_{d}}^{*}$ scores less than half of those observed in mPFC neurons ( ${M_{d}}^{*}$ = 0.352 ± 0.118). This indicates that the passive membrane filter and adaptation mechanisms included in the canonical GIF model were insufficient to capture the behavior of 5-HT neurons. By augmenting the GIF model with our experimentally constrained model of I_A, the aGIF model not only better predicted the voltage (R² = 0.301 ± 0.200, p=1.96e-4; Figure 3F1) but also the spike timing ( ${M_{d}}^{*}$ = 0.481 ± 0.148, p=0.001; Figure 3F2) of 5-HT neurons. While the more general iGIF model exhibited a similar improvement in spike timing predictions over the GIF model ( ${M_{d}}^{*}$ = 0.536 ± 0.154, p=5.89e-4), it did not significantly outperform the aGIF model (p=0.644; Figure 3F2), suggesting that accounting for additional biophysical mechanisms that regulate spiking beyond those included in the aGIF model would be unlikely to further improve performance. Repeating this process using data collected closer to physiological temperature yielded the same result (Figure 3—figure supplement 3). Thus, among the models considered, adding I_A to the subthreshold and spiking mechanisms of the GIF model best accounts for the biophysical mechanisms responsible for shaping the responses of 5-HT neurons to in vivo-like inputs.

Turning to the other main cell type of the DRN, we next analyzed the performance of each model in SOM cells (Figure 3G). In these cells, the canonical GIF model produced highly accurate predictions (R² = 0.600 ± 0.238 and ${M_{d}}^{*}$ = 0.818 ± 0.149), consistent with its high performance previously reported for cortical GABAergic neurons (Mensi et al., 2012; Teeter et al., 2018). Nonetheless, the iGIF achieved small but significant performance gains ( ${M_{d}}^{*}$ = 0.892 ± 0.094, p=0.004 vs. aGIF and p=0.003 vs. GIF; Figure 3G2), leading us to select it as our model of SOM neurons.

Multiple adaptation mechanisms in 5-HT neurons

Our model selection approach identified the most salient components required to capture the input-output functions of individual neurons and allowed us to identify functional differences across cell types. 5-HT neurons were distinguishable from SOM and mPFC cells by their long membrane time constants (Figure 4A, Figure 4—figure supplement 1A) and by the presence of conspicuously potent and protracted adaptation mechanisms (Figure 4B–D). Indeed, in addition to evoking a characteristically large and prolonged adaptation current (Figure 4C), action potential firing in 5-HT neurons produced a substantial and long-lasting increase in firing threshold (Figure 4D; but note that this effect is somewhat attenuated near physiological temperature, Figure 4—figure supplement 2). In contrast, SOM neurons most often displayed either negligible or even depolarizing spike-triggered currents (Figure 4B and C) that may underlie the burst firing patterns often observed in this cell type (Figure 1—figure supplement 2). These observations derived from the parameters of GIF models are not only consistent with our experimental characterization (Figures 1—3), but significantly expand it. Thus, 5-HT neurons are characterized by slow membrane dynamics, I_A, and particularly prominent adaptation mechanisms.

Figure 4 with 2 supplements see all

Download asset Open asset

Serotonin (5-HT) neurons are distinguished by slow membrane time constants and potent adaptation.

(A) Using the features from the best performing generalized integrate-and-fire (GIF) model variant for each cell type (legend), passive membrane time constant. (B) Spike adaptation features: potency of after-hyperpolarization potential (AHP)-mediated (spike-triggered current $η$ integral) and AHP-independent (spike-triggered threshold movement $γ$ integral) adaptation. (**C, D**) Comparison of model filters. Presented as median (lines) and interquartile range (bands). Note the long-lasting adaptation currents (C; positive values indicate hyperpolarizing current) and threshold movements (D) of 5-HT neurons. Parameters are from models fitted to $N = 18$ 5-HT, $N = 14$ somatostatin (SOM), and $N = 7$ medial prefrontal cortex (mPFC) neurons.

Preferential sensitivity of 5-HT neuron population to the onset of sudden inputs

The development and validation of accurate single-cell models allowed us to identify the population-level computations operating in the DRN. We took advantage of the one-to-one correspondence between our GIF models and real neurons to construct synthetic populations with realistic neuron-to-neuron heterogeneity by sampling from banks of single-cell models (Figure 5A). In response to step increases of synaptic-like inputs delivered to the entire population (Figure 5B left), the population firing rates (in Hz/neuron; Figure 5B right) of 5-HT, SOM, and mPFC neurons (Figure 5C) transiently increased before relaxing to a significantly lower stationary level. Strong inputs did not produce oscillations in the population firing rates, likely because of population heterogeneity (Figure 5—figure supplement 3; Naud and Gerstner, 2012; Mejias and Longtin, 2012; Mejias et al., 2014; Tripathy et al., 2013). The transient and stationary parts of the population input-output functions were approximately rectified linear functions (Figure 5—figure supplement 3) which we summarized and plotted as the time-varying slope (i.e. gain; Figure 5D). While the gain of the transient response was greater than that of the stationary response in all three cell types, the ratio of transient to stationary gain was substantially higher in 5-HT neurons (Figure 5E; ratio of 3.42 ± 0.07 vs. 1.89 ± 0.04 in SOM and 1.50 ± 0.03 in mPFC; p<0.001 in each case; but note that the gain ratio in 5-HT neurons falls to 2.13 ± 0.03 near physiological temperature [Figure 5—figure supplement 5], consistent with a smaller spike-triggered threshold movement [Figure 4—figure supplement 2]). This marked response of 5-HT cells occurred quickly, in the first 100 ms after the onset of the step. Thus, despite 5-HT neurons being characterized by slow membrane time constants, their population activity provided a remarkably strong encoding of the onset of step synaptic inputs.

Figure 5 with 5 supplements see all

Download asset Open asset

Adaptation mechanisms cause a higher gain of the transient vs stationary population response.

(A) Generation of heterogeneous population models from experimentally constrained single neuron models. (B) Schematic of population simulations. Spikes from individual neuron models in the simulated population are added together to produce a population firing rate. (C) Population responses to input step. From top in each column: stimulus (gray); sample voltage trace; spike raster of first 20 neurons; mean population firing rate across 20 independent simulations. (D) Schematic for quantifying the time-varying population input-output function for both the transient and the stationary components of the response. An input-output function is calculated for the population response at each time point after the input step. The slope of each input-output function (gain) is then plotted as a function of time since the step onset. The ratio of the maximum gain to the minimum gain is a measure of the relative amount of population adaptation. (E) Time-resolved gain of step input responses across cell types following the approach shown in D. (F) Time-resolved gain of serotonin (5-HT) populations with the adaptation parameters of somatostatin (SOM) neurons (blue) and of SOM populations with adaptation parameters of 5-HT neurons (orange). Data are presented as mean ± SD in E1 and F1. mPFC, medial prefrontal cortex.

We next considered the underlying mechanisms giving rise to the distinctive time-dependent gain of 5-HT neurons. We found that the characteristically strong spike-triggered adaptation of 5-HT neurons (spike-triggered hyperpolarizing adaptation current and threshold movement shown in Figure 4) contributed to the observed relaxation of the population response to a lower stationary level: grafting the weak adaptation from SOM neuron models onto 5-HT models dramatically reduced the ratio of transient to stationary gain, and vice-versa (Figure 5F). These findings are consistent with previous models in other cell types showing that spike-triggered adaptation reduces the sensitivity of neural populations to input changes over long timescales (Ermentrout, 1998; Benda and Herz, 2003; Naud and Gerstner, 2012). Therefore the preferential sensitivity of 5-HT neuron populations to sudden changes in synaptic inputs is a natural consequence of strong adaptation at the single neuron level.

Feedforward inhibition and I_A control 5-HT output gain of the DRN

Apart from the strong adaptation mechanisms of 5-HT neurons, two other mechanisms have the potential to dynamically modulate the 5-HT output from the DRN: I_A in 5-HT neurons and the feedforward inhibition (FFI) enacted by local DRN interneurons (Zhou et al., 2017; Geddes et al., 2016). To examine the contributions of these two mechanisms, we first connected our existing SOM population models to 5-HT population models using experimentally constrained GABA_A receptor-mediated synaptic conductances (see Methods and Figure 5—figure supplement 4).

To dissect the contribution of I_A in shaping population responses in this connected DRN network, we applied the same inputs to both 5-HT and SOM neuron populations and examined 5-HT neuron population dynamics (as in Figure 5) while varying the maximal conductance of I_A (in 5-HT neurons). The gain of the transient component of the 5-HT response increased markedly when the conductance of I_A was set to zero (Figure 6A), while increasing the potency of I_A substantially dampened and broadened the population response to fast inputs, reminiscent of I_A’s modulation of spike timing jitter observed in our toy model (Figure 2I–K). These simulations thus show that I_A substantially regulates the gain of the transient component of DRN 5-HT output evoked by sustained inputs, with negligible effects on the gain of the slower stationary component.

Figure 6

Download asset Open asset

Effect of *I_A* density, feedforward inhibition, and heterogeneity of somatostatin (SOM) neurons on the serotonin (5-HT) neuron population response.

Network input is the same set of step stimuli as in Figure 5D–F. (A) Increasing I_A reduces adaptation by selectively suppressing the early part of the response to sudden inputs, and vice-versa. (B) Gain curves with normal feedforward inhibition (blue), with reduced input strength onto the inhibitory population (green), or without inhibition (red). Reduced input strength onto the inhibitory population (green) simulates the effect of endocannabinoid input (Geddes et al., 2016). (**A2,B2**) Ratio of peak to steady-state gain. Data are presented as mean ± SD. (C) Population firing rates of SOM and 5-HT neurons in a network in which both populations are heterogenous. (D) Population firing rates of SOM and 5-HT neurons in a network in which all SOM neurons are identical. Effects of homogeneous (cyan) or heterogeneous (green) SOM populations on the population input-output functions for the transient (E) and stationary (F) components of the response (see square and circle markers in C and D). Note that the input-output function of the heterogenous SOM population is approximately linear, whereas that of the homogenous population is not (**E1, F1**). Relative to the input-output functions of a 5-HT population receiving no feed-forward inhibition, the effect of the heterogenous SOM population is divisive, but the effect of the homogenous SOM population on the transient part of the 5-HT population input-output function includes a strong subtractive component (E2).

Previous work has shown that glutamatergic excitatory inputs from the PFC make strong mono-synaptic contacts onto both DRN 5-HT and GABAergic neurons, triggering a classic FFI. Intriguingly, the PFC axonal inputs onto these two cellular elements of the DRN are functionally distinct in as much as the PFC synapses onto GABAergic neurons are far more sensitive to endocannabinoid neuromodulation than those onto 5-HT neurons (Geddes et al., 2016). The computational role of this differential sensitivity to neuromodulation is currently unknown. We began by determining the role of the DRN FFI per se by comparing the responses of 5-HT neuron population dynamics with or without SOM cells (Figure 6B). Including FFI onto 5-HT neurons substantially dampened the overall response of the 5-HT population to synaptic inputs, while still sustaining the preferential encoding of the early phase of sudden inputs (Figure 6B2). While introducing FFI did decrease the gain ratio, this decrease was quantitatively smaller than the differences between 5-HT neurons and other cell types shown in Figure 5E and the effect of changing I_A shown in Figure 6A, Figure 6B2. We next directly simulated the effects of endocannabinoid modulation of excitatory input to the DRN observed experimentally (Geddes et al., 2016) by weakening the strength of the input to SOM neuron populations by 30% while leaving that to 5-HT neurons intact. By favoring the direct monosynaptic excitation of 5-HT neurons by preferentially diminishing the glutamatergic drive of SOM neurons, this neuromodulation led to an increase in the overall gain of the DRN that was unexpectedly apparent across the entire duration of the response to step inputs (i.e. no change in the gain ratio, Figure 6B2). Thus, the target-specific endocannabinoid-mediated modulation of PFC excitatory drive in DRN exerts a normalizing role by increasing the overall gain of 5-HT output evoked by synaptic inputs without altering its preferential encoding of changes in input, which is emerging as a cardinal feature of DRN network dynamics.

Our electrophysiological recordings showed that excitability heterogeneity is a salient feature of the SOM DRN neuron population. Our modeling approach allows us to specifically examine the role of this cellular heterogeneity in shaping the output of the DRN by comparing our DRN model (Figure 6C) to an alternative homogenized version in which the parameters of SOM neurons were set to fixed values (Figure 6D). Thus, while FFI with an experimentally determined degree of heterogeneity mainly imposed a reduction of the slope of the input-output function (i.e. divisive inhibition), homogeneous FFI mainly shifted the input-output function of the transient component of the population response to the right (i.e. subtractive inhibition; Figure 6E). This subtractive feature can be traced back to a strong non-linearity in the input-output functions of homogenized SOM neuron populations (compare Figure 6E1 and Figure 6F1). In the case of the stationary component, both heterogeneous and homogenized DRN models implemented divisive inhibition (Figure 6F). Therefore, we conclude that heterogeneity among GABAergic neurons implements divisive inhibition.

5-HT neurons linearly encode the temporal derivative of inputs to the DRN

Adaptation plays a critical role in implementing temporal derivative encoding in sensory systems (Lundstrom et al., 2008; Pozzorini et al., 2013) but has not been ascribed a similar role in neuromodulatory systems such as the DRN. To determine whether the DRN also supports this computation, we parameterized the rate of change of DRN inputs by applying ramp stimuli with variable slopes (i.e. derivatives; Figure 7A and B). Remarkably, the peak 5-HT neuron population firing rate linearly reported the slope of the ramps, an effect which was enhanced by FFI (Figure 7C). We further found that this linearity was conditional on the presence of slightly depolarizing background input (≥20 pA, Figure 7D). Simulations using aGIF models fitted to data collected near physiological temperature yielded similar results; Figure 7—figure supplement 1. The potent adaptation mechanisms of 5-HT neurons play a key role in mediating this linear encoding of input derivative, since reducing the strength of adaptation reduces linearity across a wide range of input baselines (Figure 7E and F). Together, these observations suggest that the DRN signals to its brain-wide target a mixture of the intensity and temporal-derivative of its excitatory inputs, and that the derivative-encoding component dominates when the input is increasing rapidly (Figure 7—figure supplement 2).

Figure 7 with 2 supplements see all

Download asset Open asset

Dorsal raphe nucleus (DRN) serotonin (5-HT) neuron population output conditionally encodes the temporal derivative of its input.

(A) Design of simulations. A ramp stimulus with an adjustable baseline and slope (derivative) is applied to the same network models as in Figures 5 and 6, and the peak firing rate (FR) of the 5-HT neuron population is extracted. (B) Representative simulated input (top), somatostatin (SOM) neuron population activity (middle), and 5-HT neuron population output (bottom). (C) With a baseline input of 40 pA, peak 5-HT neuron population output is approximately linearly related to the derivative of the ramp input, and feed-forward inhibition by SOM neurons enhances this feature. (D) Peak FRs of 5-HT neuron populations depend on interacting effects of input baseline and slope. Panel C shows normalized data from the 40 pA row in blue. (E) 5-HT neuron adaptation and I_A dominate the DRN input-output function under different input regimes. Effect of reducing 5-HT neuron adaptation (following the approach from Figure 5F) is the most pronounced for higher levels of background input and more slowly changing inputs (red), while the effect of removing I_A (following the approach of Figure 6A) is the most pronounced for low background input and fast changing inputs (orange, blue). (F) Effect of reducing adaptation in 5-HT neuron models visualized at a 40 pA baseline. Note that 5-HT output no longer linearly encodes *dI/dt* when adaptation is reduced. (G) Effect of removing I_A from 5-HT neuron models visualized at a 0 pA baseline. Note that 5-HT output approximately linearly encodes *dI/dt* when I_A is removed.

The extent to which the output of the DRN signals the temporal derivative of its input is likely to be limited by several factors, notably: the long membrane time constants of 5-HT neurons (Table 1, Figure 4—figure supplement 1A), which cause rapidly fluctuating to be filtered out; the fact that firing rates cannot be less than zero, limiting the dynamic range available to encode negative input derivatives; the presence of I_A , which filters out inputs with a high temporal derivative (Figure 2); and the level of background input (Figure 7D). Because I_A can be partly inactivated by depolarizing background input, the effects of background input and I_A on the derivative-encoding properties of the DRN are expected to interact. Consistent with this idea, removing I_A from 5-HT neurons in our DRN network models extended the range of background input where the peak 5-HT neuron population firing rate is an approximately linear function of the slope of a ramp stimulus (Figure 7E and G). In summary, we found that the presence of strong spike-frequency adaptation in 5-HT neurons causes the DRN to signal the rate of change of its input to its brainwide targets, but that this core computation is progressively suppressed when a state of hyperpolarization engages I_A .

Discussion

Here, we sought to characterize the computational properties of the DRN using a bottom-up approach grounded in experimentally constrained models of the two most abundant cell types in this region: 5-HT and SOM GABA neurons. Consistent with, and extending, previous work, we found that 5-HT neurons were relatively homogeneous and characterized by potent spike-frequency adaptation (Figure 1) and by the presence of a strong A-type potassium current (Figure 2), while SOM neurons displayed a considerably more heterogeneous excitability profile (Figure 1 and Figure 1—figure supplement 2). Extensions to classical GIF models (Mensi et al., 2012; Pozzorini et al., 2015) to capture the non-linear subthreshold effects of I_A observed in 5-HT neurons were required to adequately capture the spiking response of 5-HT neurons to naturalistic stimuli (Figure 3). This work introduces a new approach to capturing such non-linear subthreshold effects in the form of the aGIF model, which augments the GIF model of Mensi et al., 2012 with experimentally constrained Hodgkin-Huxley style currents, improving model interpretability without compromising predictive performance. Inspecting the parameters of the best performing GIF models revealed that the substantial spike-frequency adaptation observed in 5-HT neurons is not fully explained by their distinctively large AHPs and is partly mediated by a previously undescribed dynamic spike threshold (Figure 4). This model-based approach allowed us to probe causal relationships between specific excitability features and population computations. Thus, we found that the prominent adaptation mechanisms in 5-HT neurons regulated DRN population responses to synaptic inputs (Figure 5), that I_A suppressed the response to sudden inputs, and that heterogenous FFI had a divisive rather than subtractive effect on DRN output (Figure 6). By further exploring DRN population dynamics, our simulations demonstrated that 5-HT neurons linearly reported a mixture of the intensity and temporal-derivative of their synaptic inputs (Figure 7), and that the temporal-derivative dominates DRN output when the input is increasing rapidly (Figure 7—figure supplement 2). In summary, this work points to a new computational role for the DRN in encoding the derivative of its inputs and identifies specific cellular and network mechanisms that give rise to this computation and modulate its expression. These results raise important questions about how the selective responses of the DRN to changing synaptic inputs might support its role in guiding animal’s behavior in dynamic environments.

Need for a hybrid biophysical-simplified methodology

The computational and statistical modeling methodology presented here was designed to bridge the gap between specific biophysical mechanisms and network-level computation. Closing this gap has also been the target of complex biophysical simulations, motivated by the hope to create tools for testing disease-related treatments and for untangling the computations performed by large neural networks (Markram, 2006; Billeh et al., 2020). Preserving the accuracy and identifiability of simpler approaches (Gerstner and Naud, 2009; Mensi et al., 2012; Pozzorini et al., 2013; Teeter et al., 2018), the ‘augmented GIF’ model developed here explicitly incorporates the most important biophysical features of 5-HT neurons, allowing us to probe their contributions to network-level computation by altering or removing the corresponding model components during network simulations. While the aGIF framework was developed here to capture the effects of inactivating subthreshold potassium currents in 5-HT neurons, it lends itself equally well to capturing the effects of other subthreshold voltage-gated currents. We note that, as in other methods based on linear regression of nonlinear ion channel dynamics (Huys et al., 2006; Huys and Paninski, 2009), adequate experimental estimates of the voltage-dependent gating features of the conductance at play must be available to be inserted in the aGIF model. Altogether, this expanded modeling framework adds to a toolset of computational approaches for interrogating the role of particular microcircuit motifs (e.g. FFI) or excitability features (e.g. spike-triggered adaptation) in shaping network computations, while lending itself to more elaborate inference methodologies (Gonçalves et al., 2019).

Could the dynamical features identified here have been captured by a simpler modeling framework? Two closely related approaches that we have not considered here are linear-nonlinear (LNL) and generalized linear models (GLMs), which are trained using only the spike output and external input to each cell and do not consider the subthreshold voltage (Pillow and Simoncelli, 2006; Pillow et al., 2008). Despite the fact that the GLM approach was not possible here given the very low firing rates of 5-HT neurons and the large number of action potentials required for accurate characterization in the absence of information about the subthreshold voltage, it is worth asking whether GLMs could in principle capture the network-level properties of 5-HT signaling. For instance, the role of spike-triggered adaptation in conveying preferential sensitivity to suddenly changing inputs arises in GLMs (Naud and Gerstner, 2012), but the state-dependence of the input derivative sensitivity identified in 5-HT neurons (Figure 7) could not have been captured by a GLM implementation. In summary, the GIF framework provides a more solid foundation for network modeling than LNL- or GLM-based approaches for cell types with very low firing rates or highly state-dependent output.

Does the aGIF modeling approach represent an unnecessary complication of the GIF model framework or, conversely, an oversimplification of detailed Hodgkin-Huxley models? GIF models that do not explicitly account for the effects of specific ionic conductances produce highly accurate spiketrain predictions in many cell types (Gerstner and Naud, 2009; Mensi et al., 2012; Pozzorini et al., 2013; Teeter et al., 2018); indeed, even in 5-HT neurons, the iGIF model predicts the timing of spikes with an accuracy equal to that of our aGIF model. For questions where the biophysical mechanisms that regulate spiking are not of primary interest and for systems where simpler LNL or GLM models are not able to predict the timing of spikes accurately (e.g. due to low firing rates as discussed above), non-augmented GIF models remain suitable tools. In our case, it would not have been possible to probe the effect of I_A on the network-level processing features of the DRN without the aGIF model.

Network-level role of I_A current

Previous modeling work has implicated I_A in controlling the sensitivity of the stationary response to sustained inputs (Connor and Stevens, 1971; Connor et al., 1977; Tuckwell and Penington, 2014; Drion et al., 2015). These studies contrast with our findings which implicate this current in the control of the transient component but show almost no effect on the stationary component of the response. This discrepancy can be explained by noting that the AHPs of 5-HT neurons (and thus of our computational model) do not reach the hyperpolarized potentials required to free I_A from inactivation (Figure 1 and Figure 1—figure supplement 1), in contrast to the model of Connor et al., 1977. As a result, I_A remains mostly inactivated during sustained inputs, and the stationary response is mostly regulated by the interplay between spike-triggered adaptation and the strength of the input. Other factors such as a shift in the activation and/or inactivation curves (e.g. by neuromodulators) are expected to influence how I_A controls the transient and stationary components of the response. Finally, it is interesting to note that I_A is also highly expressed in the dendrites of cortical neurons, where it may have an analogous function (Hoffman et al., 1997; Harnett et al., 2013; Ujfalussy et al., 2018; Payeur et al., 2019). Our results hint at a possible general role of I_A in suppressing transient responses to sustained inputs in the midbrain, cortex, and other systems.

5-HT neuron heterogeneity

5-HT neurons are not all alike in every respect: recent experimental work has uncovered molecular, electrophysiological (Calizo et al., 2011), developmental, and anatomical (Commons, 2015; Ren et al., 2018) differences among 5-HT neurons across raphe nuclei and within the DRN (reviewed in Okaty et al., 2019). Most relevant to our work are previously reported quantitative differences in the excitability of serotonin neurons located in the dorsomedial DRN, ventromedial DRN, and median raphe nucleus (Calizo et al., 2011). These observations suggest that the predictions made by our model, which was fitted primarily to serotonin neurons from the ventromedial DRN, may agree qualitatively but not quantitatively with the behavior of 5-HT neuron ensembles in these areas. While there is not yet any evidence that serotonin neurons in different parts of the serotonin system perform qualitatively different computational operations, this remains an intriguing possibility for future work.

Heterogeneous properties of SOM neurons ensure divisive inhibition

How the heterogeneity of excitability influences the response properties of neuronal populations depends on a number of factors. Specifically, we and others Mejias and Longtin, 2014 have argued that heterogeneity of feedback inhibition (and of principal cells) implements a divisive effect on the stationary part of the population input-output function. For FFI, a divisive effect on the gain of stationary input-output functions is expected in naturalistic conditions (Mejias and Longtin, 2014). The findings outlined here further support these theoretical results by showing that the heterogeneous FFI remains divisive on the transient part of the response. Divisive inhibition has been proposed to be essential to counteract strong excitation so as to maintain activity within an adequate dynamic range (Chance and Abbott, 2000; Ferguson and Cardin, 2020), and it is expected that brain circuits will harness cellular and circuit-level mechanisms to tune their sensitivity to relevant inputs while maintaining overall stability. This point is germane to 5-HT neurons given their position at the confluence of many excitatory input streams (Weissbourd et al., 2014; Pollak Dorocic et al., 2014; Ogawa et al., 2014; Zhou et al., 2017; Ren et al., 2018; Geddes et al., 2016). Thus, while the exact behavioral function of the 5-HT system is still unclear, uncovering important components of its gain control mechanisms might provide useful hints about how it integrates its multifold inputs.

Neuromodulation of neuromodulation

Neuromodulators can dynamically reconfigure information processing in neural circuits that are otherwise anatomically fixed (Marder, 2012; Tsuda et al., 2021). While 5-HT is considered to be a neuromodulator, the DRN network is itself under neuromodulatory influence, both from distal (e.g. locus coeruleus or ventral tegmental area) or local (e.g. endocannabinoids, 5-HT itself) sources (Baraban and Aghajanian, 1981; Aman et al., 2007; Weissbourd et al., 2014; Geddes et al., 2016; Lynn et al., 2022). Whereas previous work has outlined defined cellular metrics that are modulated by specific receptor subtypes (e.g. changes in release probability or direct membrane depolarization/hyperpolarization), the consequences of these neuromodulatory influences on higher-order network computation are only superficially understood. Here, we showcase two broad neuromodulatory mechanisms that enact different effects on population coding. Through simulations, we show that reducing the magnitude of I_A (which could be caused, for instance, in vivo by noradrenergic input the DRN [Aghajanian, 1985]) enhances the sensitivity of the raphe response to the onset of step inputs while leaving the stationary firing rate unchanged. In contrast, the cannabinoid-mediated preferential reduction of FFI onto 5-HT neurons (caused by the tonic activation of DRN endocannabinoid receptors, as expected to occur, for instance, during marijuana recreational or therapeutic use [Geddes et al., 2016]) rather causes a general reduction in the output gain of the DRN. Together with our simulations probing the temporal derivative-encoding properties of this region, these observations point to a conceptual model in which the output of the DRN represents a mixture of the intensity and temporal derivative of its input where I_A controls the relative balance of the two components, and FFI regulates the overall intensity of the output, and where these functions can be rapidly and independently tuned by neuromodulatory control.

Our heuristic model of the DRN helps to illustrate the unexpectedly multifaceted nature of the computations performed by this evolutionarily ancient region, but, like most heuristics, it remains an oversimplification. Some of the qualitative features of DRN processing emerging from our simulations are not explained by our “input intensity plus temporal derivative” heuristic (e.g. the ability of FFI to modulate the temporal derivative-encoding properties of the DRN, or the attenuation of these same coding properties by hyperpolarization [Figure 7]), presenting further opportunities to better understand the influence of neuromodulation on network computation in this region.

Role of derivative encoding in reinforcement learning

The role of 5-HT signaling in modulating behavior is increasingly conceptualized through the lens of reinforcement learning (RL) theory. Indeed, 5-HT output has been proposed to loosely encode or modulate every component of classical RL (Sutton and Barto, 2018, Dayan and Huys, 2009), including a reward signal (Li et al., 2016), state value (Cohen et al., 2015; Luo et al., 2016), bias in state-action value (Miyazaki et al., 2018), temporal discounting factor (Doya, 2002; Schweighofer et al., 2008), prediction error (Daw et al., 2002, but see Boureau and Dayan, 2011), and learning rate (Matias et al., 2017; Iigaya et al., 2018; Grossman et al., 2022), with varying degrees of experimental support. Might the derivative-like computation described here have a place in an RL-based conception of DRN function? For now, it is only possible to speculate. Existing RL models of DRN function bin time in increments of tens of seconds, obscuring the faster adaptation dynamics that are the subject of our work. How and whether the sub-second fluctuations in DRN 5-HT neuron activity that are consistently observed in reward learning experiments (Ranade and Mainen, 2009; Cohen et al., 2015; Li et al., 2016; Zhong et al., 2017; Grossman et al., 2022) should be incorporated into RL models remains unclear. Our results suggest that RL operations that can be seen as computing a temporal derivative are candidates for an RL-based account of DRN function.

If the electrophysiological features of individual 5-HT neurons directly participate in shaping the computations enacted by the DRN, the same is likely true for other neuromodulatory systems, and this work may offer overall guiding principles. For instance, dopamine neurons, well known for their reward prediction error-like coding properties (Schultz et al., 1997), bear some electrophysiological features in common with DRN 5-HT neurons, with both cell types exhibiting strong adaptation and a prominent A-type potassium current (Grace and Onn, 1989, Khaliq and Bean, 2008). Dopamine neurons have been proposed to encode reward prediction errors partly by approximating a mixture of a value signal and its temporal derivative (Kim et al., 2020), hinting at a possible role for adaptation in implementing one of the central computations of RL.

If the derivative-like operation identified here does not directly contribute to computing one of the key components of RL, what might its role in the DRN be? One possibility is that strong spike-triggered adaptation may optimize the efficiency of neural coding by filtering out temporally redundant information, a phenomenon referred to as predictive coding and that is ubiquitous in sensory systems (Brenner et al., 2000; Barlow, 2001; Ulanovsky et al., 2003; Kohn, 2007). As the search for a unified interpretation of DRN 5-HT activity continues, our results provide a new perspective on the fast component of 5-HT neuron dynamics: fluctuations in 5-HT neuron’s activity do not solely encode the intensity of their input, but rather how quickly their inputs are changing over time.

Materials and methods

Experimental methods

Animals

Experiments were performed on male and female C57/Bl6 mice aged 4–8 weeks. Slc6a4-cre::Rosa-TdTomato (SERT-Cre) and Sst-cre::Rosa-TdTomato transgenic lines were used to fluorescently label DRN 5-HT and somatostatin (SOM) GABA neurons, respectively. Animals were group-housed and kept on a 12:12 hr light/dark cycle with access to food and water ad libitum. All experiments were carried out in accordance with procedures approved by the University of Ottawa Animal Care and Veterinary Services (protocol numbers CMM-164, CMM-176, CMM-1711, CMM-1743, and CMM-2737).

Slice preparation

Request a detailed protocol

Animals were deeply anesthetized using isofluorane (Baxter Corporation) before being euthanized by decapitation. The brain was quickly removed from the skull and submerged into ice-cold dissection buffer containing the following: 119.0 mM choline chloride, 2.5 mM KCl, 4.3 mM MgSO₄, 1.0 mM CaCl₂, 1.0 mM NaH₂PO₄, 1.3 mM sodium ascorbate, 11.0 mM glucose, 26.2 mM NaHCO₃; saturated with 95% O₂/5% CO₂. A Leica VT1000S vibratome was used to cut 300-µm coronal sections of midbrain containing the DRN or of the cortex containing the mPFC in the same ice-cold choline dissection buffer. After cutting, slices were placed in a recovery chamber filled with artificial cerebrospinal fluid containing the following: 119.0 mM NaCl, 2.5 mM KCl, 1.3 mM MgSO₄, 2.5 mM CaCl₂, 1.0 mM NaH₂PO₄, 11.0 mM glucose, 26.2 mM NaHCO₃; ∼298 mOsm, maintained at 37°C, and continuously bubbled with 95% O₂/5% CO₂. The recovery chamber was allowed to equilibrate to room temperature for 1 h before beginning experiments.

In vitro whole-cell electrophysiological recording

Request a detailed protocol

Neurons were visualized using an upright microscope (Olympus BX51WI) equipped with differential interference contrast and a ×40, 0.8 NA water-immersion objective. Whole-cell recordings were obtained from fluorescently labeled DRN 5-HT and SOM neurons and unlabeled mPFC L5 pyramidal neurons using glass electrodes (Sutter Instruments; tip resistance 4–6 MOhm). For most experiments, the following potassium gluconate-based internal solution was used: 135 mM potassium gluconate, 6.98 mM KCl, 10 mM HEPES, 4 mM Mg ATP, 0.40 mM GTP, 10 mM Na phosphocreatine; adjusted to pH 7.25 with KOH, 280–290 mOsm. A subset of experiments (GABA synaptic physiology) were carried out using a cesium-based internal solution (120 mM CsMeSO₃, 10 mM EGTA, 5 mM TEA Cl, 1 mM CaCl₂, 10 mM Na HEPES, 4 mM Mg ATP, 2 mM GTP, 2 mM QX-314, and 10 mM Na phosphocreatine; adjusted to pH 7.25 with CsOH, 280–290 mOsm) and in the presence of bath-applied 100 µM (2 R)-amino-5-phosphonovaleric acid (APV) and 5 µM 2,3-dioxo-6-nitro-1,2,3,4-tetrahydrobenzo[f]quinoxaline-7-sulfonamide (NBQX). For voltage clamp experiments, whole-cell capacitance compensation was applied manually following break-in, and leak current subtraction was performed post hoc using membrane leak conductance estimated based on a –5 mV pulse at the start of each sweep. Experiments were carried out at room temperature except where noted. For current clamp experiments used to fit GIF models, access resistance was compensated using an active electrode compensation method (Pozzorini et al., 2015). For voltage clamp experiments used to characterize I_A in 5-HT neurons at room temperature, recordings had $R_{a} =$ 14.7 ± 6.2 MOhm (mean ± SD; half of recordings between 12.8 MOhm and 21.6 MOhm) after applying an access resistance cutoff of 30 MOhm (a more stringent cutoff of 20 MOhm yielded statistically indistinguishable estimates of I_A maximal conductance and kinetic parameters; compare Figure 1D and Figure 1—figure supplement 5). For voltage clamp experiments used to characterize whole-cell currents in SOM neurons, recordings had $R_{a} =$ 14.3 ± 7.0 MOhm (mean ± SD; half of recordings between 9.8 MOhm and 15.5 MOhm) after applying a similar cutoff of 30 MOhm. For synaptic electrophysiology experiments, recordings had $R_{a} =$ 5.7 ± 0.5 MOhm (mean ± SD; range 5.0 MOhm–6.1 MOhm) after applying a cutoff of 10 MOhm. Recordings were collected with an Axon MultiClamp 700B amplifier, and the analog signals were filtered at 2 kHz and digitized at 10 kHz using an Axon Digidata 1550 digitizer.

Models

GIF and related models

Request a detailed protocol

The GIF and Na-inactivation GIF (iGIF) models have been described previously in detail (Mensi et al., 2012; Pozzorini et al., 2015; Mensi et al., 2016). Briefly, the GIF and iGIF are composed of a subthreshold component which integrates input currents into voltage and a stochastic spiking rule which transforms subthreshold voltage into a series of spikes. The subthreshold dynamics of the GIF and iGIF are given by

C \frac{d V}{d t} = - g_{l} (V (t) - E_{l}) - \sum_{{\hat{t}}_{i} < t} η (t - {\hat{t}}_{i}) + I_{i n j} (t)

where ${{\hat{t}}_{i}}$ is the set of spike times and $η (t) = \sum_{j} w_{j} \exp [- t / τ_{j}^{(η)}]$ is the spike-triggered adaptation current. Here the $w_{j}$ are coefficients estimated from the data and the $τ_{j}^{(η)}$ are fixed hyperparameters; see Appendix for details. The GIF emits spikes according to an inhomogeneous Poisson process with intensity $λ (t)$ , given by

λ (t) = λ_{0} e x p [\frac{V (t) - V_{T}^{*} - \sum_{{\hat{t}}_{i} < t} γ (t - {\hat{t}}_{i})}{Δ V}]

where $V_{T}^{*}$ is the stationary threshold, $γ (t) = \sum_{j}^{} β_{j}^{(γ)} \exp [- t / τ_{j}^{(γ)}]$ is the spike-triggered threshold movement (where the $β_{j}^{(γ)}$ are coefficients estimated from the data and the $τ_{j}^{(γ)}$ are fixed; see Appendix), $Δ V$ is the threshold sharpness (mV; larger values increase the stochasticity of spiking), and $λ_{0} = 1$ Hz is a constant such that $λ (t)$ is in units of Hz. In the iGIF, an additional variable $θ (t)$ is added to the numerator of the exponentiated term in Equation 2 to account for voltage-dependent changes in threshold:

λ (t) = λ_{0} e x p [\frac{V (t) - V_{T}^{*} - \sum_{{\hat{t}}_{i} < t} γ (t - {\hat{t}}_{i}) + θ (t)}{Δ V}]

\frac{d θ}{d t} = \frac{θ_{\infty} (V) - θ}{τ^{(θ)}} .

The equilibrium voltage-dependent change in spike threshold $θ_{\infty} (V) = \sum_{j = 1}^{N_{s t e p}} β_{j}^{(θ)} rect [V; A_{j}, A_{j + 1}]$ is a piecewise constant function of voltage where each ${β_{j}}^{(θ)}$ defines the value of $θ_{\infty} (V)$ over the voltage range ${[A}_{j}, A_{j + 1})$ and $N_{s t e p} = 5$ . The locations of the steps in the piecewise constant function $A_{j}$ are selected based on the data. (See Mensi et al., 2016 for details on the iGIF model.) Our aGIF model is identical to the GIF model except that two Hodgkin-Huxley currents which together capture the voltage-gated potassium currents found in 5-HT neurons (see ‘Potassium current’, below) are added to the subthreshold dynamics given in Equation 1, yielding

C \frac{d V}{d t} = - g_{l} (V (t) - E_{l}) - I_{A} (t) - I_{K} (t) - \sum_{{\hat{t}}_{i} < t} η (t - {\hat{t}}_{i}) + I_{i n j} (t)

as the definition of the subthreshold dynamics of the aGIF model.

The procedures for fitting the GIF and iGIF models to electrophysiological data have also been described previously in detail (Mensi et al., 2012; Pozzorini et al., 2015; Mensi et al., 2016). Briefly, parameter estimation for both models occurs in two stages: first, the subthreshold parameters are estimated by regression, and second, the threshold parameters are estimated by maximizing the likelihood of the observed spiketrain as a function of the threshold parameters. The fitting procedure for the aGIF is very similar to that of the GIF, with adjustments to the subthreshold fitting procedure to accommodate the extra terms in Equation 3 (see Appendix for details). Neurons with non-stationary firing statistics (Pearson correlation between number of spikes and validation sweep number above 0.9) or highly variable spike timing (intrinsic reliability <0.1) were automatically excluded from our analysis. Exclusion criteria were fixed before comparing candidate models.

LIF neuron with an inactivating potassium current

Request a detailed protocol

Our toy model of a neuron with an inactivating potassium current is based on an LIF augmented with $I_{A} (t)$ (see ‘Potassium current’ below):

C \frac{d V}{d t} = - g_{l} (V (t) - E_{l}) - I_{A} (t) + I_{i n j} (t),

where $g_{l}$ and $E_{l}$ are the leak conductance and reversal, respectively, and $I_{i n j} (t)$ is the external input to the model. To reduce the number of free parameters, the model we used is non-dimensionalized with respect to the membrane time constant $τ_{m e m} = C / g_{l}$ and leak conductance $g_{l}$ , yielding

\frac{d V}{d t} = E_{l} - V (t) - {\bar{g}}_{A}^{'} m_{\infty} (t) h (t) (V (t) - E_{K}) + V_{i n j} (t)

where $t$ is in units of the membrane time constant, ${\bar{g}}_{A}^{'} = {\bar{g}}_{A} / g_{l}$ is the effective maximum conductance associated with $I_{A}$ , and $V_{i n j} (t) = I_{i n j} (t) / g_{l}$ is the effective external input. The gating variables $m_{\infty}$ and $h$ are described below in ‘Potassium current’.

Potassium current

Request a detailed protocol

The voltage-gated potassium currents in 5-HT neurons were modeled in terms of an inactivating current and a non-inactivating current, we refer to as I_A and I_K , respectively. These were defined as follows

\begin{array}{ll} I_{A} = {\bar{g}}_{A} m_{\infty} (V) h (t) (V (t) - E_{K}) \\ I_{K} = {\bar{g}}_{K} n_{\infty} (V) (V (t) - E_{K}) \end{array}

where $\bar{g}$ is the maximal conductance; $m$ and $h$ are the activation and inactivation gates of $I_{A}$ , respectively; $n$ is the activation gate of $I_{K}$ ; and E_K = –101 mV is the reversal potential of potassium in our recording conditions. Note that although this value is not physiological, the effect of varying this parameter is very similar to the effect of varying ${\bar{g}}_{A}$ , as we have done in the result section. For simulations involving models fitted to data collected at 29–30°C, $E_{K} = - 89.1$ mV was used. The equilibrium state of each gate $x \in {m, h, n}$ is a sigmoid function of voltage

x_{\infty} (V) = \frac{A_{x}}{1 + e^{- k_{x} (V - {V_{x}}^{*})}}

where ${V_{x^{}}}^{*}$ is the half-activation voltage (mV), $k_{x}$ is the slope (mV^{– 1}), and $A_{x}$ is a scaling factor.

To keep the number of parameters in our current model to a minimum, we assumed that the $m$ and $n$ gates have instantaneous kinetics (allowing their corresponding equilibrium gating functions $m_{\infty}$ and $n_{\infty}$ to be used directly in Equation 4), and that the $h$ gate inactivates and de-inactivates with a single time constant $τ_{h}$ (ms) that does not depend on voltage. The time dynamics of the $h$ gate are therefore given by

\frac{d h}{d t} = \frac{h_{\infty} - h}{τ_{h}} .

Quantification of single-neuron model performance

Request a detailed protocol

$R^{2}$ was calculated based on the training set $\frac{d V}{d t}$ predicted by the subthreshold component of a given GIF model (Equations 1 and 3, where the spike times $\hat{t}$ were constrained to match the data), excluding a small window around each spike (from 1.5 ms before to 6.5 ms after in 5-HT neurons, and from 1.5 ms before to 4.0 ms after in SOM and mPFC neurons). ${M_{d}}^{*}$ was calculated based on validation set data as previously described by Naud et al., 2011. This metric is defined as

{M_{d}}^{*} = \frac{2 n_{d m}}{{n_{d d}}^{*} + n_{m m}},

where $n_{d m}$ is the number of model-predicted spikes that occur within 8 ms of a spike in the validation data, and ${n_{d d}}^{*}$ and $n_{m m}$ are the corresponding numbers of coincident spikes across sweeps in the validation data and model predictions (where ${n_{d d}}^{*}$ is corrected for small sample bias). ${M_{d}}^{*}$ can be interpreted as the fraction of model-predicted spikes that occur within 8 ms of a spike emitted by a real neuron (the spike timing precision is set to 8 ms by inspecting the relationship between precision and intrinsic reliability [Jolivet et al., 2008]), corrected such that the chance level is 0 and perfect agreement between predicted and observed spikes is 1.

Population models

Request a detailed protocol

DRN network models were constructed by connecting a population of 400 SOM neuron models to a population of 600 5-HT neuron models in a feed-forward arrangement. Population models were bootstrapped by sampling with replacement from a bank of experimentally constrained GIF models. SOM neuron models were randomly connected to 5-HT neuron models with a connection probability of 2%, such that the expected number of GABAergic synapses on each 5-HT neuron model was 8. We used a conductance-based model of GABAergic synapses with a fixed reversal potential of –76.7 mV, conductance of 0.3 nS, and biexponential kinetics with $τ_{r i s e} = 1.44$ ms, $τ_{d e c a y} = 26.0$ ms, and a propagation delay of 2.0 ms.

Simulated 5-HT populations with decreased or increased I_A were generated by setting ${\bar{g}}_{A}$ in all single neuron models to 0 nS or 10 nS, respectively. DRN network models with homogenized SOM neuron populations were created by setting all SOM neuron model parameters to their respective median values from the bank of experimentally constrained single neuron models. Population models in which the adaptation mechanisms of 5-HT and SOM neuron models were swapped were generated by randomly sampling a GIF model of the opposite cell type and substituting in its adaptation filter coefficients ${β_{j}}^{(γ)}_{}$ and $w_{j}$ . This procedure is summarized in the following pseudocode:

for 5-HT_model in 5-HT_population; do

SOM_model = random_choice(SOM_models)

5-HT_model.eta.coefficients= SOM_model.eta.coefficients

5-HT_model.gamma.coefficients= SOM_model.gamma.coefficients

end for.

Numerical methods

Request a detailed protocol

Simulations were implemented in Python and C++ using custom-written extensions of the GIF Fitting Toolbox (Pozzorini et al., 2015; original code archived at https://github.com/pozzorin/GIFFittingToolbox; Pozzorini, 2016). Numerical integration was performed using the Euler method with a time step of 0.1 ms for the GIF model and related models (to match the sampling rate of electrophysiological recordings) and ${0.001 τ}_{m e m}$ for the toy model of a neuron with I_A .

Statistics

Statistical analysis was carried out using the SciPy and statannot (https://github.com/webermarcolivier/statannot; Weber, 2022) Python packages. Non-parametric tests were used for all two-sample comparisons (Mann-Whitney U test for unpaired samples and Wilcoxon signed-rank test for paired samples). Non-parametric tests were chosen because we often had reason to believe that our data did not come from a normal distribution, either due to intrinsic qualities of the data, such as being bounded between 0 and 1, or due to skewness apparent in our samples. Whenever multiple tests were performed in the same figure panel, p-values were adjusted for multiple comparisons using the Bonferroni correction. ‘*’, ‘**’, ‘***’, and ‘****’ are used in figures to denote statistical significance at the p≤0.05, 0.01, 0.001, and 0.0001 levels, respectively, and ‘o’ is used to indicate a trend toward significance (defined as 0.05<p≤0.1). Exact p-values are reported in the main text, and summary statistics are presented as mean ± SD. Sample sizes always refer to biological replicates.

Appendix 1

Definition of aGIF model subthreshold dynamics

The subthreshold dynamics of the aGIF model are given by

C \frac{d V}{d t} = - g_{l} (V (t) - E_{l}) - {\bar{g}}_{A} m_{\infty} (V (t)) h (t) D_{K} (V (t)) - {\bar{g}}_{K} n_{\infty} (V (t)) D_{K} (V (t)) - \sum_{{\hat{t}}_{i} < t} η (t - {\hat{t}}_{i}) + I_{inj} (t),

where $V (t)$ is the membrane voltage, $I_{inj} (t)$ is an externally applied current, $D_{K} = V (t) - E_{K}$ is the driving force on potassium, and $\sum_{{\hat{t}}_{i} < t} η (t - {\hat{t}}_{i})$ is the adaptation current $η$ summated over all past spikes ${{\hat{t}}_{i} \in S : {\hat{t}}_{i} < t}$ ( ${S}$ is the set of all spike times).

The adaptation current produced by a single spike $η (t - {\hat{t}}_{i})$ is implemented as a sum of $k$ exponentials given by

η (t - {\hat{t}}_{i}) = {\begin{cases} \sum_{j = 1}^{k} w_{j} e^{- (t - {\hat{t}}_{i}) / τ_{j}} & if t > {\hat{t}}_{i} \\ 0 & otherwise, \end{cases}

where the timescales $τ_{1}, τ_{2}, \dots, τ_{k}$ are treated as hyperparameters. If we substituted this implementation of $η$ back into Equation A1, the term associated with the spike-triggered current $\sum_{{\hat{t}}_{i} < t} η (t - {\hat{t}}_{i})$ would become $\sum_{{\hat{t}}_{i} < t} \sum_{j = 1}^{k} w_{j} e^{- (t - {\hat{t}}_{i}) / τ_{j}}$ . This double sum can be written more concisely as

\sum_{{\hat{t}}_{i} < t} η (t - {\hat{t}}_{i}) = \sum_{j = 1}^{k} w_{j} {\hat{η}}_{j} (t),

where ${\hat{η}}_{j} (t) = \sum_{{\hat{t}}_{i} < t} e^{- (t - {\hat{t}}_{i}) / τ_{j}}$ is a basis for the adaptation current over the timescale $τ_{j}$ .

Substituting the definition of the adaptation current from Equation A3 into Equation A1, we obtain a detailed definition of the aGIF model subthreshold dynamics as follows:

\frac{d V}{d t} = \frac{1}{C} (- g_{l} (V (t) - E_{l}) - {\bar{g}}_{A} m_{\infty} (t) h (t) D_{K} (t) - {\bar{g}}_{K} n_{\infty} (t) D_{K} (t) - w_{1} {\hat{η}}_{1} (t) - \dots - w_{k} {\hat{η}}_{k} (t) + I_{inj} (t)) .

Estimating aGIF model subthreshold parameters

Given a training dataset $D = {(V (t), I_{inj} (t)) : 1 \leq t \leq T}$ , knowledge of the equilibrium gating functions $m_{\infty}, h_{\infty}, and n_{\infty}$ , and appropriate choices of $τ_{1}, τ_{2}, \dots, τ_{k}$ in $η$ , our goal is to estimate the remaining parameters in Equation A4; namely, $g_{l}, C, E_{l}, {\bar{g}}_{A}, {\bar{g}}_{K}, w_{1}, \dots, w_{k}, and τ_{h}$ , where $τ_{h}$ is the time constant of the inactivation gate $h$ . Fortunately, all of these except for $τ_{h}$ can be estimated easily using linear regression.

We begin by rewriting Equation A4 as the product of a row vector of predictors $x$ and a column vector of coefficients $β$ as follows

\frac{d V}{d t} = x \cdot β

= {[\begin{matrix} V (t) \\ 1 \\ m_{\infty} h D_{K} \\ n_{\infty} D_{K} \\ {\hat{η}}_{1} (t) \\ ⋮ \\ {\hat{η}}_{k} (t) \\ I_{inj} (t) \end{matrix}]}^{⊤} \cdot [\begin{matrix} - g_{l} / C \\ g_{l} E_{l} / C \\ - {\bar{g}}_{A} / C \\ - {\bar{g}}_{K} / C \\ - w_{1} / C \\ ⋮ \\ - w_{k} / C \\ 1 / C \end{matrix}] .

We solve this subject to $g_{l}, C, {\bar{g}}_{A}, {\bar{g}}_{K} \geq 0$ using scipy.optimize.lsq_linear.

Next we turn to the question of calculating all of the components of $x$ . Because Equation A4 only reflects the subthreshold dynamics of the aGIF model, we begin by removing all time points in $D$ within a small window around each spike (from 1.5 ms before each spike until the end of the absolute refractory period). Given the voltages in the cleaned dataset and the set of spike times, it is simple to calculate $m_{\infty}, h_{\infty}, n_{\infty}, D_{K} and {\hat{η}}_{i} (t)$ . To calculate $h$ for each time point, we order the values of $h_{\infty}$ according to time and integrate $h$ numerically using a fixed time step, and the initial condition $h = h_{\infty}$ . This has the effect of assuming that the dynamics of the $h$ gate are paused just before each spike and resumed at the end of the refractory period.

The variance explained by the subthreshold model is a non-convex function of $τ_{h}$ . We therefore conducted a line search over plausible values of $τ_{h}$ and chose the value associated with the highest variance explained. This is equivalent to solving

arg min_{θ} ‖ \hat{\frac{d V}{d t}} - \frac{d V}{d t} ‖

where $θ = (β, τ_{h})$ .

Single-neuron model hyperparameters. For more details on the iGIF model hyperparameter , see Mensi et al., 2016.

Model	Parameter	Symbol	Cell type	Value (ms)
All	$η$ timescales	$τ_{1}, τ_{2}, \dots, τ_{k}$	All	3, 10, 30, 100, 300, 1000, 3000
All	$γ$ timescales	None	All	3, 30, 300, 3000
All	Refractory period	None	5-HT	6.5
All	Refractory period	None	SOM and mPFC	4.0
iGIF	Candidate threshold-coupling timescales	$τ_{θ}$	All	1, 2, 5, 10, 22, 46, 100
aGIF	Candidate inactivation timescales	$τ_{h}$	All	10, 13, 18, 25, 33, 45, 61, 82, 111, 150

Data availability

Raw data is available on Dryad at https://doi.org/10.5061/dryad.66t1g1k2w. Code to fit models, run simulations, and reproduce figures is available at https://github.com/nauralcodinglab/raphegif, (copy archived at swh:1:rev:0a11ab4fe19fa54ddb3f734ad9131d6789b6bed5).

The following previously published data sets were used

1. Harkin EF
2. Lynn M
3. Boucher J
4. Caya-Bissonnette L
5. Cyr D
6. Stewart C
7. Béïque J
(2021) Dryad Digital Repository
Patch-clamp recordings from dorsal raphe neurons.
https://doi.org/10.5061/dryad.66t1g1k2w

References

(2004) Anatomic and functional topography of the dorsal raphe nucleus
Annals of the New York Academy of Sciences 1018:46–57.
https://doi.org/10.1196/annals.1296.005
- PubMed
- Google Scholar
1. Aghajanian GK
2. Vandermaelen CP
(1982) Intracellular recordings from serotonergic dorsal raphe neurons: pacemaker potentials and the effect of LSD
Brain Research 238:463–469.
https://doi.org/10.1016/0006-8993(82)90124-x
- PubMed
- Google Scholar
1. Aghajanian GK
(1985) Modulation of a transient outward current in serotonergic neurones by alpha 1-adrenoceptors
Nature 315:501–503.
https://doi.org/10.1038/315501a0
- PubMed
- Google Scholar
(2007) D2-like dopamine receptors depolarize dorsal raphe serotonin neurons through the activation of nonselective cationic conductance
The Journal of Pharmacology and Experimental Therapeutics 320:376–385.
https://doi.org/10.1124/jpet.106.111690
- PubMed
- Google Scholar
1. Amo R
2. Fredes F
3. Kinoshita M
4. Aoki R
5. Aizawa H
6. Agetsuma M
7. Aoki T
8. Shiraki T
9. Kakinuma H
10. Matsuda M
11. Yamazaki M
12. Takahoko M
13. Tsuboi T
14. Higashijima S
15. Miyasaka N
16. Koide T
17. Yabuki Y
18. Yoshihara Y
19. Fukai T
20. Okamoto H
(2014) The habenulo-raphe serotonergic circuit encodes an aversive expectation value essential for adaptive active avoidance of danger
Neuron 84:1034–1048.
https://doi.org/10.1016/j.neuron.2014.10.035
- PubMed
- Google Scholar
1. Andrade R
(2011) Serotonergic regulation of neuronal excitability in the prefrontal cortex
Neuropharmacology 61:382–386.
https://doi.org/10.1016/j.neuropharm.2011.01.015
- PubMed
- Google Scholar
1. Baraban JM
2. Aghajanian GK
(1981) Noradrenergic innervation of serotonergic neurons in the dorsal raphe: demonstration by electron microscopic autoradiography
Brain Research 204:1–11.
https://doi.org/10.1016/0006-8993(81)90646-6
- PubMed
- Google Scholar
1. Barlow H
(2001) The exploitation of regularities in the environment by the brain
The Behavioral and Brain Sciences 24:602–607.
https://doi.org/10.1017/s0140525x01000024
- PubMed
- Google Scholar
(2004) Serotonergic regulation of membrane potential in developing rat prefrontal cortex: coordinated expression of 5-hydroxytryptamine (5-HT)1A, 5-HT2A, and 5-HT7 receptors
The Journal of Neuroscience 24:4807–4817.
https://doi.org/10.1523/JNEUROSCI.5113-03.2004
- PubMed
- Google Scholar
(2007) Mechanism of the 5-hydroxytryptamine 2A receptor-mediated facilitation of synaptic activity in prefrontal cortex
PNAS 104:9870–9875.
https://doi.org/10.1073/pnas.0700436104
- PubMed
- Google Scholar
1. Benda J
2. Herz AVM
(2003) A universal model for spike-frequency adaptation
Neural Computation 15:2523–2564.
https://doi.org/10.1162/089976603322385063
- PubMed
- Google Scholar
1. Billeh YN
2. Cai B
3. Gratiy SL
4. Dai K
5. Iyer R
6. Gouwens NW
7. Abbasi-Asl R
8. Jia X
9. Siegle JH
10. Olsen SR
11. Koch C
12. Mihalas S
13. Arkhipov A
(2020) Systematic integration of structural and functional data into multi-scale models of mouse primary visual cortex
Neuron 106:388–403.
https://doi.org/10.1016/j.neuron.2020.01.040
- PubMed
- Google Scholar
1. Boureau Y-L
2. Dayan P
(2011) Opponency revisited: competition and cooperation between dopamine and serotonin
Neuropsychopharmacology 36:74–97.
https://doi.org/10.1038/npp.2010.151
- PubMed
- Google Scholar
(2000) Adaptive rescaling maximizes information transmission
Neuron 26:695–702.
https://doi.org/10.1016/s0896-6273(00)81205-2
- PubMed
- Google Scholar
1. Calizo LH
2. Akanwa A
3. Ma X
4. Pan Y-Z
5. Lemos JC
6. Craige C
7. Heemstra LA
8. Beck SG
(2011) Raphe serotonin neurons are not homogenous: electrophysiological, morphological and neurochemical evidence
Neuropharmacology 61:524–543.
https://doi.org/10.1016/j.neuropharm.2011.04.008
- PubMed
- Google Scholar
1. Chance FS
2. Abbott LF
(2000) Divisive inhibition in recurrent networks
Network 11:119–129.
https://doi.org/10.1088/0954-898x_11_2_301
- PubMed
- Google Scholar
1. Cipriani A
2. Furukawa TA
3. Salanti G
4. Chaimani A
5. Atkinson LZ
6. Ogawa Y
7. Leucht S
8. Ruhe HG
9. Turner EH
10. Higgins JPT
11. Egger M
12. Takeshima N
13. Hayasaka Y
14. Imai H
15. Shinohara K
16. Tajika A
17. Ioannidis JPA
18. Geddes JR
(2018) Comparative efficacy and acceptability of 21 antidepressant drugs for the acute treatment of adults with major depressive disorder: a systematic review and network meta-analysis
Focus 16:420–429.
https://doi.org/10.1176/appi.focus.16407
- PubMed
- Google Scholar
(2015) Serotonergic neurons signal reward and punishment on multiple timescales
eLife 4:e06346.
https://doi.org/10.7554/eLife.06346
- PubMed
- Google Scholar
1. Commons KG
(2015) Two major network domains in the dorsal raphe nucleus
The Journal of Comparative Neurology 523:1488–1504.
https://doi.org/10.1002/cne.23748
- PubMed
- Google Scholar
1. Connor JA
2. Stevens CF
(1971) Prediction of repetitive firing behaviour from voltage clamp data on an isolated neurone soma
The Journal of Physiology 213:31–53.
https://doi.org/10.1113/jphysiol.1971.sp009366
- PubMed
- Google Scholar
(1977) Neural repetitive firing: modifications of the hodgkin-huxley axon suggested by experimental results from crustacean axons
Biophysical Journal 18:81–102.
https://doi.org/10.1016/S0006-3495(77)85598-7
- PubMed
- Google Scholar
(2011) Serotonin and dopamine: unifying affective, activational, and decision functions
Neuropsychopharmacology 36:98–113.
https://doi.org/10.1038/npp.2010.121
- PubMed
- Google Scholar
1. Daw ND
2. Kakade S
3. Dayan P
(2002) Opponent interactions between serotonin and dopamine
Neural Networks 15:603–616.
https://doi.org/10.1016/s0893-6080(02)00052-7
- PubMed
- Google Scholar
1. Dayan P
2. Huys QJM
(2009) Serotonin in affective control
Annual Review of Neuroscience 32:95–126.
https://doi.org/10.1146/annurev.neuro.051508.135607
- PubMed
- Google Scholar
1. Dayan P
2. Huys Q
(2015) Serotonin’s many meanings elude simple theories
eLife 4:e07390.
https://doi.org/10.7554/eLife.07390
- PubMed
- Google Scholar
1. Deakin JFW
(1991) Serotonin in panic, anxiety and depression
European Neuropsychopharmacology 1:358–360.
https://doi.org/10.1016/0924-977X(91)90566-D
- Google Scholar
1. Delgado PL
(1994) Serotonin and the neurobiology of depression
Archives of General Psychiatry 51:865.
https://doi.org/10.1001/archpsyc.1994.03950110025005
- PubMed
- Google Scholar
(2013) Genetic approaches for understanding the role of serotonin receptors in mood and behavior
Current Opinion in Neurobiology 23:399–406.
https://doi.org/10.1016/j.conb.2013.01.011
- PubMed
- Google Scholar
1. Dong Y
2. White FJ
(2003) Dopamine D1-class receptors selectively modulate a slowly inactivating potassium current in rat medial prefrontal cortex pyramidal neurons
The Journal of Neuroscience 23:2686–2695.
https://doi.org/10.1523/JNEUROSCI.23-07-02686.2003
- PubMed
- Google Scholar
1. Dong Y
2. Nasif FJ
3. Tsui JJ
4. Ju WY
5. Cooper DC
6. Hu X-T
7. Malenka RC
8. White FJ
(2005) Cocaine-induced plasticity of intrinsic membrane properties in prefrontal cortex pyramidal neurons: adaptations in potassium currents
The Journal of Neuroscience 25:936–940.
https://doi.org/10.1523/JNEUROSCI.4715-04.2005
- PubMed
- Google Scholar
1. Doya K
(2002) Metalearning and neuromodulation
Neural Networks 15:495–506.
https://doi.org/10.1016/s0893-6080(02)00044-8
- PubMed
- Google Scholar
(2015) Ion channel degeneracy enables robust and tunable neuronal firing rates
PNAS 112:E5361–E5370.
https://doi.org/10.1073/pnas.1516400112
- PubMed
- Google Scholar
1. Ermentrout B
(1998) Linearization of F-I curves by adaptation
Neural Computation 10:1721–1729.
https://doi.org/10.1162/089976698300017106
- PubMed
- Google Scholar
1. Fava M
2. Kendler KS
(2000) Major depressive disorder
Neuron 28:335–341.
https://doi.org/10.1016/s0896-6273(00)00112-4
- PubMed
- Google Scholar
1. Ferguson KA
2. Cardin JA
(2020) Mechanisms underlying gain modulation in the cortex
Nature Reviews. Neuroscience 21:80–92.
https://doi.org/10.1038/s41583-019-0253-y
- PubMed
- Google Scholar
(2015) Activation of dorsal raphe serotonergic neurons promotes waiting but is not reinforcing
Current Biology 25:306–315.
https://doi.org/10.1016/j.cub.2014.12.002
- PubMed
- Google Scholar
(2016) Target-specific modulation of the descending prefrontal cortex inputs to the dorsal raphe nucleus by cannabinoids
PNAS 113:5429–5434.
https://doi.org/10.1073/pnas.1522754113
- PubMed
- Google Scholar
1. Gerstner W
(2000) Population dynamics of spiking neurons: fast transients, asynchronous states, and locking
Neural Computation 12:43–89.
https://doi.org/10.1162/089976600300015899
- PubMed
- Google Scholar
1. Gerstner W
2. Naud R
(2009) How good are neuron models?
Science 326:379–380.
https://doi.org/10.1126/science.1181936
- PubMed
- Google Scholar
Book
(2014) Neuronal Dynamics
Cambridge University Press.
https://doi.org/10.1017/CBO9781107447615
- Google Scholar
1. Getting PA
(1983) Mechanisms of pattern generation underlying swimming in tritonia. III. intrinsic and synaptic mechanisms for delayed excitation
Journal of Neurophysiology 49:1036–1050.
https://doi.org/10.1152/jn.1983.49.4.1036
- PubMed
- Google Scholar
Preprint
1. Gonçalves P
2. Lueckmann J
3. Deistler M
4. Nonnenmacher M
5. Öca K
6. Bassetto G
7. Chintaluri C
8. Podlaski W
9. Haddad S
10. Vogels T
11. Greenberg D
12. Macke J
(2019) Training Deep Neural Density Estimators to Identify Mechanistic Models of Neural Dynamics
bioRxiv.
https://doi.org/10.1101/838383
- Google Scholar
1. Grace AA
2. Onn SP
(1989) Morphology and electrophysiological properties of immunocytochemically identified rat dopamine neurons recorded in vitro
The Journal of Neuroscience 9:3463–3481.
https://doi.org/10.1523/JNEUROSCI.09-10-03463.1989
- PubMed
- Google Scholar
(2022) Serotonin neurons modulate learning rate through uncertainty
Current Biology 32:586–599.
https://doi.org/10.1016/j.cub.2021.12.006
- PubMed
- Google Scholar
(2013) Potassium channels control the interaction between active dendritic integration compartments in layer 5 cortical pyramidal neurons
Neuron 79:516–529.
https://doi.org/10.1016/j.neuron.2013.06.005
- PubMed
- Google Scholar
(1997) K+ channel regulation of signal propagation in dendrites of hippocampal pyramidal neurons
Nature 387:869–875.
https://doi.org/10.1038/43119
- PubMed
- Google Scholar
(2006) Efficient estimation of detailed single-neuron models
Journal of Neurophysiology 96:872–890.
https://doi.org/10.1152/jn.00079.2006
- PubMed
- Google Scholar
1. Huys QJM
2. Paninski L
(2009) Smoothing of, and parameter estimation from, noisy biophysical recordings
PLOS Computational Biology 5:e1000379.
https://doi.org/10.1371/journal.pcbi.1000379
- PubMed
- Google Scholar
1. Iigaya K
2. Fonseca MS
3. Murakami M
4. Mainen ZF
5. Dayan P
(2018) An effect of serotonergic stimulation on learning rates for rewards apparent after long intertrial intervals
Nature Communications 9:2477.
https://doi.org/10.1038/s41467-018-04840-2
- PubMed
- Google Scholar
1. Jolivet R
2. Kobayashi R
3. Rauch A
4. Naud R
5. Shinomoto S
6. Gerstner W
(2008) A benchmark test for A quantitative assessment of simple neuron models
Journal of Neuroscience Methods 169:417–424.
https://doi.org/10.1016/j.jneumeth.2007.11.006
- PubMed
- Google Scholar
1. Khaliq ZM
2. Bean BP
(2008) Dynamic, nonlinear feedback regulation of slow pacemaking by A-type potassium current in ventral tegmental area neurons
The Journal of Neuroscience 28:10905–10917.
https://doi.org/10.1523/JNEUROSCI.2237-08.2008
- PubMed
- Google Scholar
1. Kim HR
2. Malik AN
3. Mikhael JG
4. Bech P
5. Tsutsui-Kimura I
6. Sun F
7. Zhang Y
8. Li Y
9. Watabe-Uchida M
10. Gershman SJ
11. Uchida N
(2020) A unified framework for dopamine signals across timescales
Cell 183:1600–1616.
https://doi.org/10.1016/j.cell.2020.11.013
- PubMed
- Google Scholar
1. Knight BW
(1972) Dynamics of encoding in a population of neurons
The Journal of General Physiology 59:734–766.
https://doi.org/10.1085/jgp.59.6.734
- PubMed
- Google Scholar
1. Kohn A
(2007) Visual adaptation: physiology, mechanisms, and functional benefits
Journal of Neurophysiology 97:3155–3164.
https://doi.org/10.1152/jn.00086.2007
- PubMed
- Google Scholar
1. Li Y
2. Zhong W
3. Wang D
4. Feng Q
5. Liu Z
6. Zhou J
7. Jia C
8. Hu F
9. Zeng J
10. Guo Q
11. Fu L
12. Luo M
(2016) Serotonin neurons in the dorsal raphe nucleus encode reward signals
Nature Communications 7:10503.
https://doi.org/10.1038/ncomms10503
- PubMed
- Google Scholar
(2005) Modulation of anxiety circuits by serotonergic systems
Stress 8:233–246.
https://doi.org/10.1080/10253890500492787
- PubMed
- Google Scholar
(2008) Fractional differentiation by neocortical pyramidal neurons
Nature Neuroscience 11:1335–1342.
https://doi.org/10.1038/nn.2212
- PubMed
- Google Scholar
1. Luo M
2. Li Y
3. Zhong W
(2016) Do dorsal raphe 5-HT neurons encode “beneficialness”?
Neurobiology of Learning and Memory 135:40–49.
https://doi.org/10.1016/j.nlm.2016.08.008
- PubMed
- Google Scholar
Preprint
1. Lynn MB
2. Geddes S
3. Chahrour M
4. Maillé S
5. Harkin E
6. Harvey-Girard É
7. Haj-Dahmane S
8. Naud R
9. Béïque JC
(2022) A Slow 5-HT1AR-Mediated Recurrent Inhibitory Network in Raphe Computes Contextual Value through Synaptic Facilitation
bioRxiv.
https://doi.org/10.1101/2022.08.31.506056
- Google Scholar
1. Maier SF
2. Watkins LR
(2005) Stressor controllability and learned helplessness: the roles of the dorsal raphe nucleus, serotonin, and corticotropin-releasing factor
Neuroscience and Biobehavioral Reviews 29:829–841.
https://doi.org/10.1016/j.neubiorev.2005.03.021
- PubMed
- Google Scholar
1. Marder E
(2012) Neuromodulation of neuronal circuits: back to the future
Neuron 76:1–11.
https://doi.org/10.1016/j.neuron.2012.09.010
- PubMed
- Google Scholar
1. Markram H
(2006) The blue brain project
Nature Reviews. Neuroscience 7:153–160.
https://doi.org/10.1038/nrn1848
- PubMed
- Google Scholar
1. Matias S
2. Lottem E
3. Dugué GP
4. Mainen ZF
(2017) Activity patterns of serotonin neurons underlying cognitive flexibility
eLife 6:e20552.
https://doi.org/10.7554/eLife.20552
- PubMed
- Google Scholar
1. Matsumoto M
2. Hikosaka O
(2007) Lateral habenula as a source of negative reward signals in dopamine neurons
Nature 447:1111–1115.
https://doi.org/10.1038/nature05860
- PubMed
- Google Scholar
1. McCormick DA
(1991) Functional properties of a slowly inactivating potassium current in guinea pig dorsal lateral geniculate relay neurons
Journal of Neurophysiology 66:1176–1189.
https://doi.org/10.1152/jn.1991.66.4.1176
- PubMed
- Google Scholar
1. Mejias JF
2. Longtin A
(2012) Optimal heterogeneity for coding in spiking neural networks
Physical Review Letters 108:228102.
https://doi.org/10.1103/PhysRevLett.108.228102
- PubMed
- Google Scholar
1. Mejias JF
2. Longtin A
(2014) Differential effects of excitatory and inhibitory heterogeneity on the gain and asynchronous state of sparse cortical networks
Frontiers in Computational Neuroscience 8:107.
https://doi.org/10.3389/fncom.2014.00107
- PubMed
- Google Scholar
1. Mejias JF
2. Payeur A
3. Selin E
4. Maler L
5. Longtin A
(2014) Subtractive, divisive and non-monotonic gain control in feedforward nets linearized by noise and delays
Frontiers in Computational Neuroscience 8:19.
https://doi.org/10.3389/fncom.2014.00019
- PubMed
- Google Scholar
1. Mensi S
2. Naud R
3. Pozzorini C
4. Avermann M
5. Petersen CCH
6. Gerstner W
(2012) Parameter extraction and classification of three cortical neuron types reveals two distinct adaptation mechanisms
Journal of Neurophysiology 107:1756–1775.
https://doi.org/10.1152/jn.00408.2011
- PubMed
- Google Scholar
(2016) Enhanced sensitivity to rapid input fluctuations by nonlinear threshold dynamics in neocortical pyramidal neurons
PLOS Computational Biology 12:e1004761.
https://doi.org/10.1371/journal.pcbi.1004761
- PubMed
- Google Scholar
1. Miyazaki KW
2. Miyazaki K
3. Tanaka KF
4. Yamanaka A
5. Takahashi A
6. Tabuchi S
7. Doya K
(2014) Optogenetic activation of dorsal raphe serotonin neurons enhances patience for future rewards
Current Biology 24:2033–2040.
https://doi.org/10.1016/j.cub.2014.07.041
- PubMed
- Google Scholar
1. Miyazaki K
2. Miyazaki KW
3. Yamanaka A
4. Tokuda T
5. Tanaka KF
6. Doya K
(2018) Reward probability and timing uncertainty alter the effect of dorsal raphe serotonin neurons on patience
Nature Communications 9:2048.
https://doi.org/10.1038/s41467-018-04496-y
- PubMed
- Google Scholar
(2016) Conditional anterograde tracing reveals distinct targeting of individual serotonin cell groups (B5-B9) to the forebrain and brainstem
Brain Structure & Function 221:535–561.
https://doi.org/10.1007/s00429-014-0924-4
- PubMed
- Google Scholar
(2008) Reward-dependent modulation of neuronal activity in the primate dorsal raphe nucleus
The Journal of Neuroscience 28:5331–5343.
https://doi.org/10.1523/JNEUROSCI.0021-08.2008
- PubMed
- Google Scholar
1. Naud R
2. Gerhard F
3. Mensi S
4. Gerstner W
(2011) Improved similarity measures for small sets of spike trains
Neural Computation 23:3016–3069.
https://doi.org/10.1162/NECO_a_00208
- PubMed
- Google Scholar
1. Naud R
2. Gerstner W
(2012) Coding and decoding with adapting neurons: A population approach to the peri-stimulus time histogram
PLOS Computational Biology 8:e1002711.
https://doi.org/10.1371/journal.pcbi.1002711
- PubMed
- Google Scholar
(2014) Organization of monosynaptic inputs to the serotonin and dopamine neuromodulatory systems
Cell Reports 8:1105–1118.
https://doi.org/10.1016/j.celrep.2014.06.042
- PubMed
- Google Scholar
(2019) Embracing diversity in the 5-HT neuronal system
Nature Reviews. Neuroscience 20:397–424.
https://doi.org/10.1038/s41583-019-0151-3
- PubMed
- Google Scholar
(2005) Comparing integrate-and-fire models estimated using intracellular and extracellular data
Neurocomputing 65–66:379–385.
https://doi.org/10.1016/j.neucom.2004.10.032
- Google Scholar
1. Paquelet GE
2. Carrion K
3. Lacefield CO
4. Zhou P
5. Hen R
6. Miller BR
(2022) Single-cell activity and network properties of dorsal raphe nucleus serotonin neurons during emotionally salient behaviors
Neuron 110:2664–2679.
https://doi.org/10.1016/j.neuron.2022.05.015
- PubMed
- Google Scholar
(2019) Classes of dendritic information processing
Current Opinion in Neurobiology 58:78–85.
https://doi.org/10.1016/j.conb.2019.07.006
- PubMed
- Google Scholar
1. Pillow JW
2. Simoncelli EP
(2006) Dimensionality reduction in neural models: an information-theoretic generalization of spike-triggered average and covariance analysis
Journal of Vision 6:414–428.
https://doi.org/10.1167/6.4.9
- PubMed
- Google Scholar
(2008) Spatio-temporal correlations and visual signalling in a complete neuronal population
Nature 454:995–999.
https://doi.org/10.1038/nature07140
- PubMed
- Google Scholar
1. Pollak Dorocic I
2. Fürth D
3. Xuan Y
4. Johansson Y
5. Pozzi L
6. Silberberg G
7. Carlén M
8. Meletis K
(2014) A whole-brain atlas of inputs to serotonergic neurons of the dorsal and median raphe nuclei
Neuron 83:663–678.
https://doi.org/10.1016/j.neuron.2014.07.002
- PubMed
- Google Scholar
1. Pozzorini C
2. Naud R
3. Mensi S
4. Gerstner W
(2013) Temporal whitening by power-law adaptation in neocortical neurons
Nature Neuroscience 16:942–948.
https://doi.org/10.1038/nn.3431
- PubMed
- Google Scholar
1. Pozzorini C
2. Mensi S
3. Hagens O
4. Naud R
5. Koch C
6. Gerstner W
(2015) Automated high-throughput characterization of single neurons by means of simplified spiking models
PLOS Computational Biology 11:e1004275.
https://doi.org/10.1371/journal.pcbi.1004275
- PubMed
- Google Scholar
Software
1. Pozzorini C
(2016) GIFFittingToolbox, version 5b2e17e
GitHub.

https://github.com/pozzorin/GIFFittingToolbox
(2004) Similar network activity from disparate circuit parameters
Nature Neuroscience 7:1345–1352.
https://doi.org/10.1038/nn1352
- PubMed
- Google Scholar
1. Ranade SP
2. Mainen ZF
(2009) Transient firing of dorsal raphe neurons encodes diverse and specific sensory, motor, and reward events
Journal of Neurophysiology 102:3026–3037.
https://doi.org/10.1152/jn.00507.2009
- PubMed
- Google Scholar
1. Ren J
2. Friedmann D
3. Xiong J
4. Liu CD
5. Ferguson BR
6. Weerakkody T
7. DeLoach KE
8. Ran C
9. Pun A
10. Sun Y
11. Weissbourd B
12. Neve RL
13. Huguenard J
14. Horowitz MA
15. Luo L
(2018) Anatomically defined and functionally distinct dorsal raphe serotonin sub-systems
Cell 175:472–487.
https://doi.org/10.1016/j.cell.2018.07.043
- PubMed
- Google Scholar
(2009) 5-ht(1a) receptor function in major depressive disorder
Progress in Neurobiology 88:17–31.
https://doi.org/10.1016/j.pneurobio.2009.01.009
- PubMed
- Google Scholar
(1997) A neural substrate of prediction and reward
Science 275:1593–1599.
https://doi.org/10.1126/science.275.5306.1593
- PubMed
- Google Scholar
1. Schweighofer N
2. Bertin M
3. Shishida K
4. Okamoto Y
5. Tanaka SC
6. Yamawaki S
7. Doya K
(2008) Low-serotonin levels increase delayed reward discounting in humans
The Journal of Neuroscience 28:4528–4532.
https://doi.org/10.1523/JNEUROSCI.4982-07.2008
- PubMed
- Google Scholar
1. Schweimer JV
2. Ungless MA
(2010) Phasic responses in dorsal raphe serotonin neurons to noxious stimuli
Neuroscience 171:1209–1215.
https://doi.org/10.1016/j.neuroscience.2010.09.058
- PubMed
- Google Scholar
1. Segal M
(1985) A potent transient outward current regulates excitability of dorsal raphe neurons
Brain Research 359:347–350.
https://doi.org/10.1016/0006-8993(85)91448-9
- PubMed
- Google Scholar
1. Soubrié P
(1986) Reconciling the role of central serotonin neurons in human and animal behavior
Behavioral and Brain Sciences 9:319–335.
https://doi.org/10.1017/S0140525X00022871
- Google Scholar
1. Storm JF
(1989) An after-hyperpolarization of medium duration in rat hippocampal pyramidal cells
The Journal of Physiology 409:171–190.
https://doi.org/10.1113/jphysiol.1989.sp017491
- PubMed
- Google Scholar
Book
1. Sutton RS
2. Barto A
(2018)
Reinforcement Learning: An Introduction

The MIT Press.
- Google Scholar
1. Teeter C
2. Iyer R
3. Menon V
4. Gouwens N
5. Feng D
6. Berg J
7. Szafer A
8. Cain N
9. Zeng H
10. Hawrylycz M
11. Koch C
12. Mihalas S
(2018) Generalized leaky integrate-and-fire models classify multiple neuron types
Nature Communications 9:709.
https://doi.org/10.1038/s41467-017-02717-4
- PubMed
- Google Scholar
1. Tops M
2. Russo S
3. Boksem MAS
4. Tucker DM
(2009) Serotonin: modulator of a drive to withdraw
Brain and Cognition 71:427–436.
https://doi.org/10.1016/j.bandc.2009.03.009
- PubMed
- Google Scholar
(2013) Intermediate intrinsic diversity enhances neural population coding
PNAS 110:8248–8253.
https://doi.org/10.1073/pnas.1221214110
- PubMed
- Google Scholar
1. Trulson ME
2. Jacobs BL
(1979) Raphe unit activity in freely moving cats: correlation with level of behavioral arousal
Brain Research 163:135–150.
https://doi.org/10.1016/0006-8993(79)90157-4
- PubMed
- Google Scholar
Preprint
1. Tsuda B
2. Pate SC
3. Tye KM
4. Hava T
(2021) Neuromodulators Enable Overlapping Synaptic Memory Regimes and Nonlinear Transition Dynamics in Recurrent Neural Networks
bioRxiv.
https://doi.org/10.1101/2021.05.31.446462
- Google Scholar
1. Tuckwell HC
2. Penington NJ
(2014) Computational modeling of spike generation in serotonergic neurons of the dorsal raphe nucleus
Progress in Neurobiology 118:59–101.
https://doi.org/10.1016/j.pneurobio.2014.04.001
- PubMed
- Google Scholar
(2018) Global and multiplexed dendritic computations under in vivo-like conditions
Neuron 100:579–592.
https://doi.org/10.1016/j.neuron.2018.08.032
- PubMed
- Google Scholar
(2003) Processing of low-probability sounds by cortical neurons
Nature Neuroscience 6:391–398.
https://doi.org/10.1038/nn1032
- PubMed
- Google Scholar
1. Vandermaelen CP
2. Aghajanian GK
(1983) Electrophysiological and pharmacological characterization of serotonergic dorsal raphe neurons recorded extracellularly and intracellularly in rat brain slices
Brain Research 289:109–119.
https://doi.org/10.1016/0006-8993(83)90011-2
- PubMed
- Google Scholar
1. Varga V
2. Losonczy A
3. Zemelman BV
4. Borhegyi Z
5. Nyiri G
6. Domonkos A
7. Hangya B
8. Holderith N
9. Magee JC
10. Freund TF
(2009) Fast synaptic subcortical control of hippocampal circuits
Science 326:449–453.
https://doi.org/10.1126/science.1178307
- PubMed
- Google Scholar
1. Warden MR
2. Selimbeyoglu A
3. Mirzabekov JJ
4. Lo M
5. Thompson KR
6. Kim S-Y
7. Adhikari A
8. Tye KM
9. Frank LM
10. Deisseroth K
(2012) A prefrontal cortex-brainstem neuronal projection that controls response to behavioural challenge
Nature 492:428–432.
https://doi.org/10.1038/nature11617
- PubMed
- Google Scholar
Software
1. Weber M
(2022) Statannot, version 8de4059
GitHub.

https://github.com/webermarcolivier/statannot
1. Weissbourd B
2. Ren J
3. DeLoach KE
4. Guenthner CJ
5. Miyamichi K
6. Luo L
(2014) Presynaptic partners of dorsal raphe serotonergic and gabaergic neurons
Neuron 83:645–662.
https://doi.org/10.1016/j.neuron.2014.06.024
- PubMed
- Google Scholar
Conference
(2011) A Spiking Neuronal Network Model of the Dorsal Raphe Nucleus
2011 International Joint Conference on Neural Networks.
https://doi.org/10.1109/IJCNN.2011.6033414
- Google Scholar
1. Young SN
2. Smith SE
3. Pihl RO
4. Ervin FR
(1985) Tryptophan depletion causes a rapid lowering of mood in normal males
Psychopharmacology 87:173–177.
https://doi.org/10.1007/BF00431803
- PubMed
- Google Scholar
1. Zhong W
2. Li Y
3. Feng Q
4. Luo M
(2017) Learning and stress shape the reward response patterns of serotonin neurons
The Journal of Neuroscience 37:8863–8875.
https://doi.org/10.1523/JNEUROSCI.1181-17.2017
- PubMed
- Google Scholar
1. Zhou L
2. Liu M-Z
3. Li Q
4. Deng J
5. Mu D
6. Sun Y-G
(2017) Organization of functional long-range circuits controlling the activity of serotonergic neurons in the dorsal raphe nucleus
Cell Reports 18:3018–3032.
https://doi.org/10.1016/j.celrep.2017.02.077
- PubMed
- Google Scholar

Decision letter

Naoshige Uchida

Reviewing Editor; Harvard University, United States
John R Huguenard

Senior Editor; Stanford University School of Medicine, United States
Jochen Roeper

Reviewer; Goethe University Frankfurt, Germany
Hitoshi Morikawa

Reviewer; University of Texas at Austin, United States
Paul Miller

Reviewer; Brandeis University, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Temporal derivative computation in the dorsal raphe network revealed by an experimentally-driven augmented integrate-and-fire modeling framework" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and John Huguenard as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Jochen Roeper (Reviewer #1); Hitoshi Morikawa (Reviewer #2); Paul Miller (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

In this study, Harkin and colleagues examined the biophysical properties of dorsal raphe (DR) serotonin neurons using a combination of whole-cell patch clamp recording in brain slices and biophysical computational modeling. The authors adopted the approach called a generalized integrate-and-fire [aGIF] model, which incorporates a relatively small number of salient biophysical properties (ionic currents) of a specific neuron, and whose parameters are optimized based on voltage dynamics obtained experimentally. First, the authors observe two prominent features in serotonin neurons. The first is a prominent spike frequency adaptation, resulting from after-hyperpolarization potentials. The second is a "kink" in the voltage trace leading up to the first spike, which can be attributed to the effect of A-type potassium currents (I-A). Using the data, the authors constructed aGIF model of serotonin neurons which contains I-A, two Hodgkin-Huxley currents and a stochastic spike generation process. Next, the authors constructed a simple network model using aGIF serotonin and GABA neurons, and show that heterogeneous feedforward inhibition due to heterogeneous electrophysiological properties of local GABA neurons lead to divisive inhibition of serotonin neuron firing (i.e., change in the slope of input-output function). Finally, using a ramp depolarization, the authors found that serotonin neurons encode the temporal derivative of depolarization, i.e., the slope of ramp depolarization. This property can be ascribed to the prominent spike-frequency adaptation observed in serotonin neurons.

The reviewers found the conclusions of this study very interesting and functionally relevant. The reviewers also thought that the experiments and modeling are generally sound and rigorous, and the results are presented clearly. However, the reviewers have identified various concerns and points to be clarified.

1. The slice experiments were carried out in room temperature, which distorts the kinetics of ion channels and pumps relevant for in vivo electrical activity. Ideally, some of the main conclusions should be confirmed at room temperature. In case this is difficult, please explicitly discuss this caveat and potential impacts on conclusions.

2. The series resistance up to 30 Mohm remains uncompensated and in combination with potassium currents in the range of 1 nA generates voltage errors up to 30 mV. In addition, gating properties of A-type currents will be distorted under these recording conditions. Do the main experimental results hold with a more stringent cutoff of the series resistance? Conversely, how robust are the modeling results with respect to the specific gating properties of A-type currents including the inactivation kinetics?

3. Some of the model parameters appear to be non-physiological. For example, the reversal potential of potassium is given as at least 20mV lower than that used typically in the literature. Please justify the use of -101mV rather than the range -75mV to -85mV as is more common. Please clarify whether the results depend on this very low reversal potential.

4. It is unlikely that the derivative-like computation holds in all of potential input regimes. It is important to clarify in what conditions derivative-like computation is a good approximation of the input-output relationship of serotonin neurons, and when it breaks down.

5. The derivative-like computation and its role in reinforcement learning are emphasized in the manuscript (e.g. abstract). However, as discussed above, the derivative-like computation likely occurs only in a limited input regimen, and the relevance to reinforcement learning is very speculative. The reviewers thought that these points need to be a little toned down.

Reviewer #1 (Recommendations for the authors):

The electrophysiological characterization of 5-HT neurons has a number of shortcomings that should be addressed by additional experiments:

1. Experiments are carried out in room temperature. This distorts the kinetics of ion channels and pumps relevant for in vivo electrical activity.

3. An intracellular solution is used without calcium buffering. Given the key role of adaptation and relevant calcium-dependent conductances in 5-HT neurons this is not state-of-the-art

4. 4-AP is not a selective Kv4 (IA) channel blocker. There are more-selective, commercially available alternatives like Heteropodatoxin-2 or AmmTx3.

5. Author show that augmentation of the GIF-model with IA-currents does not specifically improve the model in comparison to a generic iGIF model (Figure 3F2). Doesn´t this call in into question the whole strategy of using experimentally validated biophysical enhancements of GIFs? Does the experimentally augmented GIF outperform the standard computational model of 5-HT neurons (see Tuckwell & Penington (2014) Computational modeling of spike generation in serotonergic neurons of the dorsal raphe nucleus. Prog Neurobiol 118:59-101).

6. Authors should experimentally study the biophysics of adaptation in 5-HT neurons as this is most related to temporal derivative computation, the key insight of this ms. For experimental conditions see 3.

Reviewer #2 (Recommendations for the authors):

1) Electrophysiological recordings are conducted at room temperature. It is not clear why the authors did not conduct recordings at more physiological temperature, as temperature should affect properties on many ion channels, and thus the excitability of neurons.

2) The authors may want to briefly mention the potential behavioral conditions in which endocannabinoids could be released in the dorsal raphe.

3) Dopamine neurons are also known to display low frequency firing, large afterhyperpolarizations, and A-type K currents. It would be helpful to briefly discuss whether similar computations could be applied to these neurons, especially in light of their role in reinforcement learning,

4) P. 14: Figure 1C3 is not present in Figure 1. There might be other places where figure panels are not correctly referred to.

Reviewer #3 (Recommendations for the authors):

The paper is very well written and clear, with a solid approach.

p.10 the reversal potential of potassium is given as at least 20mV lower than any I have come across. Please justify the use of -101mV rather than the range -75mV to -85mV as is more common.

p.11 "within 8 ms": this seems like an arbitrary time range and a rather large one to predict spike timing given typical membrane time constants and jitter. Also, it is unclear if this is related to the "1.5 ms before to 6.5 ms after" a spike where the dV/dt is constrained – and again, why this 8 ms window? And why the smaller window for SOM neurons for the fitting of dV/dt, though not the validation of spike times?

p.21-p.22 The authors state that the aGIF is more parsimonious than the iGIF but it appears that the aGIF has more parameters so the reverse would be true? Perhaps because the aGIF includes a known, well-characterized current, the number of free parameters is smaller even if the total number is greater? Or am I missing something? For sure if the current is shown to be present in the neuron then it is reasonable to use a model with that current rather than the iGIF, but "more parsimonious" would not be the reason.

The authors then on p.22 state the aGIF "best accounts" for the data, but again the iGIF appears to fit the data just as well, and no other model currents or neural were tested. I think the valid conclusion is that adding the extra current(s) improves the fit, no more. It would be good to see a criterion like AIC used to justify the improvement to fit is greater than that bound to arise with extra free parameters.

Figure 5: Given the 5-HT neuron has a sustained response that depends on step height, its output is not the differential of its input. This is an important proviso in any discussion of the properties of the circuit as a differentiator, since it is not one.

p.32, Figure 6 caption. "includes a strong subtractive component (F2)". I think "E2" may be meant as I see a shift between "Het SOM" and "Hom SOM" curves in E2 but not in F2. Also it is important to be clear when discussing the "input-output function" of a neuron whether the transient peak response is being discussed, or the sustained response, or both. Normally sustained firing rate to sustained input would be meant, so if it is the peak of the transient, that should be added in the description.

p.33 Figure 7 caption. "output is … linearly related": again, only the peak of the transient response appears to have a linear relationship. This is not the same as the neuron's output in general.

It would be nice to see that the output of the neuron responding to some range of fluctuating inputs, I(t), has higher linear correlation with dI/dt than with I or to show the neuron responds more to dI/dt than I(t) via some other statistical test. In particular, if dI/dt is more and more negative, does the rate decrease more and more (or is it only a positive slope detector)?

Figure 7F is too opaque for me to understand. It seems like two different manipulations on one plot, but then it is unclear why the different manipulations and color scales correspond to different regions of the x-y plane. I am confused, so I suggest a bit more explanation, or 2 figures if 2 manipulations are carried out.

p.39 temporal difference learning for example is not the same as producing the time-derivative of an input. It is the difference between two inputs (one being expected reward the other being actual reward) – or one can calculate it as the difference across trials from some average reward to the current reward, but that difference between the current input and trial-averaged input is over a far slower timescale than that of the peak responses to a change in input demonstrated here.

https://doi.org/10.7554/eLife.72951.sa1

Author response

Essential revisions:

In this study, Harkin and colleagues examined the biophysical properties of dorsal raphe (DR) serotonin neurons using a combination of whole-cell patch clamp recording in brain slices and biophysical computational modeling. The authors adopted the approach called a generalized integrate-and-fire [aGIF] model, which incorporates a relatively small number of salient biophysical properties (ionic currents) of a specific neuron, and whose parameters are optimized based on voltage dynamics obtained experimentally. First, the authors observe two prominent features in serotonin neurons. The first is a prominent spike frequency adaptation, resulting from after-hyperpolarization potentials. The second is a "kink" in the voltage trace leading up to the first spike, which can be attributed to the effect of A-type potassium currents (I-A). Using the data, the authors constructed aGIF model of serotonin neurons which contains I-A, two Hodgkin-Huxley currents and a stochastic spike generation process. Next, the authors constructed a simple network model using aGIF serotonin and GABA neurons, and show that heterogeneous feedforward inhibition due to heterogeneous electrophysiological properties of local GABA neurons lead to divisive inhibition of serotonin neuron firing (i.e., change in the slope of input-output function). Finally, using a ramp depolarization, the authors found that serotonin neurons encode the temporal derivative of depolarization, i.e., the slope of ramp depolarization. This property can be ascribed to the prominent spike-frequency adaptation observed in serotonin neurons.

The reviewers found the conclusions of this study very interesting and functionally relevant. The reviewers also thought that the experiments and modeling are generally sound and rigorous, and the results are presented clearly. However, the reviewers have identified various concerns and points to be clarified.

1. The slice experiments were carried out in room temperature, which distorts the kinetics of ion channels and pumps relevant for in vivo electrical activity. Ideally, some of the main conclusions should be confirmed at room temperature. In case this is difficult, please explicitly discuss this caveat and potential impacts on conclusions.

This is a valid and admittedly important point. Because this issue was presented as perhaps the chief concern related to our experimental work, we have spent considerable time and effort to experimentally address it. Experiments were initially carried out at room temperature simply because of our difficulty in obtaining high quality and stable recordings at a higher temperature. To address the possibility that some of our results could be an artifact of this experimental limitation, we optimized some slicing procedures and had a more experienced electrophysiologist in our group carry out key recordings closer to physiological temperature. We managed to carry out recordings at 29°C to 30°C measured in the recording bath. In brief, we found that the voltage-dependence and kinetics of IA were indeed slightly altered at increased temperatures, and that the very strong adaptation observed in serotonin neurons was somewhat attenuated due to a reduction in the spike-triggered threshold movement. These quantitative changes prompted us to revisit our simulations seeking to determine the output of the DRN in response to inputs of different derivatives. However, these simulations revealed that increasing temperature was not accompanied by substantive qualitative changes to our main conclusions (the simulations are different though). Since we believe other readers may share this important concern, we have added several supplementary figures to our manuscript that directly illustrate this point (i.e., Figures 1S5, 2S2, 2S3, 3S3, 4S2, 5S5, and 7S1, results reported on pages 7, 15, 18, 20, and 26).

2. The series resistance up to 30 Mohm remains uncompensated and in combination with potassium currents in the range of 1 nA generates voltage errors up to 30 mV. In addition, gating properties of A-type currents will be distorted under these recording conditions. Do the main experimental results hold with a more stringent cutoff of the series resistance? Conversely, how robust are the modeling results with respect to the specific gating properties of A-type currents including the inactivation kinetics?

We thank the reviewers for explicitly noting this point and we wholeheartedly agree that the series resistance cutoff used in our initial analysis was too charitable (especially for the present case where we are trying to clamp large and fast currents). This cutoff was used for early and still expanding datasets and somehow got carried forward. Upon closer inspection prompted by this comment, we realized that the overwhelming majority of the voltage-clamp recordings used to characterize IA in 5-HT neurons had a significantly lower series resistance (14.7 +/- 6.2 MOhm, mean +/- SD). A standard series resistance cutoff of 20 MOhm ended up excluding only three of the thirteen original recordings (series resistance of the ten remaining recordings: 12.1 +/- 4.0 MOhm). We have added a note about these values on page 37 of our revised manuscript.

Do the main experimental results hold using the more stringent cutoff?

Yes. The estimates of IA peak amplitude, time to peak, and decay time constant using this more rigorous inclusion criterion are only marginally different from the original estimates (see table below). We agree with the reviewer that a more stringent cutoff would be preferable. Unfortunately, this would involve rerunning all of our simulations involving the aGIF (which is a significant portion of this study). The absence of a significant change in the kinetics strongly indicates that rerunning the simulations would not produce different results. We have made a note of the series resistance values (above) and their effects on the estimated parameters of IA ( Author response table 1) in the manuscript and kindly propose to keep the original series resistance cutoff.

Author response table 1

Parameter	Original estimate	Updated estimate
Amplitude (pA)	929±69	937±79
Time to peak (ms)	7.46±0.21	7.42±0.27
Time constant (ms)	42.9±2.6	41.0±3.1

Conversely, how robust are the modeling results with respect to the specific gating properties of A-type currents including the inactivation kinetics?

The modeling results do not depend on our experimental measurements of IA’s inactivation kinetics and are expected to be robust to Ra-related distortion of the equilibrium gating curves. The features of IA most important to our modeling results are that the current both activates and inactivates below spike threshold and that the inactivation time constant is not very small relative to the membrane time constant. If the ‘true’ activation and inactivation gating curves were located further to the left or rise more rapidly than the ones we report due to series resistance-related voltage errors, this current would still activate and inactivate below threshold, leaving unaltered our main conclusions. Moreover, and importantly, because the inactivation time constants of IA used in our simulations are estimated directly from the fluctuating input data as part of the GIF model fitting procedure, these values represent effective quantities that take into account the impact of IA on the subthreshold dynamics of 5-HT neurons and are not sensitive to the same sources of error as our equilibrium gating curves. (The values of the inactivation time constant estimated by the aGIF model are 53 +/- 42 ms (mean +/- SD) with an interquartile range of 27 to 77 ms (N=18), consistent with the values obtained using voltage clamp experiments). Lastly, the simulations shown in Figure 2S1 show empirically that reducing the effective inactivation time constant of IA by more than 40% does not qualitatively alter our results. Altogether, the modeling results robustly hold up to variability of gating parameters significantly beyond that expected for series resistance-mediated distortion.

3. Some of the model parameters appear to be non-physiological. For example, the reversal potential of potassium is given as at least 20mV lower than that used typically in the literature. Please justify the use of -101mV rather than the range -75mV to -85mV as is more common. Please clarify whether the results depend on this very low reversal potential.

Thank you for the opportunity to clarify this aspect of our methodology. The low potassium reversal was initially chosen to match our experimental conditions to ensure accurate model parameter estimates (we have noted this point on page 41 of our revised manuscript). Since the most important results that emerge from our simulations are related to spike timing and firing rate, which occur within the voltage range above the physiological reversal potential of potassium, raising the reversal potential of potassium would have the exact same effect as reducing the maximal conductance associated with IA (by approx. 25% to 50% depending on the reversal potential used). We have extensively explored the effects of changing this parameter in Figures 2B, 6A, 7E, and 7G, and therefore changing the reversal potential of potassium would not significantly affect our conclusions. Nonetheless, to empirically support this idea, we adjusted the potassium reversal potential to -89.1 mV (calculated using the Nernst formula and physiological potassium concentrations of [K]o = 5 mM, [K]i = 140 mM) in our aGIF models fitted to the data collected near physiological temperature. As noted in our response to point #1, the results obtained using these models are consistent with our main conclusions.

4. It is unlikely that the derivative-like computation holds in all of potential input regimes. It is important to clarify in what conditions derivative-like computation is a good approximation of the input-output relationship of serotonin neurons, and when it breaks down.

This is an important point, especially considering that 5-HT neurons can exhibit sustained firing even when the derivative of the input is zero (e.g., Figure 1S1 and 5). Our results suggest that the DRN is best conceptualized as encoding both the intensity and temporal-derivative of its input. We have added a toy rate-based model of the DRN to capture this idea and to clarify its limitations (Figure 7S2). Reassuringly, we found that the derivative-like computation dominates DRN output precisely when the derivative of the input is large (and positive), but we have amended the results and Discussion sections to clarify when and why this computation breaks down (pages 26, 27 and 34). We have also changed our wording of this result throughout the manuscript to emphasize that the DRN encodes a mixture of the intensity and temporal derivative of its input (e.g., abstract, page 5).

5. The derivative-like computation and its role in reinforcement learning are emphasized in the manuscript (e.g. abstract). However, as discussed above, the derivative-like computation likely occurs only in a limited input regimen, and the relevance to reinforcement learning is very speculative. The reviewers thought that these points need to be a little toned down.

Thank you for this feedback: we entirely agree that the role of this computation in RL is at this time purely speculative. After deliberation, we removed the explicit link to RL from the abstract and significantly rewrote and toned down the Discussion section on this speculative point (see pages 34, 35).

Reviewer #1 (Recommendations for the authors):

The electrophysiological characterization of 5-HT neurons has a number of shortcomings that should be addressed by additional experiments:

1. Experiments are carried out in room temperature. This distorts the kinetics of ion channels and pumps relevant for in vivo electrical activity.

We have experimentally addressed this point. Please see our response to essential point #1 above.

2. The series resistance up to 30Mohm remains uncompensated and in combination with potassium currents in the range of 1 nA generates voltage errors up to 30 mV. In addition, gating properties of IA-currents will be distorted under these recording conditions

Please see our response to essential revision #2.

3. An intracellular solution is used without calcium buffering. Given the key role of adaptation and relevant calcium-dependent conductances in 5-HT neurons this is not state-of-the-art

This is an important point. 5-HT neurons are indeed endowed with several Calcium-dependent conductances that are bound to regulate their firing behaviors. As a mere example, 5-HT neurons are well known to exhibit a prominent calcium- and SK channel dependent AHP (Scuvée-Moreau et al., Br. J. Pharm, 2004; Sargin et al., eLife 2016). There are likely other parameters of ionic conductances expressed in 5-HT neurons that can be regulated by intracellular Calcium. As such, it is by design that we chose to omit Calcium chelators from our intracellular solution in order to have recording conditions that are as physiological as possible. In this line, the GIF formalism is designed to extract excitability parameters from spiking in such naturalistic conditions, and is thus aptly applied in conditions without Calcium buffering. Moreover, because Calcium buffering would alter the contribution of a potentially unknown set of (calcium-dependent) conductances to these spiking behaviors, these experiments would be fundamentally challenging to interpret.

4. 4-AP is not a selective Kv4 (IA) channel blocker. There are more-selective, commercially available alternatives like Heteropodatoxin-2 or AmmTx3.

We were initially primarily concerned with an algorithmic description of this conductance, and not so much with its molecular identification. That said, we followed this reviewer’s suggestion and found that Heteropodatoxin-2 and AmmTx3 block most of the transient component of the outward current readily observed in 5-HT neurons. These results are now in the revised manuscript (Figure 1S4). It is noteworthy that the delay in obtaining these particular drugs are in large part responsible for delaying the submission of this revised manuscript.

5. Author show that augmentation of the GIF-model with IA-currents does not specifically improve the model in comparison to a generic iGIF model (Figure 3F2). Doesn´t this call in into question the whole strategy of using experimentally validated biophysical enhancements of GIFs? Does the experimentally augmented GIF outperform the standard computational model of 5-HT neurons (see Tuckwell & Penington (2014) Computational modeling of spike generation in serotonergic neurons of the dorsal raphe nucleus. Prog Neurobiol 118:59-101).

Thank you for the opportunity to clarify these points. Concerning aGIF vs iGIF, although both models predict the spike trains of 5-HT neurons with comparable accuracy, we would like to emphasize that the aGIF model predicts the subthreshold voltage of these cells significantly more accurately than the iGIF model (Figure 3F1). Although accurately predicting the subthreshold voltage is certainly not always necessary (for example, the GLM and LNL-based approaches discussed on page 31, which are not viable for 5-HT neurons due to their very low firing rates, do not model subthreshold voltage at all), we contend that this represents an important feature in cases where the primary interest lies in exploring the causal relationships between mechanisms operating in subtreshold regimes and spiking behaviors. We have briefly noted this point on page 31 of our revised manuscript.

Beyond this point, the fact that aGIF and iGIF models achieve similar spike timing prediction using squarely different mechanisms has two possible interpretations: either both models capture the same variance (they explain the same phenomenon using different mechanisms) or the models complement one another (they capture different phenomena). Although we believe that the flexibility of the iGIF model makes the former more likely, if it were the latter, this would suggest that GIF models should be gradually augmented with various mechanisms, in a manner indeed reminiscent of Hodgkin-Huxely-type modelling.

There remains, however, an important drawback to biophysically-detailed Hodgkin-Huxley-style models, such as the model of Tuckwell and Penington (2014). Off-the-shelf versions of these models typically predict about half of the variance predicted by GIF models (Gerstner and Naud, Science 2009). Fitting the maximal conductances to data improves the prediction, typically reaching a performance on par with GIF models. Importantly, however, the fitting of biophysical models is exceedingly difficult and may be impossible on limited amounts of data. The fitting process for biophysically-detailed models also uses about 1000 times more computer time. As a result, even if GIF-model-based approaches appear more and more complex, they remain fundamentally tractable such that they can be calibrated on typical datasets (as we have done in this paper).

6. Authors should experimentally study the biophysics of adaptation in 5-HT neurons as this is most related to temporal derivative computation, the key insight of this ms. For experimental conditions see 3.

One of the main contributions of our GIF approach was to estimate the functional amplitude and kinetics of adaptation in 5-HT neurons directly from their spiking and subthreshold responses to in vivo-like fluctuating input (Figure 4). Somewhat surprisingly, this approach demonstrated that an increase in threshold contributed to adaptation in these neurons, in addition to the AHP that has previously been well described in these neurons. The procedure used by GIF-type models to estimate adaptation is well-established (e.g., Mensi, Naud, et al., 2012, Pozzorini, Naud, et al., 2013) and does not require independent characterization of the underlying biophysical mechanisms. This is entirely different from the situation that arose with IA, where augmenting the GIF with this current required separate experiments to characterize its voltage-dependence (as briefly noted in our discussion on page 30). In an independent study, we have carried out voltage-clamp recordings that characterize IAHP in some depth (e.g., modulation by apamin), but similar experiments have previously been reported and it is unclear to us how these results would be integrated into our approach (since the GIF fitting procedure would not use the experimentally-determined parameter values). Furthermore, it is not immediately apparent what type of experiments would be required to fully parameterize the dynamic change in threshold that also contributes to adaptation in these cells. Altogether, we respectfully submit that while the proposed experiments are certainly interesting, it is not clear how they would be incorporated into our model-based approach and they are therefore not required to support the main message of our work.

Reviewer #2 (Recommendations for the authors):

1) Electrophysiological recordings are conducted at room temperature. It is not clear why the authors did not conduct recordings at more physiological temperature, as temperature should affect properties on many ion channels, and thus the excitability of neurons.

We have experimentally addressed this point; please see our response to essential point #1.

2) The authors may want to briefly mention the potential behavioral conditions in which endocannabinoids could be released in the dorsal raphe.

There is an expanding literature that is outlining with increasing mechanistic and molecular detail how endocannabinoids are release and retrogradely transmitted in the DRN (e.g., Haj-Dahmane, PNAS 2018) to affect synaptic transmission (e.g., Haj-Dahmane et al., JPET 2009). We have participated in this research effort by showing that glutamatergic synapses onto DRN GABAergic neurons are far more sensitive to cannabinoid-mediated inhibition than their counterpart onto 5-HT neuron (Geddes et al., PNAS 2016; this finding was one of the important motivating factor for the present study). Despite these advances, the conditions that lead to endocannabinoids release in the DRN during behavior are far less understood. For instance, while it has been shown that 5-HT neuronal activity per se regulates tonic release of endocannabinoid in the DRN, it is unknown whether this is also true for GABAergic neurons. Thus, beyond these broad strokes of guiding principles governing endocannabinoid release in the DRN, our knowledge is very superficial and remains mainly speculative. That said, the homogeneous rise of cannabinoid receptor activation induced by recreational cannabis use is likely to be broadly recapitulated by our modeling. We made a mention on page 34 of the revised manuscript.

3) Dopamine neurons are also known to display low frequency firing, large afterhyperpolarizations, and A-type K currents. It would be helpful to briefly discuss whether similar computations could be applied to these neurons, especially in light of their role in reinforcement learning,

We thank the reviewer for this point, this is indeed the case and increases the scope of our work. We have added a remark to that effect in the Discussion section on reinforcement learning (page 35).

4) P. 14: Figure 1C3 is not present in Figure 1. There might be other places where figure panels are not correctly referred to.

Thank you for pointing this out! We have corrected this error and reviewed the manuscript for similar issues.

Reviewer #3 (Recommendations for the authors):

The paper is very well written and clear, with a solid approach.

p.10 the reversal potential of potassium is given as at least 20mV lower than any I have come across. Please justify the use of -101mV rather than the range -75mV to -85mV as is more common.

Please see our response to essential revision #3.

p.11 "within 8 ms": this seems like an arbitrary time range and a rather large one to predict spike timing given typical membrane time constants and jitter. Also, it is unclear if this is related to the "1.5 ms before to 6.5 ms after" a spike where the dV/dt is constrained – and again, why this 8 ms window? And why the smaller window for SOM neurons for the fitting of dV/dt, though not the validation of spike times?

The 8ms precision used to calculate the Md* spiketrain similarity metric is admittedly somewhat arbitrary. Previous work carried out in cortical neurons used a precision of 4ms. This number was based on a systematic comparison of the spike timing precision and the intrinsic reliability (the consistency of spike timing across repetitions for a given spike timing precision) (Jolivet et al. J Neurosci Meth 2008). We have chosen Δ according to a similar assay in DRN neurons and added a note to this effect on page 42 of our revised manuscript.

The precision used in the Md* metric is unrelated to the “1.5 ms before to 6.5 ms after” window around each spike during which the voltage is ignored. Because integrate-and-fire models including GIF-class models do not explicitly model the membrane voltage during a spike, it was necessary to remove spikes from our data before attempting to fit the subthreshold components of our models. We define the timing of a spike as the instant that the membrane voltage rises past 0 mV, and creating a window of several milliseconds around this point ensures that the entire action potential is removed. The duration of this window varies slightly across cell types due to differences in action potential shape.

p.21-p.22 The authors state that the aGIF is more parsimonious than the iGIF but it appears that the aGIF has more parameters so the reverse would be true? Perhaps because the aGIF includes a known, well-characterized current, the number of free parameters is smaller even if the total number is greater? Or am I missing something? For sure if the current is shown to be present in the neuron then it is reasonable to use a model with that current rather than the iGIF, but "more parsimonious" would not be the reason.

Thank you for the opportunity for clarification. The “non-linear coupling term” that differentiates the GIF and iGIF models mentioned in our results text actually consists of a non-parametric function of past voltage. The non-parametric definition of this term introduces roughly ten additional free parameters compared with the GIF and aGIF models (our aGIF model has only three additional free parameters compared with the GIF model), and, as with non-parametric models in general, the fitted values of these parameters can be difficult to interpret. We believe that this makes the aGIF model more parsimonious and have added clarification to the methods and Results sections to emphasize this point (pages 39 and 14 in the revised manuscript).

The authors then on p.22 state the aGIF "best accounts" for the data, but again the iGIF appears to fit the data just as well, and no other model currents or neural were tested. I think the valid conclusion is that adding the extra current(s) improves the fit, no more. It would be good to see a criterion like AIC used to justify the improvement to fit is greater than that bound to arise with extra free parameters.

Please see point no. 5 from reviewer 1 for changes related to this point. In writing the manuscript, we were conscious of the equal predictive performance of the aGIF and iGIF (specifically with respect to spike-timing; the aGIF better predicts the subthreshold dynamics of 5-HT neurons), which led us to attempt to make the more limited claim that the aGIF model best accounts for the biophysical mechanisms responsible for shaping the spiking behaviour of 5-HT neurons rather than the data (see page 15). Regarding the AIC, however, we respectfully disagree as this methodology was invented for cases where data that was not used for model fitting is not available (i.e., test or validation data). Our data collection approach was designed to collect a separate validation dataset and obviate the need for statistical tools such as AIC (Figure 3C and D). Cross-validation remains a more accurate approach to model selection that takes into account overfitting.

Figure 5: Given the 5-HT neuron has a sustained response that depends on step height, its output is not the differential of its input. This is an important proviso in any discussion of the properties of the circuit as a differentiator, since it is not one.

We believe that DRN is indeed best conceptualized as encoding a mixture of the intensity and temporal derivative of its input. We have added a toy model to formalize this intuition (Figure 7S2) and have adjusted the manuscript text accordingly. (See our responses to other points.)

p.32, Figure 6 caption. "includes a strong subtractive component (F2)". I think "E2" may be meant as I see a shift between "Het SOM" and "Hom SOM" curves in E2 but not in F2. Also it is important to be clear when discussing the "input-output function" of a neuron whether the transient peak response is being discussed, or the sustained response, or both. Normally sustained firing rate to sustained input would be meant, so if it is the peak of the transient, that should be added in the description.

Thank you for pointing this out, we have corrected the issue.

p.33 Figure 7 caption. "output is … linearly related": again, only the peak of the transient response appears to have a linear relationship. This is not the same as the neuron's output in general.

Agreed. We have adjusted the wording of the caption to specify the peak firing rate and added Figure 7S2 to support our broader point.

It would be nice to see that the output of the neuron responding to some range of fluctuating inputs, I(t), has higher linear correlation with dI/dt than with I or to show the neuron responds more to dI/dt than I(t) via some other statistical test. In particular, if dI/dt is more and more negative, does the rate decrease more and more (or is it only a positive slope detector)?

Thank you for the suggestion. We had initially considered taking a similar approach, but were concerned that quantifying the correlation between the output and dI/dt across multiple timescales and levels of background input would introduce sufficient analytical complexity to alienate part of our target readership (and perhaps ourselves). To add clearer support for our derivative-encoding argument, we have fitted a toy rate model of DRN output as a function of the input and its derivative to the data presented in Figure 7 and compared it with a restricted model that does not take the derivative of the input into account (Figure 7S2). We hope that the stark differences between the predictions made by the two toy models will be an adequate compromise to address this point.

Figure 7F is too opaque for me to understand. It seems like two different manipulations on one plot, but then it is unclear why the different manipulations and color scales correspond to different regions of the x-y plane. I am confused, so I suggest a bit more explanation, or 2 figures if 2 manipulations are carried out.

The goal of this figure was to show that two of the biophysical features of 5-HT neurons contribute to shaping the input-output function of the DRN under distinct input regimes. We have updated the caption to emphasize this point and hope that it is now clearer.

p.39 temporal difference learning for example is not the same as producing the time-derivative of an input. It is the difference between two inputs (one being expected reward the other being actual reward) – or one can calculate it as the difference across trials from some average reward to the current reward, but that difference between the current input and trial-averaged input is over a far slower timescale than that of the peak responses to a change in input demonstrated here.

Thank you for raising this point. We have made changes to that effect, see essential point 5.

https://doi.org/10.7554/eLife.72951.sa2

Article and author information

Author details

Emerson F Harkin

Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada

Contribution
Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0698-5894
Michael B Lynn

Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada

Contribution
Conceptualization, Investigation, Methodology

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0760-4555
Alexandre Payeur
1. Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada
2. Department of Physics, University of Ottawa, Ottawa, Canada
Present address
Mila, Université de Montréal, Montréal, Canada

Contribution
Conceptualization, Formal analysis, Validation, Investigation, Methodology

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-2437-8249
Jean-François Boucher

Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada

Contribution
Data curation, Investigation

Competing interests
No competing interests declared
Léa Caya-Bissonnette

Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada

Contribution
Investigation

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-2893-6949
Dominic Cyr

Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada

Contribution
Investigation

Competing interests
No competing interests declared
Chloe Stewart

Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada

Contribution
Investigation

Competing interests
No competing interests declared
André Longtin
1. Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada
2. Department of Physics, University of Ottawa, Ottawa, Canada
Contribution
Conceptualization, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0678-9893
Richard Naud
1. Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada
2. Department of Physics, University of Ottawa, Ottawa, Canada
Contribution
Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Methodology, Writing – review and editing

For correspondence
rnaud@uottawa.ca

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7383-3095
Jean-Claude Béïque

Brain and Mind Research Institute, Centre for Neural Dynamics, Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, Canada

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Validation, Writing – review and editing

For correspondence
jbeique@uottawa.ca

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7278-4906

Funding

Canadian Institutes of Health Research (175325)

Richard Naud
Jean-Claude Béïque

Canadian Institutes of Health Research (175319)

Richard Naud
Jean-Claude Béïque

Natural Sciences and Engineering Research Council of Canada (06972)

Richard Naud

Canada Foundation for Innovation, Brain Canada (Canadian Neurophotonic Platform)

Jean-Claude Béïque

Krembil Foundation

Jean-Claude Béïque

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was carried out on the unceded and unsurrendered land of the Algonquin Anishinaabe people. We thank Dr. Simon Chen for providing SOM-Cre::Rosa-TdTomato mice and Liwen Cai and David Lemelin for technical assistance with mouse lines. We would also like to thank Sébastien Maillé for contributions to pilot-testing of the GIF model in 5-HT neurons and all other members of the J-CB. and RN labs for helpful discussions. EH is thankful to have received graduate scholarships from the Canadian Institutes of Health Research and the Natural Sciences and Engineering Research Council of Canada. This work was supported by grants from the Canadian Institutes of Health Research, the Natural Sciences and Engineering Research Council of Canada, the Canada Foundation for Innovation, Brain Canada (Canadian Neurophotonic Platform), and the Krembil Foundation.

Ethics

All experiments were carried out in accordance with procedures approved by the University of Ottawa Animal Care and Veterinary Services (protocol numbers CMM-164, CMM-176, CMM-1711, CMM-1743, and CMM-2737). At the beginning of each experiment, animals were deeply anesthetized using isofluorane to minimize suffering before being euthanized.

Senior Editor

John R Huguenard, Stanford University School of Medicine, United States

Reviewing Editor

Naoshige Uchida, Harvard University, United States

Reviewers

Jochen Roeper, Goethe University Frankfurt, Germany
Hitoshi Morikawa, University of Texas at Austin, United States
Paul Miller, Brandeis University, United States

Version history

Preprint posted: June 27, 2021 (view preprint)
Received: August 10, 2021
Accepted: December 19, 2022
Accepted Manuscript published: January 19, 2023 (version 1)
Version of Record published: February 27, 2023 (version 2)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.