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
Article 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.
Version history
 Preprint posted: June 27, 2021 (view preprint)
 Received: August 10, 2021
 Accepted: December 19, 2022
 Accepted Manuscript published: January 19, 2023 (version 1)
 Version of Record published: February 27, 2023 (version 2)
Copyright
© 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

 1,246
 views

 159
 downloads

 4
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
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

 Cell Biology
 Computational and Systems Biology
Bats have unique characteristics compared to other mammals, including increased longevity and higher resistance to cancer and infectious disease. While previous studies have analyzed the metabolic requirements for flight, it is still unclear how bat metabolism supports these unique features, and no study has integrated metabolomics, transcriptomics, and proteomics to characterize bat metabolism. In this work, we performed a multiomics data analysis using a computational model of metabolic fluxes to identify fundamental differences in central metabolism between primary lung fibroblast cell lines from the black flying fox fruit bat (Pteropus alecto) and human. Bat cells showed higher expression levels of Complex I components of electron transport chain (ETC), but, remarkably, a lower rate of oxygen consumption. Computational modeling interpreted these results as indicating that Complex II activity may be low or reversed, similar to an ischemic state. An ischemiclike state of bats was also supported by decreased levels of central metabolites and increased ratios of succinate to fumarate in bat cells. Ischemic states tend to produce reactive oxygen species (ROS), which would be incompatible with the longevity of bats. However, bat cells had higher antioxidant reservoirs (higher total glutathione and higher ratio of NADPH to NADP) despite higher mitochondrial ROS levels. In addition, bat cells were more resistant to glucose deprivation and had increased resistance to ferroptosis, one of the characteristics of which is oxidative stress. Thus, our studies revealed distinct differences in the ETC regulation and metabolic stress responses between human and bat cells.

 Computational and Systems Biology
 Neuroscience
Normal aging leads to myelin alterations in the rhesus monkey dorsolateral prefrontal cortex (dlPFC), which are positively correlated with degree of cognitive impairment. It is hypothesized that remyelination with shorter and thinner myelin sheaths partially compensates for myelin degradation, but computational modeling has not yet explored these two phenomena together systematically. Here, we used a twopronged modeling approach to determine how agerelated myelin changes affect a core cognitive function: spatial working memory. First, we built a multicompartment pyramidal neuron model fit to monkey dlPFC empirical data, with an axon including myelinated segments having paranodes, juxtaparanodes, internodes, and tight junctions. This model was used to quantify conduction velocity (CV) changes and action potential (AP) failures after demyelination and subsequent remyelination. Next, we incorporated the single neuron results into a spiking neural network model of working memory. While complete remyelination nearly recovered axonal transmission and network function to unperturbed levels, our models predict that biologically plausible levels of myelin dystrophy, if uncompensated by other factors, can account for substantial working memory impairment with aging. The present computational study unites empirical data from ultrastructure up to behavior during normal aging, and has broader implications for many demyelinating conditions, such as multiple sclerosis or schizophrenia.