Temporal derivative computation in the dorsal raphe network revealed by an experimentally driven augmented integrateandfire modeling framework
Abstract
By means of an expansive innervation, the serotonin (5HT) 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 5HT and somatostatin (SOM) neurons. We next developed a singleneuron modeling framework that combines the realism of HodgkinHuxley models with the simplicity and predictive power of generalized integrateandfire models. We found that feedforward inhibition of 5HT neurons by heterogeneous SOM neurons implemented divisive inhibition, while endocannabinoidmediated modulation of excitatory drive to the DRN increased the gain of 5HT 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 5HT neurons, including a previously undescribed dynamic threshold. By applying a bottomup 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 integrateandfire [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 afterhyperpolarization and Atype potassium currents, in combination with heterogeneous feedforward inhibition from local GABA neurons, give rise to a derivativelike inputoutput relationship in serotonin neurons.
https://doi.org/10.7554/eLife.72951.sa0Introduction
The forebrainprojecting serotonin (5HT) 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 5HT neurons into anatomical submodules 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 5HT 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 5HT’s actions. These anatomical and dynamical perspectives on 5HT 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 5HT regulates behavior.
The spiking statistics of 5HT neurons necessarily shape and constrain their computational role. For instance, the slow firing rate (~5 Hz) of 5HT neurons, in large part attributable to a large afterhyperpolarization 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 ensemblerate codes (Knight, 1972; Gerstner, 2000). Consistent with this idea, the in vivo population activity of 5HT neurons has been observed to track impending rewards over second to subsecond timescales (Zhong et al., 2017), and the trialaveraged ensemble rates of individual 5HT neurons can track environmental changes over the millisecond timescale (Ranade and Mainen, 2009; Cohen et al., 2015). In addition, the fact that 5HT receptor subtypes can regulate the excitability of target neurons over different timescales, including ionotropic 5HT3 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 5HT 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 networklevel function emerges from excitability features identified at the singlecell level. In spite of their conceptual utility, the most detailed single cell models, including those of DRN neurons (Tuckwell and Penington, 2014; WongLin et al., 2011), do not lend themselves with ease to bottomup 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 integrateandfire (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 higherorder network function.
In this study, we developed and validated for DRN neurons a hybrid modeling approach that lies between reductionist GIF and biophysical HodgkinHuxleytype 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 5HT 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 bestperforming models of 5HT neurons featured slow membrane time constants, an Atype potassium current, and strong adaptation mechanisms. Network simulations of optimized GIF models of both 5HT and GABAergic SOM neurons organized in a feedforward inhibitory circuit revealed that 5HT neuron populations contextdependently 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 singleneuron models of the two main cell types found in the DRN: 5HT and SOM GABA neurons. We performed wholecell electrophysiological recordings from genetically identified 5HT (Figure 1A1; SERTCre::RosaTdTomato mice) and SOM (Figure 1A2; Table 1; SOMCre::RosaTdTomato mice) neurons in slices. In keeping with previous descriptions (e.g. Vandermaelen and Aghajanian, 1983; Calizo et al., 2011), in the majority of the 5HT 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 5HT 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 5HT neurons and responded to weaker inputs (Figure 1B). The gain showed greater variability in SOM neurons than in 5HT neurons (BrownForsythe equality of variance test p=0.001 on N = 17 5HT 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 5HT neurons fired at 2.81 ± 2.22 Hz vs 8.16 ± 5.70 Hz; BrownForsythe test p=0.005 in $N=17$ 5HT neurons and $N=14$ SOM cells). Together, these observations outlined three salient cellularlevel features of DRN neurons, namely the strong AHP and voltage kink of 5HT neurons as well as noticeable heterogeneous excitability of SOM neurons.
The characteristic kink in the voltage leading up to the first spike in 5HT neurons in principle may be caused by nearthreshold activation of voltagegated potassium channels (VGKCs; Connor and Stevens, 1971; Connor et al., 1977; Drion et al., 2015). We therefore examined wholecell currents evoked by voltage steps (from –90 mV to –20 mV) in both 5HT and SOM neurons to look for evidence of such a VGKC. In 5HT cells, these experiments revealed a large (peak amplitude 928 ± 249 pA, leaksubtracted), partly inactivating (steadystate amplitude 142 ± 45 pA, leaksubtracted) outward current (Figure 1C1) that was sensitive to K_{v}4selective 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 ${\tau}_{h}=42.9\pm 9.4$ ms; kinetics are similar at nearphysiological 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 5HT neurons are broadly similar to those expected of the Atype 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 steadystate component as ${I}_{K}$. Thus, an I_{A} like inactivating VGKC is a consistent feature of DRN 5HT neurons.
The same voltageclamp 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 5HT neurons (Figure 1E3). Together, these results show that the expression of this subthreshold voltagegated current is substantially more variable in SOM neurons than in 5HT 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 5HT populations, we first created a toy leaky integrateandfire (LIF) model that captured the effect of this conductance on singlecell 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 Atype 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 5HT 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 wholecell recordings from identified 5HT 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 4AP ( p=3.2e6 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 initialvoltage/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 5HT cells.
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 timelocked 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 timevarying 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, $\kappa $, 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 spiketriggered current mediating the commonly observed AHP, $\eta $, and change in firing threshold, $\gamma $ (Figure 3A1, see Methods). These components are described by parameters, the values of which are inferred from the electrophysiological data using a combination of leastsquares multilinear 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.
Our results outlined in Figure 2 show that I_{A} regulates spike timing in 5HT 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 HodgkinHuxleytype model of the subthreshold voltagedependent currents we recorded in 5HT 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 channelinactivation GIF model (iGIF; Figure 3A3), which extends the GIF model of Mensi et al., 2012 by adding a nonparametric 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 nonparametric 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 wholecell electrophysiological recordings not only from DRN 5HT and SOM cells but also from canonical deeplayer pyramidal neurons of the medial prefrontal cortex (mPFC). For each recording, we applied two distinct instantiations of noisy in vivolike 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^{2}, 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^{2} = 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^{2} = 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 5HT neurons. As previously intuited, the canonical GIF model performed rather poorly in 5HT neurons (Figure 3F), predicting <15% of the variance of the subthreshold voltage (R^{2} = 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 5HT neurons. By augmenting the GIF model with our experimentally constrained model of I_{A}, the aGIF model not only better predicted the voltage (R^{2} = 0.301 ± 0.200, p=1.96e4; Figure 3F1) but also the spike timing (${{M}_{d}}^{*}$ = 0.481 ± 0.148, p=0.001; Figure 3F2) of 5HT 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.89e4), 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 5HT neurons to in vivolike 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^{2} = 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 5HT neurons
Our model selection approach identified the most salient components required to capture the inputoutput functions of individual neurons and allowed us to identify functional differences across cell types. 5HT 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 5HT neurons produced a substantial and longlasting 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 spiketriggered 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, 5HT neurons are characterized by slow membrane dynamics, I_{A}, and particularly prominent adaptation mechanisms.
Preferential sensitivity of 5HT neuron population to the onset of sudden inputs
The development and validation of accurate singlecell models allowed us to identify the populationlevel computations operating in the DRN. We took advantage of the onetoone correspondence between our GIF models and real neurons to construct synthetic populations with realistic neurontoneuron heterogeneity by sampling from banks of singlecell models (Figure 5A). In response to step increases of synapticlike inputs delivered to the entire population (Figure 5B left), the population firing rates (in Hz/neuron; Figure 5B right) of 5HT, 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 inputoutput functions were approximately rectified linear functions (Figure 5—figure supplement 3) which we summarized and plotted as the timevarying 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 5HT 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 5HT neurons falls to 2.13 ± 0.03 near physiological temperature [Figure 5—figure supplement 5], consistent with a smaller spiketriggered threshold movement [Figure 4—figure supplement 2]). This marked response of 5HT cells occurred quickly, in the first 100 ms after the onset of the step. Thus, despite 5HT neurons being characterized by slow membrane time constants, their population activity provided a remarkably strong encoding of the onset of step synaptic inputs.
We next considered the underlying mechanisms giving rise to the distinctive timedependent gain of 5HT neurons. We found that the characteristically strong spiketriggered adaptation of 5HT neurons (spiketriggered 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 5HT models dramatically reduced the ratio of transient to stationary gain, and viceversa (Figure 5F). These findings are consistent with previous models in other cell types showing that spiketriggered 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 5HT 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 5HT output gain of the DRN
Apart from the strong adaptation mechanisms of 5HT neurons, two other mechanisms have the potential to dynamically modulate the 5HT output from the DRN: I_{A} in 5HT 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 5HT population models using experimentally constrained GABA_{A} receptormediated 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 5HT and SOM neuron populations and examined 5HT neuron population dynamics (as in Figure 5) while varying the maximal conductance of I_{A} (in 5HT neurons). The gain of the transient component of the 5HT 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 5HT output evoked by sustained inputs, with negligible effects on the gain of the slower stationary component.
Previous work has shown that glutamatergic excitatory inputs from the PFC make strong monosynaptic contacts onto both DRN 5HT 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 5HT 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 5HT neuron population dynamics with or without SOM cells (Figure 6B). Including FFI onto 5HT neurons substantially dampened the overall response of the 5HT 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 5HT 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 5HT neurons intact. By favoring the direct monosynaptic excitation of 5HT 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 targetspecific endocannabinoidmediated modulation of PFC excitatory drive in DRN exerts a normalizing role by increasing the overall gain of 5HT 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 inputoutput function (i.e. divisive inhibition), homogeneous FFI mainly shifted the inputoutput 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 nonlinearity in the inputoutput 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.
5HT 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 5HT 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 5HT 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 brainwide target a mixture of the intensity and temporalderivative of its excitatory inputs, and that the derivativeencoding component dominates when the input is increasing rapidly (Figure 7—figure supplement 2).
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 5HT 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 derivativeencoding properties of the DRN are expected to interact. Consistent with this idea, removing I_{A} from 5HT neurons in our DRN network models extended the range of background input where the peak 5HT 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 spikefrequency adaptation in 5HT 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 bottomup approach grounded in experimentally constrained models of the two most abundant cell types in this region: 5HT and SOM GABA neurons. Consistent with, and extending, previous work, we found that 5HT neurons were relatively homogeneous and characterized by potent spikefrequency adaptation (Figure 1) and by the presence of a strong Atype 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 nonlinear subthreshold effects of I_{A} observed in 5HT neurons were required to adequately capture the spiking response of 5HT neurons to naturalistic stimuli (Figure 3). This work introduces a new approach to capturing such nonlinear subthreshold effects in the form of the aGIF model, which augments the GIF model of Mensi et al., 2012 with experimentally constrained HodgkinHuxley style currents, improving model interpretability without compromising predictive performance. Inspecting the parameters of the best performing GIF models revealed that the substantial spikefrequency adaptation observed in 5HT neurons is not fully explained by their distinctively large AHPs and is partly mediated by a previously undescribed dynamic spike threshold (Figure 4). This modelbased approach allowed us to probe causal relationships between specific excitability features and population computations. Thus, we found that the prominent adaptation mechanisms in 5HT 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 5HT neurons linearly reported a mixture of the intensity and temporalderivative of their synaptic inputs (Figure 7), and that the temporalderivative 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 biophysicalsimplified methodology
The computational and statistical modeling methodology presented here was designed to bridge the gap between specific biophysical mechanisms and networklevel computation. Closing this gap has also been the target of complex biophysical simulations, motivated by the hope to create tools for testing diseaserelated 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 5HT neurons, allowing us to probe their contributions to networklevel 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 5HT neurons, it lends itself equally well to capturing the effects of other subthreshold voltagegated 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 voltagedependent 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. spiketriggered 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 linearnonlinear (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 5HT 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 networklevel properties of 5HT signaling. For instance, the role of spiketriggered adaptation in conveying preferential sensitivity to suddenly changing inputs arises in GLMs (Naud and Gerstner, 2012), but the statedependence of the input derivative sensitivity identified in 5HT 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 GLMbased approaches for cell types with very low firing rates or highly statedependent output.
Does the aGIF modeling approach represent an unnecessary complication of the GIF model framework or, conversely, an oversimplification of detailed HodgkinHuxley 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 5HT 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), nonaugmented GIF models remain suitable tools. In our case, it would not have been possible to probe the effect of I_{A} on the networklevel processing features of the DRN without the aGIF model.
Networklevel 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 5HT 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 spiketriggered 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.
5HT neuron heterogeneity
5HT 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 5HT 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 5HT 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 inputoutput function. For FFI, a divisive effect on the gain of stationary inputoutput 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 circuitlevel mechanisms to tune their sensitivity to relevant inputs while maintaining overall stability. This point is germane to 5HT 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 5HT 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 5HT 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, 5HT 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 higherorder 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 cannabinoidmediated preferential reduction of FFI onto 5HT 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 derivativeencoding 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 derivativeencoding 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 5HT signaling in modulating behavior is increasingly conceptualized through the lens of reinforcement learning (RL) theory. Indeed, 5HT 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 stateaction 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 derivativelike computation described here have a place in an RLbased 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 subsecond fluctuations in DRN 5HT 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 RLbased account of DRN function.
If the electrophysiological features of individual 5HT 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 errorlike coding properties (Schultz et al., 1997), bear some electrophysiological features in common with DRN 5HT neurons, with both cell types exhibiting strong adaptation and a prominent Atype 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 derivativelike 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 spiketriggered 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 5HT activity continues, our results provide a new perspective on the fast component of 5HT neuron dynamics: fluctuations in 5HT 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. Slc6a4cre::RosaTdTomato (SERTCre) and Sstcre::RosaTdTomato transgenic lines were used to fluorescently label DRN 5HT and somatostatin (SOM) GABA neurons, respectively. Animals were grouphoused 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 CMM164, CMM176, CMM1711, CMM1743, and CMM2737).
Slice preparation
Request a detailed protocolAnimals were deeply anesthetized using isofluorane (Baxter Corporation) before being euthanized by decapitation. The brain was quickly removed from the skull and submerged into icecold dissection buffer containing the following: 119.0 mM choline chloride, 2.5 mM KCl, 4.3 mM MgSO_{4}, 1.0 mM CaCl_{2}, 1.0 mM NaH_{2}PO_{4}, 1.3 mM sodium ascorbate, 11.0 mM glucose, 26.2 mM NaHCO_{3}; saturated with 95% O_{2}/5% CO_{2}. 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 icecold 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_{4}, 2.5 mM CaCl_{2}, 1.0 mM NaH_{2}PO_{4}, 11.0 mM glucose, 26.2 mM NaHCO_{3}; ∼298 mOsm, maintained at 37°C, and continuously bubbled with 95% O_{2}/5% CO_{2}. The recovery chamber was allowed to equilibrate to room temperature for 1 h before beginning experiments.
In vitro wholecell electrophysiological recording
Request a detailed protocolNeurons were visualized using an upright microscope (Olympus BX51WI) equipped with differential interference contrast and a ×40, 0.8 NA waterimmersion objective. Wholecell recordings were obtained from fluorescently labeled DRN 5HT 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 gluconatebased 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 cesiumbased internal solution (120 mM CsMeSO_{3}, 10 mM EGTA, 5 mM TEA Cl, 1 mM CaCl_{2}, 10 mM Na HEPES, 4 mM Mg ATP, 2 mM GTP, 2 mM QX314, and 10 mM Na phosphocreatine; adjusted to pH 7.25 with CsOH, 280–290 mOsm) and in the presence of bathapplied 100 µM (2 R)amino5phosphonovaleric acid (APV) and 5 µM 2,3dioxo6nitro1,2,3,4tetrahydrobenzo[f]quinoxaline7sulfonamide (NBQX). For voltage clamp experiments, wholecell capacitance compensation was applied manually following breakin, 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 5HT 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 wholecell 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 protocolThe GIF and Nainactivation 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
where $\left\{{\hat{t}}_{i}\right\}$ is the set of spike times and $\eta (t)=\sum _{j}{w}_{j}\mathrm{exp}\left[t/{\tau}_{j}^{(\eta )}\right]$ is the spiketriggered adaptation current. Here the $w}_{j$ are coefficients estimated from the data and the ${\tau}_{j}^{\left(\eta \right)}$ are fixed hyperparameters; see Appendix for details. The GIF emits spikes according to an inhomogeneous Poisson process with intensity $\lambda \left(t\right)$ , given by
where ${V}_{T}^{*}$ is the stationary threshold, $\gamma \left(t\right)=\sum _{j}^{}{\beta}_{j}^{\left(\gamma \right)}\mathrm{exp}\left[t/{\tau}_{j}^{\left(\gamma \right)}\right]$ is the spiketriggered threshold movement (where the ${\beta}_{j}^{\left(\gamma \right)}$ are coefficients estimated from the data and the ${\tau}_{j}^{\left(\gamma \right)}$ are fixed; see Appendix), $\Delta V$ is the threshold sharpness (mV; larger values increase the stochasticity of spiking), and ${\lambda}_{0}=1$ Hz is a constant such that $\lambda \left(t\right)$ is in units of Hz. In the iGIF, an additional variable $\theta \left(t\right)$ is added to the numerator of the exponentiated term in Equation 2 to account for voltagedependent changes in threshold:
The equilibrium voltagedependent change in spike threshold ${\theta}_{\mathrm{\infty}}(V)=\sum _{j=1}^{{N}_{step}}{\beta}_{j}^{(\theta )}\text{rect}\left[V;\text{}{A}_{j},\text{}{A}_{j+1}\right]$ is a piecewise constant function of voltage where each ${{\beta}_{j}}^{\left(\theta \right)}$ defines the value of ${\theta}_{\infty}\left(V\right)$ over the voltage range ${[A}_{j},{A}_{j+1})$ and ${N}_{step}=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 HodgkinHuxley currents which together capture the voltagegated potassium currents found in 5HT neurons (see ‘Potassium current’, below) are added to the subthreshold dynamics given in Equation 1, yielding
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 nonstationary 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 protocolOur toy model of a neuron with an inactivating potassium current is based on an LIF augmented with ${I}_{A}\left(t\right)$ (see ‘Potassium current’ below):
where ${g}_{l}$ and ${E}_{l}$ are the leak conductance and reversal, respectively, and ${I}_{inj}\left(t\right)$ is the external input to the model. To reduce the number of free parameters, the model we used is nondimensionalized with respect to the membrane time constant $\tau}_{mem}=C/{g}_{l$ and leak conductance ${g}_{l}$ , yielding
where $t$ is in units of the membrane time constant, $\overline{g}}_{A}^{\prime}={\overline{g}}_{A}/{g}_{l$ is the effective maximum conductance associated with ${I}_{A}$ , and ${V}_{inj}\left(t\right)={I}_{inj}\left(t\right)/{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 protocolThe voltagegated potassium currents in 5HT neurons were modeled in terms of an inactivating current and a noninactivating current, we refer to as I_{A} and I_{K} , respectively. These were defined as follows
where $\overline{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 $\overline{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
where ${{V}_{{x}^{}}}^{*}$ is the halfactivation 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 deinactivates with a single time constant ${\tau}_{h}$ (ms) that does not depend on voltage. The time dynamics of the $h$ gate are therefore given by
Quantification of singleneuron model performance
Request a detailed protocol${R}^{2}$ was calculated based on the training set $\frac{dV}{dt}$ predicted by the subthreshold component of a given GIF model (Equations 1 and 3, where the spike times $\widehat{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 5HT 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
where ${n}_{dm}$ is the number of modelpredicted spikes that occur within 8 ms of a spike in the validation data, and ${{n}_{dd}}^{*}$ and ${n}_{mm}$ are the corresponding numbers of coincident spikes across sweeps in the validation data and model predictions (where ${{n}_{dd}}^{*}$ is corrected for small sample bias). ${{M}_{d}}^{*}$ can be interpreted as the fraction of modelpredicted 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 protocolDRN network models were constructed by connecting a population of 400 SOM neuron models to a population of 600 5HT neuron models in a feedforward arrangement. Population models were bootstrapped by sampling with replacement from a bank of experimentally constrained GIF models. SOM neuron models were randomly connected to 5HT neuron models with a connection probability of 2%, such that the expected number of GABAergic synapses on each 5HT neuron model was 8. We used a conductancebased model of GABAergic synapses with a fixed reversal potential of –76.7 mV, conductance of 0.3 nS, and biexponential kinetics with ${\tau}_{rise}=1.44$ ms, ${\tau}_{decay}=26.0$ ms, and a propagation delay of 2.0 ms.
Simulated 5HT populations with decreased or increased I_{A} were generated by setting ${\stackrel{\u0304}{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 5HT 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 ${{{\beta}_{j}}^{\left(\gamma \right)}}_{}$ and $w}_{j$. This procedure is summarized in the following pseudocode:
for 5HT_model in 5HT_population; do
SOM_model = random_choice(SOM_models)
5HT_model.eta.coefficients= SOM_model.eta.coefficients
5HT_model.gamma.coefficients= SOM_model.gamma.coefficients
end for.
Numerical methods
Request a detailed protocolSimulations were implemented in Python and C++ using customwritten 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\mathit{}\mathit{\tau}}_{mem}$ 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. Nonparametric tests were used for all twosample comparisons (MannWhitney U test for unpaired samples and Wilcoxon signedrank test for paired samples). Nonparametric 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, pvalues 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 pvalues 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
where $V(t)$ is the membrane voltage, ${I}_{\text{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}\eta (t{\hat{t}}_{i})$ is the adaptation current $\eta $ summated over all past spikes $\{{\hat{t}}_{i}\in \mathcal{S}:{\hat{t}}_{i}<t\}$ ($\{\mathcal{S}\}$ is the set of all spike times).
The adaptation current produced by a single spike $\eta (t{\widehat{t}}_{i})$ is implemented as a sum of $k$ exponentials given by
where the timescales ${\tau}_{1},{\tau}_{2},\mathrm{\dots},{\tau}_{k}$ are treated as hyperparameters. If we substituted this implementation of $\eta $ back into Equation A1, the term associated with the spiketriggered current $\sum _{{\hat{t}}_{i}<t}\eta (t{\hat{t}}_{i})$ would become $\sum _{{\hat{t}}_{i}<t}\sum _{j=1}^{k}{w}_{j}{e}^{(t{\hat{t}}_{i})/{\tau}_{j}}$. This double sum can be written more concisely as
where $\hat{\eta}}_{j}(t)=\sum _{{\hat{t}}_{i}<t}{e}^{(t{\hat{t}}_{i})/{\tau}_{j}$ is a basis for the adaptation current over the timescale ${\tau}_{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:
Estimating aGIF model subthreshold parameters
Given a training dataset $\mathcal{D}=\{(V(t),{I}_{\text{inj}}(t)):1\le t\le T\}$, knowledge of the equilibrium gating functions ${m}_{\mathrm{\infty}},{h}_{\mathrm{\infty}},\text{and}{n}_{\mathrm{\infty}}$, and appropriate choices of ${\tau}_{1},{\tau}_{2},\mathrm{\dots},{\tau}_{k}$ in $\eta $, our goal is to estimate the remaining parameters in Equation A4; namely, ${g}_{l},C,{E}_{l},{\overline{g}}_{A},{\overline{g}}_{K},{w}_{1},\mathrm{\dots},{w}_{k},\text{and}{\tau}_{h}$, where ${\tau}_{h}$ is the time constant of the inactivation gate $h$. Fortunately, all of these except for ${\tau}_{h}$ can be estimated easily using linear regression.
We begin by rewriting Equation A4 as the product of a row vector of predictors $\mathbf{x}$ and a column vector of coefficients $\beta$ as follows
We solve this subject to ${g}_{l},C,{\overline{g}}_{A},{\overline{g}}_{K}\ge 0$ using scipy.optimize.lsq_linear.
Next we turn to the question of calculating all of the components of $\mathbf{x}$. Because Equation A4 only reflects the subthreshold dynamics of the aGIF model, we begin by removing all time points in $\mathcal{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}_{\mathrm{\infty}},{h}_{\mathrm{\infty}},{n}_{\mathrm{\infty}},{D}_{K}\text{and}{\widehat{\eta}}_{i}(t)$. To calculate $h$ for each time point, we order the values of ${h}_{\mathrm{\infty}}$ according to time and integrate $h$ numerically using a fixed time step, and the initial condition $h={h}_{\mathrm{\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 nonconvex function of ${\tau}_{h}$. We therefore conducted a line search over plausible values of ${\tau}_{h}$ and chose the value associated with the highest variance explained. This is equivalent to solving
where $\theta =(\beta ,{\tau}_{h})$.
Singleneuron model hyperparameters. For more details on the iGIF model hyperparameter , see Mensi et al., 2016.
Model  Parameter  Symbol  Cell type  Value (ms) 

All  $\eta $ timescales  ${\tau}_{1},{\tau}_{2},\mathrm{\dots},{\tau}_{k}$  All  3, 10, 30, 100, 300, 1000, 3000 
All  $\gamma $ timescales  None  All  3, 30, 300, 3000 
All  Refractory period  None  5HT  6.5 
SOM and mPFC  4.0  
iGIF  Candidate thresholdcoupling timescales  ${\tau}_{\theta}$  All  1, 2, 5, 10, 22, 46, 100 
aGIF  Candidate inactivation timescales  ${\tau}_{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).

Dryad Digital RepositoryPatchclamp recordings from dorsal raphe neurons.https://doi.org/10.5061/dryad.66t1g1k2w
References

Anatomic and functional topography of the dorsal raphe nucleusAnnals of the New York Academy of Sciences 1018:46–57.https://doi.org/10.1196/annals.1296.005

D2like dopamine receptors depolarize dorsal raphe serotonin neurons through the activation of nonselective cationic conductanceThe Journal of Pharmacology and Experimental Therapeutics 320:376–385.https://doi.org/10.1124/jpet.106.111690

The exploitation of regularities in the environment by the brainThe Behavioral and Brain Sciences 24:602–607.https://doi.org/10.1017/s0140525x01000024

A universal model for spikefrequency adaptationNeural Computation 15:2523–2564.https://doi.org/10.1162/089976603322385063

Opponency revisited: competition and cooperation between dopamine and serotoninNeuropsychopharmacology 36:74–97.https://doi.org/10.1038/npp.2010.151

Two major network domains in the dorsal raphe nucleusThe Journal of Comparative Neurology 523:1488–1504.https://doi.org/10.1002/cne.23748

Prediction of repetitive firing behaviour from voltage clamp data on an isolated neurone somaThe Journal of Physiology 213:31–53.https://doi.org/10.1113/jphysiol.1971.sp009366

Serotonin and dopamine: unifying affective, activational, and decision functionsNeuropsychopharmacology 36:98–113.https://doi.org/10.1038/npp.2010.121

Opponent interactions between serotonin and dopamineNeural Networks 15:603–616.https://doi.org/10.1016/s08936080(02)000527

Serotonin in affective controlAnnual Review of Neuroscience 32:95–126.https://doi.org/10.1146/annurev.neuro.051508.135607

Serotonin in panic, anxiety and depressionEuropean Neuropsychopharmacology 1:358–360.https://doi.org/10.1016/0924977X(91)90566D

Serotonin and the neurobiology of depressionArchives of General Psychiatry 51:865.https://doi.org/10.1001/archpsyc.1994.03950110025005

Genetic approaches for understanding the role of serotonin receptors in mood and behaviorCurrent Opinion in Neurobiology 23:399–406.https://doi.org/10.1016/j.conb.2013.01.011

Metalearning and neuromodulationNeural Networks 15:495–506.https://doi.org/10.1016/s08936080(02)000448

Linearization of FI curves by adaptationNeural Computation 10:1721–1729.https://doi.org/10.1162/089976698300017106

Mechanisms underlying gain modulation in the cortexNature Reviews. Neuroscience 21:80–92.https://doi.org/10.1038/s415830190253y

Serotonin neurons modulate learning rate through uncertaintyCurrent Biology 32:586–599.https://doi.org/10.1016/j.cub.2021.12.006

Efficient estimation of detailed singleneuron modelsJournal of Neurophysiology 96:872–890.https://doi.org/10.1152/jn.00079.2006

Smoothing of, and parameter estimation from, noisy biophysical recordingsPLOS Computational Biology 5:e1000379.https://doi.org/10.1371/journal.pcbi.1000379

A benchmark test for A quantitative assessment of simple neuron modelsJournal of Neuroscience Methods 169:417–424.https://doi.org/10.1016/j.jneumeth.2007.11.006

Dynamic, nonlinear feedback regulation of slow pacemaking by Atype potassium current in ventral tegmental area neuronsThe Journal of Neuroscience 28:10905–10917.https://doi.org/10.1523/JNEUROSCI.223708.2008

Dynamics of encoding in a population of neuronsThe Journal of General Physiology 59:734–766.https://doi.org/10.1085/jgp.59.6.734

Visual adaptation: physiology, mechanisms, and functional benefitsJournal of Neurophysiology 97:3155–3164.https://doi.org/10.1152/jn.00086.2007

Serotonin neurons in the dorsal raphe nucleus encode reward signalsNature Communications 7:10503.https://doi.org/10.1038/ncomms10503

Fractional differentiation by neocortical pyramidal neuronsNature Neuroscience 11:1335–1342.https://doi.org/10.1038/nn.2212

Do dorsal raphe 5HT neurons encode “beneficialness”?Neurobiology of Learning and Memory 135:40–49.https://doi.org/10.1016/j.nlm.2016.08.008

Stressor controllability and learned helplessness: the roles of the dorsal raphe nucleus, serotonin, and corticotropinreleasing factorNeuroscience and Biobehavioral Reviews 29:829–841.https://doi.org/10.1016/j.neubiorev.2005.03.021

Optimal heterogeneity for coding in spiking neural networksPhysical Review Letters 108:228102.https://doi.org/10.1103/PhysRevLett.108.228102

Differential effects of excitatory and inhibitory heterogeneity on the gain and asynchronous state of sparse cortical networksFrontiers in Computational Neuroscience 8:107.https://doi.org/10.3389/fncom.2014.00107

Subtractive, divisive and nonmonotonic gain control in feedforward nets linearized by noise and delaysFrontiers in Computational Neuroscience 8:19.https://doi.org/10.3389/fncom.2014.00019

Parameter extraction and classification of three cortical neuron types reveals two distinct adaptation mechanismsJournal of Neurophysiology 107:1756–1775.https://doi.org/10.1152/jn.00408.2011

Rewarddependent modulation of neuronal activity in the primate dorsal raphe nucleusThe Journal of Neuroscience 28:5331–5343.https://doi.org/10.1523/JNEUROSCI.002108.2008

Improved similarity measures for small sets of spike trainsNeural Computation 23:3016–3069.https://doi.org/10.1162/NECO_a_00208

Coding and decoding with adapting neurons: A population approach to the peristimulus time histogramPLOS Computational Biology 8:e1002711.https://doi.org/10.1371/journal.pcbi.1002711

Embracing diversity in the 5HT neuronal systemNature Reviews. Neuroscience 20:397–424.https://doi.org/10.1038/s4158301901513

Classes of dendritic information processingCurrent Opinion in Neurobiology 58:78–85.https://doi.org/10.1016/j.conb.2019.07.006

Temporal whitening by powerlaw adaptation in neocortical neuronsNature Neuroscience 16:942–948.https://doi.org/10.1038/nn.3431

Automated highthroughput characterization of single neurons by means of simplified spiking modelsPLOS Computational Biology 11:e1004275.https://doi.org/10.1371/journal.pcbi.1004275

Similar network activity from disparate circuit parametersNature Neuroscience 7:1345–1352.https://doi.org/10.1038/nn1352

Transient firing of dorsal raphe neurons encodes diverse and specific sensory, motor, and reward eventsJournal of Neurophysiology 102:3026–3037.https://doi.org/10.1152/jn.00507.2009

5ht(1a) receptor function in major depressive disorderProgress in Neurobiology 88:17–31.https://doi.org/10.1016/j.pneurobio.2009.01.009

Lowserotonin levels increase delayed reward discounting in humansThe Journal of Neuroscience 28:4528–4532.https://doi.org/10.1523/JNEUROSCI.498207.2008

Reconciling the role of central serotonin neurons in human and animal behaviorBehavioral and Brain Sciences 9:319–335.https://doi.org/10.1017/S0140525X00022871

An afterhyperpolarization of medium duration in rat hippocampal pyramidal cellsThe Journal of Physiology 409:171–190.https://doi.org/10.1113/jphysiol.1989.sp017491

Serotonin: modulator of a drive to withdrawBrain and Cognition 71:427–436.https://doi.org/10.1016/j.bandc.2009.03.009

Processing of lowprobability sounds by cortical neuronsNature Neuroscience 6:391–398.https://doi.org/10.1038/nn1032

ConferenceA Spiking Neuronal Network Model of the Dorsal Raphe Nucleus2011 International Joint Conference on Neural Networks.https://doi.org/10.1109/IJCNN.2011.6033414

Tryptophan depletion causes a rapid lowering of mood in normal malesPsychopharmacology 87:173–177.https://doi.org/10.1007/BF00431803

Learning and stress shape the reward response patterns of serotonin neuronsThe Journal of Neuroscience 37:8863–8875.https://doi.org/10.1523/JNEUROSCI.118117.2017
Decision letter

Naoshige UchidaReviewing Editor; Harvard University, United States

John R HuguenardSenior Editor; Stanford University School of Medicine, United States

Jochen RoeperReviewer; Goethe University Frankfurt, Germany

Hitoshi MorikawaReviewer; University of Texas at Austin, United States

Paul MillerReviewer; 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 experimentallydriven augmented integrateandfire 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 wholecell patch clamp recording in brain slices and biophysical computational modeling. The authors adopted the approach called a generalized integrateandfire [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 afterhyperpolarization potentials. The second is a "kink" in the voltage trace leading up to the first spike, which can be attributed to the effect of Atype potassium currents (IA). Using the data, the authors constructed aGIF model of serotonin neurons which contains IA, two HodgkinHuxley 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 inputoutput 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 spikefrequency 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 Atype 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 Atype currents including the inactivation kinetics?
3. Some of the model parameters appear to be nonphysiological. 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 derivativelike computation holds in all of potential input regimes. It is important to clarify in what conditions derivativelike computation is a good approximation of the inputoutput relationship of serotonin neurons, and when it breaks down.
5. The derivativelike computation and its role in reinforcement learning are emphasized in the manuscript (e.g. abstract). However, as discussed above, the derivativelike 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 5HT 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.
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 IAcurrents will be distorted under these recording conditions
3. An intracellular solution is used without calcium buffering. Given the key role of adaptation and relevant calciumdependent conductances in 5HT neurons this is not stateoftheart
4. 4AP is not a selective Kv4 (IA) channel blocker. There are moreselective, commercially available alternatives like Heteropodatoxin2 or AmmTx3.
5. Author show that augmentation of the GIFmodel with IAcurrents 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 5HT neurons (see Tuckwell & Penington (2014) Computational modeling of spike generation in serotonergic neurons of the dorsal raphe nucleus. Prog Neurobiol 118:59101).
6. Authors should experimentally study the biophysics of adaptation in 5HT 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 Atype 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.21p.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, wellcharacterized 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 5HT 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 "inputoutput 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 xy 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 timederivative 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 trialaveraged 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.sa1Author response
Essential revisions:
In this study, Harkin and colleagues examined the biophysical properties of dorsal raphe (DR) serotonin neurons using a combination of wholecell patch clamp recording in brain slices and biophysical computational modeling. The authors adopted the approach called a generalized integrateandfire [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 afterhyperpolarization potentials. The second is a "kink" in the voltage trace leading up to the first spike, which can be attributed to the effect of Atype potassium currents (IA). Using the data, the authors constructed aGIF model of serotonin neurons which contains IA, two HodgkinHuxley 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 inputoutput 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 spikefrequency 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 voltagedependence 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 spiketriggered 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 Atype 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 Atype 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 voltageclamp recordings used to characterize IA in 5HT 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.
Conversely, how robust are the modeling results with respect to the specific gating properties of Atype 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 Rarelated 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 resistancerelated 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 5HT 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 resistancemediated distortion.
3. Some of the model parameters appear to be nonphysiological. 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 derivativelike computation holds in all of potential input regimes. It is important to clarify in what conditions derivativelike computation is a good approximation of the inputoutput relationship of serotonin neurons, and when it breaks down.
This is an important point, especially considering that 5HT 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 temporalderivative of its input. We have added a toy ratebased model of the DRN to capture this idea and to clarify its limitations (Figure 7S2). Reassuringly, we found that the derivativelike 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 derivativelike computation and its role in reinforcement learning are emphasized in the manuscript (e.g. abstract). However, as discussed above, the derivativelike 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 5HT 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 IAcurrents 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 calciumdependent conductances in 5HT neurons this is not stateoftheart
This is an important point. 5HT neurons are indeed endowed with several Calciumdependent conductances that are bound to regulate their firing behaviors. As a mere example, 5HT neurons are well known to exhibit a prominent calcium and SK channel dependent AHP (ScuvéeMoreau et al., Br. J. Pharm, 2004; Sargin et al., eLife 2016). There are likely other parameters of ionic conductances expressed in 5HT 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 (calciumdependent) conductances to these spiking behaviors, these experiments would be fundamentally challenging to interpret.
4. 4AP is not a selective Kv4 (IA) channel blocker. There are moreselective, commercially available alternatives like Heteropodatoxin2 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 Heteropodatoxin2 and AmmTx3 block most of the transient component of the outward current readily observed in 5HT 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 GIFmodel with IAcurrents 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 5HT neurons (see Tuckwell & Penington (2014) Computational modeling of spike generation in serotonergic neurons of the dorsal raphe nucleus. Prog Neurobiol 118:59101).
Thank you for the opportunity to clarify these points. Concerning aGIF vs iGIF, although both models predict the spike trains of 5HT 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 LNLbased approaches discussed on page 31, which are not viable for 5HT 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 HodgkinHuxelytype modelling.
There remains, however, an important drawback to biophysicallydetailed HodgkinHuxleystyle models, such as the model of Tuckwell and Penington (2014). Offtheshelf 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 biophysicallydetailed models also uses about 1000 times more computer time. As a result, even if GIFmodelbased 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 5HT 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 5HT neurons directly from their spiking and subthreshold responses to in vivolike 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 GIFtype models to estimate adaptation is wellestablished (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 voltagedependence (as briefly noted in our discussion on page 30). In an independent study, we have carried out voltageclamp 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 experimentallydetermined 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 modelbased 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., HajDahmane, PNAS 2018) to affect synaptic transmission (e.g., HajDahmane 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 cannabinoidmediated inhibition than their counterpart onto 5HT 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 5HT 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 Atype 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 integrateandfire models including GIFclass 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.21p.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, wellcharacterized 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 “nonlinear coupling term” that differentiates the GIF and iGIF models mentioned in our results text actually consists of a nonparametric function of past voltage. The nonparametric 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 nonparametric 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 spiketiming; the aGIF better predicts the subthreshold dynamics of 5HT 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 5HT 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). Crossvalidation remains a more accurate approach to model selection that takes into account overfitting.
Figure 5: Given the 5HT 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 "inputoutput 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 derivativeencoding 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 xy 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 5HT neurons contribute to shaping the inputoutput 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 timederivative 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 trialaveraged 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.sa2Article and author information
Author details
Funding
Canadian Institutes of Health Research (175325)
 Richard Naud
 JeanClaude Béïque
Canadian Institutes of Health Research (175319)
 Richard Naud
 JeanClaude Béïque
Natural Sciences and Engineering Research Council of Canada (06972)
 Richard Naud
Canada Foundation for Innovation, Brain Canada (Canadian Neurophotonic Platform)
 JeanClaude Béïque
Krembil Foundation
 JeanClaude 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 SOMCre::RosaTdTomato 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 pilottesting of the GIF model in 5HT neurons and all other members of the JCB. 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 CMM164, CMM176, CMM1711, CMM1743, and CMM2737). 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
Publication 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
© 2023, Harkin et al.
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.
Metrics

 817
 Page views

 104
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
Different strains of a microorganism growing in the same environment display a wide variety of growth rates and growth yields. We developed a coarsegrained model to test the hypothesis that different resource allocation strategies, corresponding to different compositions of the proteome, can account for the observed rateyield variability. The model predictions were verified by means of a database of hundreds of published rateyield and uptakesecretion phenotypes of Escherichia coli strains grown in standard laboratory conditions. We found a very good quantitative agreement between the range of predicted and observed growth rates, growth yields, and glucose uptake and acetate secretion rates. These results support the hypothesis that resource allocation is a major explanatory factor of the observed variability of growth rates and growth yields across different bacterial strains. An interesting prediction of our model, supported by the experimental data, is that high growth rates are not necessarily accompanied by low growth yields. The resource allocation strategies enabling highrate, highyield growth of E. coli lead to a higher saturation of enzymes and ribosomes, and thus to a more efficient utilization of proteomic resources. Our model thus contributes to a fundamental understanding of the quantitative relationship between rate and yield in E. coli and other microorganisms. It may also be useful for the rapid screening of strains in metabolic engineering and synthetic biology.

 Computational and Systems Biology
 Neuroscience
Biological motor control is versatile, efficient, and depends on proprioceptive feedback. Muscles are flexible and undergo continuous changes, requiring distributed adaptive control mechanisms that continuously account for the body's state. The canonical role of proprioception is representing the body state. We hypothesize that the proprioceptive system could also be critical for highlevel tasks such as action recognition. To test this theory, we pursued a taskdriven modeling approach, which allowed us to isolate the study of proprioception. We generated a large synthetic dataset of human arm trajectories tracing characters of the Latin alphabet in 3D space, together with muscle activities obtained from a musculoskeletal model and modelbased muscle spindle activity. Next, we compared two classes of tasks: trajectory decoding and action recognition, which allowed us to train hierarchical models to decode either the position and velocity of the endeffector of one's posture or the character (action) identity from the spindle firing patterns. We found that artificial neural networks could robustly solve both tasks, and the networks'units show tuning properties similar to neurons in the primate somatosensory cortex and the brainstem. Remarkably, we found uniformly distributed directional selective units only with the actionrecognitiontrained models and not the trajectorydecodingtrained models. This suggests that proprioceptive encoding is additionally associated with higherlevel functions such as action recognition and therefore provides new, experimentally testable hypotheses of how proprioception aids in adaptive motor control.