Degeneracy and stability in neural circuits of dopamine and serotonin neuromodulators: A theoretical consideration

Degenerate neural circuits perform the same function despite being structurally different. However, it is unclear whether neural circuits with interacting neuromodulator sources can themselves degenerate while maintaining the same neuromodulatory function. Here, we address this by computationally modeling the neural circuits of neuromodulators serotonin and dopamine, local glutamatergic and GABAergic interneurons, and their possible interactions, under reward/punishment-based conditioning tasks. The neural modeling is constrained by relevant experimental studies of the VTA or DRN system using, e.g., electrophysiology, optogenetics, and voltammetry. We first show that a single parsimonious, sparsely connected neural circuit model can recapitulate several separate experimental findings that indicated diverse, heterogeneous, distributed, and mixed DRNVTA neuronal signaling in reward and punishment tasks. The inability of this model to recapitulate all observed neuronal signaling suggests potentially multiple circuits acting in parallel. Then using computational simulations and dynamical systems analysis, we demonstrate that several different stable circuit architectures can produce the same observed network activity profile, hence demonstrating degeneracy. Due to the extensive D2-mediated connections in the investigated circuits, we simulate the D2 receptor agonist by increasing the connection strengths emanating from the VTA DA neurons. We found that the simulated D2 agonist can distinguish among sub-groups of the degenerate neural circuits based on substantial deviations in specific neural populations’ activities in reward and punishment conditions. This forms a testable model prediction using pharmacological means. Overall, this theoretical work suggests the plausibility of degeneracy within neuromodulator circuitry and has important implications for the stable and robust maintenance of neuromodulatory functions.

1. Introduction

The nervous system can be modulated by endogenous chemical messengers, called neuromodulators (Levitan, 1987). Neurons or synapses, and hence neural circuits, that succumb to neuromodulation can change circuit configuration and function (Marder, 2012; Marder et al., 2014). It is also known that neural circuits can degenerate, that is, circuits can comprise different elements and/or structure while performing the same function or yielding the same output (Cropper et al., 2016). This can allow robust maintenance of functions and behavior in the face of changes in the underlying structure (Edelman and Gally, 2001; Whitacre, 2010). It should be noted that the degeneracy defined here should be distinguished from other definitions e.g., having different states with the same energy level in quantum physics or a decline in health condition (Griffiths and Schroeter, 2018; Hocking et al., 2019). Although it has been shown that neuromodulators can selectively regulate degenerate neural circuits (Marder et al., 2014; Cropper et al., 2016), it is unclear whether neural circuits with neuromodulator-containing neurons can themselves be degenerate, which could in turn provide stable widespread neuromodulator influences on targeted neural circuits (Figure 1A).

FIGURE 1

Figure 1. Degenerate neuromodulator circuits are constrained by observed neuronal signaling. (A) Multiple neuromodulators that influence neural circuits, cognition, and behavior, may be embedded within degenerate neural circuits. Neuromod: Specific neuromodulator type. (B) Schematic of DRN and VTA activity profiles in reward and punishment tasks. Activities (firing rates) aligned to the timing of unexpected punishment outcome (left, vertical red dashed lines) and learned reward-predictive cue (right, vertical green dashed lines) and reward outcome (right, vertical red dashed lines). Top-to-bottom: VTA DA neural activity exhibits phasic excitation (inhibition) at reward-predictive cue (punishment) outcome (e.g., Cohen and Uchida, 2012; Tan et al., 2012). VTA GABAergic neural activity shows phasic excitation upon punishment (e.g., Tan et al., 2012; Eshel et al., 2015) while exhibiting post-cue tonic activity which is not modulated by the presence/absence of actual outcome (e.g., Cohen and Uchida, 2012). DRN Type-I 5-HT neurons show phasic activation by a reward-predicting cue (right) but not punishment (left). DRN Type-II 5-HT neurons signal punishment outcome (left) and sustained activity toward expected reward outcome (right) (e.g., Cohen et al., 2015). DRN GABAergic neurons have phasic activation upon punishment but have tonic inhibition during waiting and reward delivery (e.g., Li et al., 2016). DRN glutamatergic neurons were deduced to be excited by reward-predicting cues, in line with VTA DA neural activation (McDevitt et al., 2014), and assumed not to respond to punishment outcome. Baseline activity for Type-I 5-HT DRN neurons is higher in reward than punishment tasks (e.g., Cohen et al., 2015).

In this theoretical study, we investigate the plausibility of degenerate and stable neuromodulator circuits by focusing on the neural circuits in the midbrain which are the source of ascending pathways of two highly studied monoaminergic neuromodulators, serotonin (5-hydroxytryptamine; 5-HT) and dopamine (DA). These neuromodulators have major roles in modulating cognition, emotion, and behavior, and are linked to the pathogenesis and pharmacological treatment of many common neuropsychiatric and neurological disorders (Muller and Cunningham, 2020). The majority of 5-HT-producing neurons are found in the dorsal and median raphe nuclei (DRN and MRN), while most DA-producing neurons reside in the ventral tegmental area (VTA) and substantia nigra compacta (SNc) (Muller and Cunningham, 2020).

At the functional level, DA and 5-HT are known to play important roles in reward- and punishment-based learning and decision-making (Hu, 2016). For example, there is strong evidence that the DA neuronal activity signals reward prediction error (difference between predicted and actual reward outcome) to guide reinforcement learning (Doya, 2002). Specifically, DA neurons are phasically excited upon unexpected reward outcomes or reward-predictive cues and inhibited upon unexpected reward omission or punishment (Schultz et al., 1997; Cohen and Uchida, 2012; Watabe-Uchida et al., 2017), although there is heterogeneity among DA neurons in this regard (Cohen and Uchida, 2012; Cohen et al., 2015).

In comparison to DA neurons, DRN 5-HT neurons exhibit greater complexity in function, with studies reporting that 5-HT neuronal activity encodes both reward and punishment. For instance, it has been shown that certain 5-HT neurons (labeled “Type I”) were phasically activated only by reward-predicting cues, but not punishment in a classical conditioning paradigm (Cohen et al., 2015). Yet, in the same study, a different population of 5-HT neurons (“Type II”) signaled both expected reward and punishment with sustained elevated activity toward reward outcome (Cohen et al., 2015). The latter study also found that, unlike DA neurons, baseline firing of Type-I 5-HT neurons was generally higher in rewarding than punishment trials, and this effect lasted across many trials, suggesting information processing over a long timescale, possibly used to perform meta-learning via learning rate modulation (Grossman et al., 2022). Similar 5-HT neuronal responses to reward and punishment were reported in other rodents (Liu et al., 2014; Li et al., 2016; Matias et al., 2017; Zhong et al., 2017) and non-human primates (Hayashi et al., 2015) studies. These differential responses of DA and 5-HT neurons to reward and punishment cannot be trivially reconciled with a simple model, e.g., using two opposing neuromodulatory systems as previously proposed (Boureau and Dayan, 2010).

Other studies reveal further complexity in reward/punishment processing, specifically in the form of altered activity of non-5-HT/DA midbrain neurons. For example, DRN neurons utilizing gamma-aminobutyric (GABA) were tonically inhibited during reward waiting with further inhibition during reward acquisition but phasically activated by aversive stimuli (Li et al., 2016). On the contrary, GABAergic neuronal activity in the VTA exhibited sustained activity upon rewarding cue onset but no response to the presence or absence of an actual reward outcome (Cohen and Uchida, 2012). Furthermore, other studies found that VTA GABAergic neuronal activity was potently and phasically activated by punishment outcome, which in turn inhibited the VTA DA neuronal activity (Tan et al., 2012; Eshel et al., 2015). Another study showed that glutamatergic (Glu) neurons in the DRN reinforced instrumental behavior through VTA DA neurons (McDevitt et al., 2014).

This complexity of signaling within the DRN-VTA system in response to reward and punishment may reflect the DRN and VTA having shared afferent inputs (Watabe-Uchida et al., 2012, 2017; Ogawa et al., 2014; Pollak Dorocic et al., 2014; Beier et al., 2015; Tian et al., 2016; Ogawa and Watabe-Uchida, 2018). Another possibility is that the DRN and VTA interact with each other. Indeed, a growing number of studies have suggested that there are direct and indirect interactions among distinctive neuronal types between and within the DRN and VTA (di Giovanni et al., 2008; Tan et al., 2012; Watabe-Uchida et al., 2012; McDevitt et al., 2014; Ogawa et al., 2014; Beier et al., 2015; de Deurwaerdère and di Giovanni, 2017; Valencia-Torres et al., 2017; Xu et al., 2017; Li et al., 2019; Wang et al., 2019).

Taken together, information on reward and punishment signaled by neuronal activities within the DRN and VTA seems to be diverse, heterogeneous, distributed, and mixed. Some of these signaling responses are illustrated in Figure 1B and summarized in Supplementary Table 1. This led to the following questions (Figure 1A). First, can these experimental findings from separate studies be reconciled and understood in terms of a single, parsimonious neural circuit model encompassing both the DRN and VTA? Second, can there be several degenerate DRN-VTA circuits, which are stable, that can signal reward and punishment information as observed in experiments? Third, if there are degenerate and stable DRN-VTA circuits, can some of them be distinguishable from others, for example, through perturbative means?

To address these questions, we developed a biologically plausible DRN-VTA computational neural circuit model based on our previous multiscale modeling framework for neuromodulators (Joshi et al., 2017; Wong-Lin et al., 2017), constrained by some of the observed neuronal signaling profiles under reward and punishment tasks. The modeling considers known direct and indirect pathways between DRN 5-HT and VTA DA neurons. Upon simulating the model under reward and punishment conditions, we found that many, but not all of the experimental findings, could be captured in a single, parsimonious DRN-VTA model. Furthermore, several distinct model architectures could replicate the same neural circuit activity response profile, hence reflecting the possibility of degeneracy. Applying dynamical systems theory, we found that all these circuits were dynamically stable. To distinguish among these degenerated models, we simulated the drug effects of DA D2-receptor-based agonists and were able to distinguish sub-groups of these degenerate model architectures. Overall, our study provides theoretical support for the degeneracy and stability of neural circuits of neuromodulators.

2. Materials and methods

2.1. DRN-VTA network modeling

To develop the DRN-VTA neural circuit models, we made use of our dynamic mean-field (neuronal population-based) modeling framework (Joshi et al., 2017) for neuromodulator interactions. The modeling approach was constrained by data from known electrophysiological, neuropharmacological, and voltammetry parameters [see (Joshi et al., 2017) and below]. Each neural circuit model architecture investigated consisted of DRN 5-HT, VTA DA, VTA GABA, DRN GABA, and DRN Glu neuronal populations. Direct and indirect interactions among these five neuronal populations were then explored. The main aim of this work was to evaluate the plausibility of neuromodulator circuit degeneracy and stability rather than replicate every neuronal population in these brain regions.

Directed interaction between two neuronal populations, in which the source was a neuromodulator, was mathematically described by the neuromodulator neuronal population firing (rate) activity, followed by the release-and-uptake dynamics of the neuromodulator, which in turn induced certain population-averaged currents on the targeted neuronal population (Joshi et al., 2017). An induced current, I_x, which could be effectively excitatory or inhibitory depending on experimental findings, can be described by:

$${\mathrm{\tau }}_x\frac{d{I}_x}{dt}=-I_x+\frac{k_x}{1+e^{-g_x(\left[x\right]-{[x]}_0)}}(1)$$

where x was some neuromodulator (5-HT or DA), τ_x the associated time constant, k_x some constant that determined the current amplitude, and constants g_x and [x]₀ that controlled the slope and offset of the function on the right-hand-side of Eq. (1). The release-and-uptake dynamics for a neuromodulator x are described by:

$$\frac{d[x]}{dt}={[x]}_p{R}_x-\frac{V_{max,x}[x]}{K_{m,x}+[x]}(2)$$

where [x] was the concentration of x, [x]_p the release per neural firing frequency (Joshi et al., 2017), and constants V_max,x and K_m,x were constants determined from voltammetry measurements.

If the source of the interaction came from GABAergic (or Glu) neurons, then for simplicity, we assumed an instantaneous inhibitory (or excitatory) current-based synaptic influence on the targeted neuronal populations. This reflected the faster ionotropic receptor-based synaptic dynamics compared to metabotropic receptor-based neuromodulatory effects, while also reducing the number of free model parameters. Similarly, we also ignored the relatively faster neuronal membrane dynamics. Threshold-linear input–output function for each neuronal population was used (Joshi et al., 2017) and described by:

$$F=g{[I-I_0]}_{+}(3)$$

where F was the neural population firing rate (output), I was the total input current into a neural population, g was the slope, and I₀ was some threshold current, and with [z]₊ = z if z ≥ 0, and 0 otherwise. Here, I is a function of the summed currents, including the neuromodulator-induced currents $I_x^{\prime }$s and ionotropic receptor-based currents (proportional to presynaptic neural firing rates). For simplicity, fast co-transmission of neurotransmitters was only considered in one modeling instance (co-release of 5-HT and glutamate via fast 5-HT₃ and ionotropic receptors) based on findings by Wang et al. (2019). From a computational modeling perspective, the DRN (Type I) 5-HT and Glu neuronal populations, which have rather similar activity profiles (Figure 1B, 3rd and 6th rows) could also be effectively grouped and considered as a single 5-HT neuronal population that “co-transmit” both 5-HT and Glu to DA neurons (Wang et al., 2019).

For each model circuit architecture that we investigated, we considered separately the inclusion of either Type I or II 5-HT neurons in the circuit (Cohen et al., 2015), and the possibility of excitatory and inhibitory projections from 5-HT to DRN Glu/GABA and DA neurons, and from DA to GABA neurons in the VTA. To allow tractability in the search for the many possible connectivity structures, the models’ connections were largely based on known experimental evidence. For example, the connections between VTA GABAergic and DRN Glu neurons were not considered as, to date, there is little experimental support. We also focused only on learned rewards (with reward-predicting cue followed by reward outcome) and unexpected punishment conditions, simulated using a combination of tonic and/or phasic afferent inputs. It should be noted that we did not consider other conditions and network learning effects (e.g., Hu, 2016) as the main aim was to demonstrate the plausibility of DRN-VTA circuit degeneracy and stability, constrained by specific observed neuronal signaling in reward and punishment tasks. In particular, we investigated various DRN-VTA circuit architectures with network activity profiles that closely resembled those shown in Figure 1B. Note also that all activity profiles in Figure 1B were based on experimental studies except VTA Glu neural activity, which we deduced to be similar to that of DA neural activity in the rewarding task (McDevitt et al., 2014) and assumed to be non-responsive in the punishment task. See below for further details regarding the modeling procedures, parameters, simulations, and analyses.

2.2. Input-output functions of neural population firing rates

The computational models developed were based on our previous mean-field, neural population-based modeling framework for neuromodulator circuits (Joshi et al., 2017), in which the averaged concentration releases of neuromodulators (e.g., 5-HT) were monotonic functions of the averaged firing rate of (e.g., 5-HT) neuronal populations via some neuromodulator induced slow currents. All five neural populations’ firing rates were described by threshold-linear functions [general form in Eq. (3)] (Jalewa et al., 2014; Joshi et al., 2017):

$$F_{5-HT}=g_{5-HT}{[I_{5-HT}-I_{0,5-HT}]}_{+}(4)$$

$$F_{DA}=g_{DA}{[I_{DA}-I_{0,DA}]}_{+}(5)$$

$$F_{Glu}=g_{Glu}{[I_{Glu}-I_{0,Glu}]}_{+}(6)$$

$$F_{GABA-DRN}=g_{GABA-DRN}{[I_{GABA-DRN}-I_{0,GABA-DRN}]}_{+}(7)$$

$$F_{GABA-VTA}=g_{GABA-VTA}{[I_{GABA-VTA}-I_{0,GABA-VTA}]}_{+}(8)$$

where [z]₊ = z if z ≥ 0, and 0 otherwise. The threshold input values for I_0,DA was –10 (a.u.) for DA neurons, and I_0,5–HT was 0.13 (a.u.) for 5-HT neurons, to allow spontaneous activities mimicking in vivo conditions. A 5-HT neurons had a threshold-linear function with a gain value g_5–HT of about 1.7 times higher than that for DA neurons, and so we set there for DA and 5-HT neurons to be 0.019 and 0.033 (Hz), respectively (e.g., Shepard and Bunney, 1991; Richards et al., 1997; Crawford et al., 2010; Wong-Lin et al., 2012; Challis et al., 2013). For simplicity, we assumed the same current-frequency or input-output function in either tonic or phasic activity mode (Jalewa et al., 2014; Joshi et al., 2017).

2.3. Afferent currents and connectivity

The afferent current, I, for a neural population consisted of summed contributions from external excitatory inputs I_ext including those induced by reward or aversive stimuli, and recurrent interactions with other neural populations (see below) e.g., I_5–HT,DA for effective DA-induced currents in 5-HT neurons. Additionally, for a neuromodulator population, auto receptor-induced current, I_auto, was included.

Due to limited experimental evidence, and to reduce the parameter search space, we did not consider the following connections from (i) DRN GABA to VTA DA neurons; (ii) DRN GABA to VTA GABA neurons; (iii) VTA GABA to DRN Glu neurons; (iv, v) DRN Glu to DRN GABA neurons, and vice versa; and (vi) VTA DA to DRN Glu neurons. Then the total (population-averaged) afferent input currents to DA and 5-HT neurons were, respectively, described by:

$$I_{DA}=-I_{auto,DA}\pm I_{DA,5-HT}+I_{DA,Glu}-I_{DA,GABA-VTA}+I_{DA,ext}(9)$$

$$\begin{array}{cc}{I}_{5-HT}=\hfill & -I_{auto,5-HT}+I_{5-HT,DA}+I_{5-HT,Glu}-I_{5-HT,GABA-DRN}\hfill \\ \hfill & -I_{5-HT,GABA-VTA}+I_{5-HT,ext}\hfill \end{array}(10)$$

where the first terms on the right-hand sides of Eqs. (6) and (7) were autoreceptor-induced self-inhibitory currents, and the second terms were the 5-HT-to-DA (labeled with subscript DA,5−HT) and DA-to-5-HT (with subscript 5−HT,DA) interactions, the third terms were excitatory interactions from DRN Glu neurons, the fourth/fifth terms were inhibitory interactions from local DRN/VTA GABAergic neurons, and the last terms were additional external constant biased inputs from the rest of the brain and the influence of behaviorally relevant stimuli (due to rewards or punishments; see below) 5-HT neurons have a possible additional negative interaction from VTA GABA neurons (second last term on the right-hand-side) (Li et al., 2019). Negative or positive signs in front of each term indicated whether the interaction was effectively inhibitory or excitatory. The ± sign indicated effective excitatory (+) or inhibitory (−) interactions which we investigated (Figures 2–5), given their mixed findings in the literature [see also Eqs. (11–13)]. This form of summed currents was consistent with some experimental evidence that showed different afferents modulating the tonic and phasic activation (e.g., Floresco et al., 2003; Tian et al., 2016).

FIGURE 2

Figure 2. A parsimonious, sparsely connected DRN-VTA circuit model. Gray: Brain region. Colored circle: Neuronal population. Legend: network’s afferent inputs. Model architecture implicitly encompasses either Type I or II 5-HT neurons with two different inputs for reward/punishment task (bright red arrows if Type I; blue arrows if Type II; black arrows denote common inputs for reward/punishment task for both Types). Circuit connections: triangular-end arrows (excitatory); circle-end arrows (inhibitory). Thicker arrows: Stronger connection weights. Constant long-term reward inputs simultaneously to 5-HT and DA neurons to alter baseline activities. Sustained activity for the expectation of reward outcome implemented with tonic input between cue and reward outcome. All other inputs are brief, at cue or reward/punishment outcome, producing phasic excitations/inhibitions. Note: Self-inhibitory (self-excitatory) connections within GABAergic (Glu) neurons, and auto-receptor inhibitions within 5-HT or DA neurons were implemented but not shown here (see section “2. Materials and methods” and Supplementary Figure 1). This is the most basic, sparsely connected model architecture considered.

FIGURE 3

Figure 3. DRN-VTA model replicates signaling patterns and suggests multiple parallel circuits. (A,B) Model with reward (black dashed lines) and punishment (orange bold lines) tasks with 5-HT neurons that are of Type I (A) or Type II (B). This will be used as a standard activity profile template to test for degeneracy. Time label from cue onset. Green (red) vertical dashed-dotted lines: Cue (outcome) onset time (as in Figure 1B). Top-to-bottom: VTA DA, DRN 5-HT, DRN GABAergic, DRN Glu, and VTA GABAergic neural populations. (C) Hypothesis for multiple different DRN-VTA circuits operating in parallel, which may consist of different clusters of neuronal sub-populations and different sets of afferent inputs, leading to different outputs. Vertical dots denote the potential of having more than two distinctive circuits.

FIGURE 4

Figure 4. Neural circuit model architectures with similar network activity profiles. The activity profiles were similar to that in Figures 2A, B, satisfying the inclusion criterion. Model architectures (A–K) denote architectures of decreasing connectivity, with Figure 2 as architecture (K). Model architecture (L) has an asterisk to denote that it was the only model with fast 5-HT to VTA DA connection, simulating fast 5-HT₃ or Glu receptor-mediated connection or their combination (co-transmission). Architectures (A,L) have additional inhibitory input to DRN GABA neurons in reward tasks. All labels, connections, and nomenclature have the same meaning as that in Figure 3, except that, for simplicity, the relative connection weights (thickness) are not illustrated, and the diamond-end arrows denote connections that are either excitatory or inhibitory, with both explored. Self-connectivity is not illustrated [see the general example in Supplementary Figure 1 for a detailed version of model architecture (A)]. Each architecture consists of several distinct model types (with a total of 84 types) with different 5-HT neuronal or excitatory/inhibitory connectivity types (see Supplementary Tables 2–4 and Supplementary Figure 5).

FIGURE 5

Figure 5. Negative real eigenvalues at steady states of degenerate models. (A) Complete set of the real part of the eigenvalues for model #1 (Supplementary Tables 2–4) with architecture “A” in Figure 4. Horizontal axis: Eigenvalues ranked from the largest to the smallest (magnitude wise). Blue (red): More negative eigenvalues with phasic (blue) than tonic (red) input. Asterisk: Maximal eigenvalue (largest magnitude) for each condition. (B) For each of the 84 models, only the real part of the eigenvalue with the largest magnitude is plotted under phasic (blue cross) and tonic (red circle) input conditions. Model architectures “A” to “L” refer to the different architectures as in Figure 4, in which each has its own distinctive model types (e.g., different 5-HT neuronal or excitatory/inhibitory connectivity types). Eigenvalues for all model types have negative real parts, indicating dynamically stable. For most models, the eigenvalues are generally more negative during phasic than tonic activities.

Similarly, the total (population-averaged) afferent current to the glutamatergic (Glu), and VTA and DRN GABAergic neurons can, respectively, be described by:

$$I_{Glu}=I_{self,Glu}\pm I_{Glu,5-HT}+I_{Glu,ext}(11)$$

$$\begin{array}{cc}{I}_{GABA-DRN}\hfill & =-I_{self,GABA-DRN}\pm I_{GABA-DRN,5-HT}+I_{GABA-DRN,DA}\hfill \\ \hfill & -I_{GABA-DRN,GABA-VTA}+I_{GABA-DRN,ext}\hfill \end{array}(12)$$

$$\begin{array}{cc}{I}_{GABA-VTA}\hfill & =-I_{self,GABA-VTA}+I_{GABA-VTA,5-HT}\hfill \\ \hfill & \pm I_{GABA-VTA,DA}+I_{GABA-VTA,ext}\hfill \end{array}(13)$$

where the subscript Glu denoted DRN Glu neural population, and subscripts GABA-VTA and GABA-DRN denoted GABAergic neural populations in the VTA and DRN, respectively. The subscript self denoted self-connection.

The averaged synaptic currents of non-5-HT/DA ionotropic glutamatergic and GABAergic neurons, namely, I_DA,Glu, I_5–HT,Glu, I_{DA,GABA–VTA}, I_{5–HT,GABA–VTA}, I_{5–HT,GABA–DRN}, I_self,Glu, I_{self,GABA–DRN}, I_{self,GABA–VTA}, and I_{GABA–DRN,GABA–VTA} were typically faster than currents induced by (metabotropic) 5-HT or DA currents. Thus, we assumed the former currents to reach quasi-steady states and described (Jalewa et al., 2014) and represented by I_e/i = ±J_e/iF_e/i, where the subscript e/i denoted excitatory/inhibitory synaptic current, J the connectivity coupling strength, F the presynaptic firing rate for neural population e/i, and the sign ± for excitatory or inhibitory currents. Furthermore, dimensionless coefficients or relative connectivity weights, W′s (with values 0), were later multiplied to the above neuromodulator-induced current terms [right-hand-side of terms in Eqs. (9–13); see Eqs. (20–24)]. Both the J′s and W′s were allowed to vary to fit the network activity profiles of Figure 1B within certain tolerance ranges (see below) while exploring different neural circuit architectures (Figure 4 and see Supplementary Table 2 for specific values). The self-connection weights J′s within the DRN Glu, DRN GABA, and VTA GABA neurons were set at 0.5, 0.5, and 10, respectively, for all network activity’s response profiles.

Autoreceptor-induced currents were known to trigger relatively slow G protein-coupled inwardly-rectifying potassium (GIRK) currents (Tuckwell and Penington, 2014). For 5-HT_1A auto-receptors, the inhibitory current I_5−HT,autowas described by Ritter et al. (2008), Tuckwell and Penington (2014), Joshi et al. (2017).

$$\begin{array}{cc}{\mathrm{\tau }}_{auto,5-HT}\frac{d{I}_{auto,5-HT}}{dt}=\hfill & -I_{auto,5-HT}\hfill \\ \hfill & +\frac{k_{auto,5-HT}}{1+e^{-g_{auto,5-HT}(\left[5-HT\right]-{[5-HT]}_0)}}\hfill \end{array}(14)$$

and similarly, for DA auto-receptor-induced inhibitory current I_auto,DA:

$${\mathrm{\tau }}_{auto,DA}\frac{d{I}_{auto,DA}}{dt}=-I_{auto,DA}+\frac{k_{auto,DA}}{1+e^{-g_{auto,DA}(\left[DA\right]-{[DA]}_0)}}(15)$$

where τ_auto,5−HT was set at 500 ms (Joshi et al., 2017) and τ_auto,DA at 150 ms (Benoit-Marand et al., 2001; Courtney et al., 2012; Cullen and Wong-Lin, 2015). The threshold values [5−HT]₀ and [DA]₀ were set at 0.1μM. These parameters can be varied to mimic the effects of auto-receptor antagonists/agonists (Joshi et al., 2017). The gains g_auto,5–HT and g_auto,DA were set at 10 μM^–1 each, and k_auto,5−HT=k_auto,DA=80 a.u. These values were selected to allow reasonable spontaneous neural firing activities and baseline neuromodulator concentration levels (see below), similar to those observed in experiments.

Similarly, we assumed the sigmoid-like influence of [5−HT] ([DA]) on DA (5-HT) neural firing activities between the DRN and VTA populations such that the induced current dynamics could be described by Wang and Wong-Lin (2013), Joshi et al. (2017):

$${\mathrm{\tau }}_{DA,5-HT}\frac{d{I}_{DA,5-HT}}{dt}=-I_{DA,5-HT}+\frac{k_{DA,5-HT}}{1+e^{-g_{DA}(\left[5-HT\right]-{[5-HT]}_1)}}(16)$$

$${\mathrm{\tau }}_{5-HT,DA}\frac{d{I}_{5-HT,DA}}{dt}=-I_{5-HT,DA}+\frac{k_{5-HT,DA}}{1+e^{-g_{5-HT}(\left[DA\right]-{[DA]}_1)}}(17)$$

with the time constants τ_5−HT,DA = 1s and τ_DA,5−HT = 1.2s (Haj-Dahmane, 2001; Aman et al., 2007). We set k_5−HT,DA = k_DA,5−HT = 0.03 a.u. and g_5−HT=g_DA = 20μM^–1, [5−HT]₁ = 0.3nM,[DA]₁=0.1 nM such that the neural firing activities and baseline neuromodulator concentration levels were at reasonable values (see below) and similar to those in experiments (e.g., Bunin et al., 1998; Hashemi et al., 2011). For simplicity, we assumed Eqs. (16) and (17) to be applied equally to all targeted neural populations, but with their currents multiplied by their appropriate weights w′s (see above).

2.4. Release-and-reuptake dynamics of neuromodulators

The release-and-reuptake dynamics of 5-HT followed the form of a Michaelis–Menten equation (Bunin et al., 1998; Hashemi et al., 2011; Joshi et al., 2011, 2017):

$$\frac{d[5-HT]}{dt}={[5-HT]}_p{R}_{5-HT}-\frac{V_{max,5-HT}[5-HT]}{K_{m,5-HT}+[5-HT]}(18)$$

where [5−HT]_p = 0.08 nM was defined as the release per firing frequency (Joshi et al., 2011, 2017; Flower and Wong-lin, 2014) [value selected to fit reasonable baseline activities (Hashemi et al., 2011); see below], and the Michaelis–Menten constants V_max,5−HT = 1.3μM/s (maximum uptake rate) and K_m,5−HT = 0.17μM (substrate concentration where uptake proceeds at half of the maximum rate) were adopted from voltammetry measurements (Hashemi et al., 2011).

Similarly, the release-and-reuptake dynamics for DA were described by:

$$\frac{d[DA]}{dt}={[DA]}_p{R}_{DA}-\frac{V_{max,DA}[DA]}{K_{m,DA}+[DA]}(19)$$

where V_max,DA = −0.004μM/s and K_m,DA = 0.15μM (May et al., 1988). We set [DA]_p = 0.1 nM to constrain the ratio [DA]_p/[5−HT]_p = 1.25 (May et al., 1988; Bunin et al., 1998; Hashemi et al., 2011). For simplicity, we assumed Eqs. (18) and (19) to be applied equally to all targeted neural populations.

2.5. Reward and punishment conditions with type I and type II 5-HT neurons

To limit the size of the parameter search space, we focused on only the classical, fully learned reward conditioning task, and unexpected punishment task. For each simulated trial or realization within a set of conditions (reward/punishment, excitatory/inhibitory connection), we set the cue onset time at 4.5 s to allow the network time to stabilize and reach a steady state. The within-trial protocol for the external input current, I_ext, was implemented as a function of time t as followed, depending on the simulated conditions. Note that, for simplicity, all external input currents were assumed to be excitatory, regardless of reward or punishment task, unless stated.

For reward tasks with Type-I 5-HT neurons (i) constant input currents I_5–HT,ext and I_DA,ext to 5-HT and VTA DA neurons, respectively, of amplitude 50 a.u. to simulate long-term reward; (ii) brief 0.2 s pulse I_Glu,extto DRN Glu neurons of amplitude 1,000 a.u. at 4.5 s; (iii) constant step input current I_{GABA–VTA,ext}to VTA GABAergic neurons with amplitude 200 a.u. from 4.5 to 5.7s to simulate sustained activity; and (iv) no external input I_{GABA–DRN,ext}to DRN GABAergic neurons. For punishment task with Type-I 5-HT neurons (i) no external input to VTA DA, 5-HT neurons, and DRN Glu neurons; (ii) brief 0.2 s pulse to VTA (I_app) and DRN GABAergic (I_app5i) neurons with amplitude 1,000 a.u. at 5.7 s.

For reward task with Type II 5-HT neurons (i) constant input current to 5-HT and VTA DA neurons with amplitude 50 a.u. to simulate long-term reward; (ii) additional step input current to 5-HT neurons with amplitude 100 a.u. from 4.5 to 5.7 s to simulate sustained activity; (iii) brief 0.2 s pulse to DRN Glu neurons of 1,000 a.u. at 4.5 s; and (iv) no external input to VTA and DRN GABAergic neurons. For punishment task with Type II 5-HT neurons (i) no external input to VTA DA and GABAergic neurons, and DRN Glu neurons; (ii) brief 0.2 s pulse to DRN GABAergic and 5-HT neurons with amplitude of 1,000 a.u. at 5.7 s.

In addition, for transient inputs, we have used multiplicative exponential factors exp (-t/τ) with τ of 50 ms to smooth out activity time courses (e.g., afferent synaptic filtering), but they do not affect the overall results. To simulate long-term, across-trial reward/punishment signaling, we assumed a higher constant excitatory input to both 5-HT and DA neurons in reward than punishment trials. When searching for the neural circuit architecture using either Type I or II 5-HT neurons, we limited ourselves as much as possible to the same internal DRN-VTA circuit structure. This also reduced the complexity of the parameter search space.

2.6. Baseline neural activities and inclusion criterion

We define the neural circuit activities under baseline conditions (right before cue onset) to follow that in Figures 3A, B. Namely, the baseline firing rates for 5-HT, DRN GABA, DRN Glu, DA, and VTA GABA neurons in the punishment task were 3.0, 21.5, 4.1, 4.8, and 13.5 Hz, and those in the reward task were 4.5, 19.4, 4.1, 4.8, and 16.3 Hz, respectively (compared with Supplementary Table 3). Baseline [5−HT] and [DA] levels were constrained to be at 10 and 1.5 nM, respectively. However, it is known that these activities can vary widely across subjects, species, and studies (Supplementary Table 3).

While searching for degenerate neural circuit architecture or investigating substantial changes due to simulated D₂ agonists, we had to define acceptable ranges of neural activities to evaluate whether the variant neural circuit still behaved similarly to that of the model in Figures 3A, B. Specifically, an inclusion criterion was set that evaluated the time-averaged percentage change in the neural population activities of any new model with respect to that of the template activity profiles (Figures 3A, B); the time-averaged percentage change in activities of the DA, 5-HT, DRN GABA, VTA GABA, and Glu neural populations for any qualified (degenerate) circuits had to be less than 10, 10, 16, 16, and 10% from those of the template activity profiles, respectively. These mean percentage changes were calculated within a 3 s time duration, from 1 s before cue onset to 1 s after outcome onset, encompassing both baseline and stimulus-evoked activities. Any neural circuit architecture which did not satisfy the inclusion criterion for any neural population activity was discarded when in search of degenerate neural circuits or was highlighted (embedded in non-black region in Figure 6) when investigating simulated D₂ agonist (see below). As this part of the study was to prove a concept, specific values of the inclusion criterion were theoretical and not critical–the general results, i.e., overall trends, remained even if these values were altered.

FIGURE 6

Figure 6. D₂ receptor agonists can distinguish subsets of DRN-VTA neural circuits. Simulated drug administered during punishment and reward tasks with efficacy factor X increments of 1, 10, 40, 70, and 100 times. Baseline condition: X = 1. Colors other than black in the pie chart denote that at least one neural population activity in specific model architecture(s) (labeled in yellow, and as in Figure 4) has deviated (increased or decreased) beyond the inclusion criterion. Specific colors denote changes in specific neural population activities (see Supplementary Tables 5–9 for details.).

2.7. Simulating the effects of D2 agonist

DA and 5-HT induced currents can lead to overall excitatory or inhibitory effects, depending on receptor subtype(s) and the targeted neurons. In particular, DA enhances VTA GABAergic neuronal activity via D₂ receptors and depolarizes the membrane of 5-HT neurons (e.g., Haj-Dahmane, 2001; Aman et al., 2007; Ludlow et al., 2009; Courtney et al., 2012; Ford, 2014). DA also regulates the activity of other DA neurons via D₂ auto-inhibitory receptors (Adell and Artigas, 2004). To study how D₂ receptor-mediated drugs can affect DRN-VTA architecture differently, we simulated different D₂ agonist dosage levels simultaneously by multiplying the connection weights of D₂ receptor-mediated currents (see above; Figure 4; orange connections in Supplementary Figure 1) by a factor “X” (Figure 6) of 10, 40, 70, and 100 times (default of X=1). Then for each dosage, we separately observed the deviation in activity profiles for each neural population with respect to the allowed range. Again, we considered the time-averaged percentage changes due to the drug to be substantial if the mean percentage change in at least one neural population activity had violated the inclusion criterion (see above, and Supplementary Tables 3–8).

2.8. Network stability analysis

The five neural population firing rates [Eqs. (4–8)], when combined with their associated afferent currents [Eqs. (9–13)] with explicit parameter values and relative connectivity weights, W′s and J′s (Supplementary Figure 1 and Supplementary Table 2), can be rewritten, respectively, as:

$$\begin{array}{cc}{F}_{5-HT}=\hfill & g_{5-HT}[-W_{auto,5-HT}{I}_{auto,5-HT}+W_{5d}{I}_{5-HT,DA}\hfill \\ \hfill & +J_{55e}{I}_{5-HT,Glu}-J_{55i}{I}_{5-HT,GABA-DRN}\hfill \\ \hfill & -J_{5ii}{I}_{5-HT,GABA-VTA}+99.87+I_{5-HT,ext}]{}_{+}\hfill \end{array}(20)$$

$$\begin{array}{cc}{F}_{DA}=\hfill & g_{DA}[-W_{auto,DA}{I}_{auto,DA}\pm W_{d5}{I}_{DA,5-HT}\hfill \\ \hfill & +J_{d5e}{I}_{DA,Glu}-J_{di}{I}_{DA,VTA-GABA}+210+I_{DA,ext}]{}_{+}\hfill \end{array}(21)$$

$$F_{Glu}=g_{Glu}{[J_{self,Glu}{F}_{Glu}\pm W_{5e5}{I}_{Glu,5-HT}+100+I_{Glu,ext}]}_{+}(22)$$

$$\begin{array}{cc}{F}_{GABA-DRN}=\hfill & g_{GABA-DRN}[-J_{self,GABA-DRN}{F}_{GABA-DRN}\pm \hfill \\ \hfill & W_{5i5}{I}_{GABA-DRN,5-HT}+W_{5id}{I}_{GABA-DRN,DA}\hfill \\ \hfill & -W_{5ii}{I}_{GABA-DRN,GABA-VTA}+450\hfill \\ \hfill & +I_{GABA-DRN,ext}]{}_{+}\hfill \end{array}(23)$$

$$\begin{array}{cc}{F}_{GABA-VTA}=\hfill & g_{GABA-VTA}[-J_{self,GABA-VTA}{F}_{GABA-VTA}\hfill \\ \hfill & +W_{i5}{I}_{GABA-VTA,5-HT}\pm W_{id}{I}_{GABA-VTA,DA}+\hfill \\ \hfill & 200+I_{GABA-VTA,ext}]{}_{+}\hfill \end{array}(24)$$

We also inserted the explicit parameter values to the dynamical equations [Eqs. (14–19)] to obtain:

$$500\frac{d{I}_{auto,5-HT}}{dt}=-I_{auto,5-HT}+\frac{80}{1+e^{-10(\left[5-HT\right]-0.1)}}(25)$$

$$150\frac{d{I}_{auto,DA}}{dt}=-I_{auto,DA}+\frac{80}{1+e^{-10(\left[DA\right]-0.1)}}(26)$$

$$1200\frac{d{I}_{DA,5-HT}}{dt}=-I_{DA,5-HT}+\frac{0.03}{1+e^{-20(\left[5-HT\right]-0.1)}}(27)$$

$$1000\frac{d{I}_{5-HT,DA}}{dt}=-I_{5-HT,DA}+\frac{0.03}{1+e^{-20(\left[DA\right]-0.3)}}(28)$$

$$\frac{d\left[5-HT\right]}{dt}=\frac{0.08F_{5-HT}}{1000}-\frac{0.0013\left[5-HT\right]}{0.17+\left[5-HT\right]}(29)$$

$$\frac{d\left[DA\right]}{dt}=\frac{0.1F_{DA}}{1000}-\frac{0.004[DA]}{0.15+[DA]}(30)$$

To check for network stability for each of the considered degenerate neural circuits, we first find each network’s steady state (or fixed point) by setting the rate of change for all the dynamical equations to zero, i.e., $\frac{d{I}_{auto,5-HT}}{dt}=\frac{d{I}_{auto,DA}}{dt}=\frac{d{I}_{DA,5-HT}}{dt}=\frac{d{I}_{5-HT,DA}}{dt}=\frac{d\left[5-HT\right]}{dt}=\frac{d[DA]}{dt}=0$, and then solving them algebraically. The solution of these equations will give the steady-state value (equilibrium point) of the system. Specifically, the currents (dynamical variables) from Eqs. (25–28) (e.g., $I_{auto,5-HT}=\frac{80}{1+e^{-10(\left[5-HT\right]-0.1)}}$) were substituted into Eqs. (20–24). See Supplementary material, for more detailed mathematical derivation.

Using only the linear parts of the above threshold-linear functions (which were validated post hoc), and after some algebraic manipulations, we obtained the following system of equations, in matrix form:

$$\begin{array}{c}\left(\begin{array}{c}\hfill \begin{array}{c}\hfill F_{Glu}\hfill \\ \hfill F_{GABA-VTA}\hfill \end{array}\hfill \\ \hfill F_{GABA-DRN}\hfill \end{array}\right)\hfill \\ ={\left(\begin{array}{ccc}1-g_{Glu}{J}_{self\text{-}Glu}\mathrm{\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0}\hfill & & \\ \mathrm{0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}1}+g_{GABA\text{-}VTA}\mathrm{\hspace{1em}\hspace{1em}0}\hfill & & \\ {\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}J}_{self,GABA\text{-}VTA}\mathit{\hspace{1em}}\hfill & & \\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{g}_{GABA\text{-}DRN}{J}_{5ii}\hspace{1em}\hspace{1em}1+g_{GABA\text{-}DRN}\hfill & & \\ {\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}J}_{self,GABA\text{-}DRN}\hfill & & \end{array}\right)}^{-1}\hfill \\ \left(\begin{array}{c}\frac{0.03W_{5e5}{g}_{Glu}}{1+e^{-20([5-HT]-0.1)}}+\frac{0}{1+e^{-16([DA]-0.1)}}\hfill \\ +g_{Glu}(100+I_{Glu,ext})\hfill \\ \frac{0.03W_{i5}{g}_{GABA-VTA}}{1+e^{-20([5-HT]-0.1)}}+\frac{0.03W_{id}{g}_{GABA\text{-}VTA}}{1+a{e}^{-20([DA]-0.1)}}\hfill \\ +g_{GABA\text{-}VTA}(200+I_{GABA\text{-}VTA})\hfill \\ \frac{0.03W_{5i5}{g}_{GABA\text{-}DRN}}{1+e^{-20([5-HT]-0.1)}}+\frac{0.03W_{5id}{g}_{GABA\text{-}DRN}}{1+a{e}^{-20([DA]-0.1)}}\hfill \\ +g_{GABA\text{-}DRN}(450+I_{GABA\text{-}DRN})\hfill \end{array}\right)\hfill \end{array}(31)$$

$$\begin{array}{c}\mathbf{\left(}\begin{array}{c}\hfill F_{5-HT}\hfill \\ \hfill F_{DA}\hfill \end{array}\mathbf{\right)}\hfill \\ =\mathbf{\left(}\begin{array}{c}{g}_{5-HT}\mathbf{(}\frac{-80}{1+e^{-10\mathbf{(}[5-HT]-0.1\mathbf{)}}}+\frac{0.03W_{5d}}{1+a{e}^{-20\mathbf{(}[DA]-0.1\mathbf{)}}}\hfill \\ \hfill +99.87+I_{5-HT,ext}\mathbf{)}\hfill \\ -g_{DA}\mathbf{(}\frac{0.03W_{d5}}{1+e^{-20\mathbf{(}[5-HT]-0.1\mathbf{)}}}+\frac{80}{1+e^{-10\mathbf{(}[DA]-0.1\mathbf{)}}}\hfill \\ \hfill -210-I_{DA,ext}\mathbf{)}\hfill \end{array}\mathbf{\right)}\hfill \\ +\mathbf{\left(}\begin{array}{cc}\hfill g_{5-HT}\hspace{1em}0\mathit{\hspace{1em}}\hfill & \\ \hfill \mathrm{\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}}{\mathit{g}}_{DA}\hfill & \end{array}\mathbf{\right)}\mathbf{\left(}\begin{array}{ccc}\hfill J_{55e}\hspace{1em}-J_{5i}\hspace{1em}-J_{55i}\hfill & & \\ \hfill J_{d5e}\hspace{1em}-J_{di}\mathrm{\hspace{1em}\hspace{1em}0}\hfill & & \end{array}\mathbf{\right)}\mathbf{\left(}\begin{array}{c}\hfill F_{Glu}\hfill \\ \hfill F_{GABA-VTA}\hfill \\ \hfill F_{GABA-DRN}\hfill \end{array}\mathbf{\right)}\hfill \end{array}(32)$$

But by setting Eqs. (29) and (30) to be zero, we obtained $F_{DA}=\frac{40\left[DA\right]}{0.15+\left[DA\right]}$ and $F_{5-HT}=\frac{16.25\left[5-HT\right]}{0.17+\left[5-HT\right]}$ at steady state. Substituting these two into Eqs. (31) and (32), we can solve for the values of [DA] and [5-HT] at steady state.

Next, to check whether the system is stable at a steady state, we compute the 5-by-5 Jacobian matrix M_Jacobian (Strogatz, 2018) for the five dynamical equations [Eqs. (25–30)], which can be computed using partial derivatives on the right-hand-side of these equations:

$$\begin{array}{c}{M}_{Jacobian}\hfill \\ =(\begin{array}{cccc}\hspace{1em}-1/500\hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}0\hfill & & & \\ \hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}-1/150\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}0\hfill & & & \\ \hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} 0\hspace{1em}-1/1200 0\hfill & & & \\ \hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em} 0\hspace{1em}-1/1000\hfill & & & \\ \frac{0.008}{1000}{g}_{5-HT}\hspace{1em} \hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em} 0\hspace{1em}\frac{0.08}{1000}{g}_{5-HT}\hfill & & & \\ \hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\frac{0.1}{1000}{g}_{DA}\hspace{1em}-\frac{0.1}{1000}{g}_{DA}\hspace{1em}0\hfill & & & \end{array}\hfill \\ \hfill \begin{array}{c}\hfill \frac{80}{500}\frac{10e^{-10([5-HT]-0.1)}}{{(1+e^{-10([5-HT]-0.1)})}^2}\hfill \\ \hfill 0\hfill \\ \hfill \frac{0.03}{1200}\frac{20e^{-20([5-HT]-0.1)}}{{(1+e^{-20([5-HT]-0.1)})}^2}\hfill \\ \hfill 0\hfill \\ \hfill \frac{-0.0013}{0.17+[5-HT]}+\frac{0.0013[5-HT]}{{(0.17+[5-HT])}^2}\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill \begin{array}{c}\hfill 0\hfill \\ \hfill \frac{80}{150}\frac{10e^{-10([DA]-0.1)}}{{(1+e^{-10([DA]-0.1)})}^2}\hfill \\ \hfill 0\hfill \\ \hfill \frac{0.03}{1000}\frac{20e^{-20([DA]-0.3)}}{{(1+e^{-20([DA]-0.3)})}^2}\hfill \\ \hfill 0\hfill \\ \hfill -\frac{0.004}{0.15+[DA]}+\frac{0.004[DA]}{{(0.15+[DA])}^2}\hfill \end{array})\hfill \end{array}(33)$$

The eigenvalues of this Jacobian matrix were computed for each steady state for each model under each simulated condition (e.g., reward/punishment task). If the real parts of all the eigenvalues of the Jacobian matrix were negative at a given steady state, then the model was considered to be dynamically stable at that steady state (Strogatz, 2018).

3. Results

3.1. A DRN-VTA model can reconcile many signaling patterns

We began with a parsimonious, sparsely connected DRN-VTA model architecture and adjusted the strength of the afferent inputs and the internal connection weights of the network (section “2. Materials and methods”) such that the network activity profiles attained (Figure 3) were readily comparable to the stereotypical profiles as illustrated in Figure 1B. We started the simulations using a minimal network configuration, the most sparsely connected DRN-VTA circuit model investigated in this work (Figure 2). All other subsequent architectures considered would henceforth be derived from this basic architecture. This minimal network architecture readily recapitulated many of the neuronal signaling changes in the DRN and VTA in separate experimental studies, both for reward (Figure 3, black dashed) and punishment (Figure 3, orange bold) tasks. The neural activity generated from the model (Figures 3A, B) will be used as a standard activity profile template to later test for degeneracy. Deviations from this template will be later used to distinguish different subgroups of the DRN-VTA degenerate networks (see below).

In particular, in the punishment task with Type I 5-HT neurons, brief excitatory inputs to GABAergic neurons in the DRN and VTA led to their punishment-based phasic activation (Figure 3A, 3rd, and 5th rows, orange) (Tan et al., 2012; Li et al., 2016). To replicate the activity profiles for a circuit with Type II 5-HT neurons, brief excitatory inputs to the DRN GABAergic and 5-HT neurons were implemented instead, leading to punishment-based phasic activation for both (Cohen and Uchida, 2012; Tan et al., 2012; Li et al., 2016; Figure 3B, 2nd and 3rd rows, orange). With both Type I and II 5-HT neurons, there was also a phasic inhibition of the VTA DA activity via VTA GABAergic neurons (Figure 3, 1st row, orange), in line with findings from previous studies (Schultz et al., 1997; Cohen and Uchida, 2012; Tan et al., 2012; Watabe-Uchida et al., 2017). It should be noted that with Type II 5-HT neurons, the phasic activation of VTA GABAergic neurons was not as potent due to the filtering by the slow excitatory connection from 5-HT to VTA GABA neurons. Alternatively, a phasic input might instead be sent to VTA GABAergic neurons, leading to more potent activation of the latter.

In the reward task, to replicate a sustained Type II 5-HT neuronal reward signaling between cue onset and reward outcome (Figure 3B, 2nd row, black dashed) (Liu et al., 2014; Cohen et al., 2015; Li et al., 2016), tonic excitatory input to DRN Type II 5-HT was implemented. The resulting sustained 5-HT activity led to gradual suppression of DRN GABAergic activity (Figure 3B, 3rd row, black dashed) (Liu et al., 2014; Li et al., 2016) but also a gradual rise in VTA GABAergic activity regardless of reward outcome (Figure 3B, last row, black dashed) (Cohen and Uchida, 2012; Eshel et al., 2015). The effects on these two GABAergic populations were, respectively, due to 5-HT’s inhibitory connection to DRN GABAergic neurons and excitatory connection to VTA GABAergic neurons. With regard to the latter, evidence of such excitatory influence via 5-HT_2C receptors was reported (Valencia-Torres et al., 2017).

Excitatory connection from DRN Glu to VTA DA neurons was particularly strong in the models. This, together with no (or weak) 5-HT-to-DA connection, was consistent with previous work (McDevitt et al., 2014). In the rewarding task, brief excitatory input to the DRN Glu neuronal population led to its reward-sensitive phasic activation (Figures 3A, B, 4th row, black dashed). Alternatively, we could have directly implemented phasic excitatory input to the VTA DA neurons. In any case, the phasic activation of DA activity was consistent with the reported DA neuronal response to a fully learned reward-predicting cue (Schultz et al., 1997; Cohen and Uchida, 2012; Watabe-Uchida et al., 2017).

Variations in some of the model parameters could still recapitulate these profiles, exhibiting model robustness (Supplementary Table 2). Moreover, only modifications to the afferent inputs to DRN 5-HT and VTA GABA neurons (section “2. Materials and methods”) were needed to replicate the signaling of Type I or II 5-HT neurons (Figures 3A, B, respectively), while maintaining the same internal connectivity structure. For example, a lack of sustained reward-based activity of Type I 5-HT activity (Figure 3B, 2nd row, black dashed) required additional external input to sustain VTA GABAergic neural activity (Figure 3A, bottom row, black dashed; Figure 2). It should be noted that this is assumed that the activity profiles for these non-5-HT neurons were qualitatively similar, regardless of the 5-HT neuronal types (Figure 2).

To understand across-trial reward vs. punishment effects in the model, we implemented a higher constant excitatory input into both DRN 5-HT and VTA DA neurons under reward compared to punishment conditions (Figure 2, long black arrows). This particular model required differential inputs to 5-HT and DA neurons (section “2. Materials and methods”) such that the overall tonic 5-HT neural activity was higher for reward than punishment trials, while DA neural activity remained unchanged (Figure 3, black dashed vs. orange bold lines in top two rows), again consistent with experimental observation (e.g., Cohen et al., 2015). In the model, although both 5-HT and DA neurons directly received constant across-trial reward-based excitatory inputs, the indirect inhibitory pathway from 5-HT neurons through VTA GABA neurons onto VTA DA neurons nullified the overall effects on DA neurons (Figures 2, 3). In other words, increased firing of 5-HT neurons could be activating VTA GABAergic neurons to a level sufficient to inhibit VTA DA neurons and thereby canceling out the net long-term reward signals (Figure 2). Note that this is just one out of several possible models; e.g., a more complex case (with more model parameters) which we are not considering could be that the VTA DA, VTA, GABA, and DRN 5-HT neurons individually receive very different inputs.

However, under these conditions, with both Type I and II 5-HT neurons, the baseline DRN and VTA GABAergic activities in the rewarding task were slightly different than those in the punishment task (Figure 3, 3rd and 5th rows, black dashed vs. orange bold), which have yet to be observed in experiments. The model’s inability to recapitulate this specific phenomenon might perhaps suggest that there are more complex features in the system, such as further division of neuronal subgroups. This also holds for other model architectures (see below, Figure 4). For example, it could be possible that a subpopulation of DRN GABAergic neurons is directly connected to 5-HT neurons (as in Figure 1B), but not another DRN GABAergic neuronal sub-population, such that across-trial reward signal inputs are distributed differently than that of Figure 1B. Moreover, high chemical and functional diversity among DRN 5-HT neurons are now well recognized (Okaty et al., 2019). Hence, it might be possible that there could exist multiple neural circuits with different circuit architectures operating in parallel, as illustrated in Figure 3C. Another possibility could be due to merely different inputs from e.g., different sources.

In summary, we have shown that, under reward and punishment conditions, many of the observed signaling patterns in different DRN and VTA neuronal types could readily be reconciled within a single sparsely connected DRN-VTA circuit model. However, not all the signaling patterns could be captured, suggesting that multiple different neuronal sub-populations and circuits may be operating in parallel within the DRN-VTA system. Next, we shall investigate whether various different DRN-VTA neural circuits could produce the same output, i.e., be degenerate.

3.2. Multiple degenerate DRN-VTA circuits

To search for degenerate DRN-VTA circuit models, we evaluated various combinations of connections within and between the DRN and VTA, and where necessary, adjusted any afferent inputs (section “2. Materials and methods”). We used the network activity profile template as illustrated in Figures 3A, B as output “targets” to check whether other different circuits could replicate similar activity profiles.

Given the variability of neuronal firing rates reported in the literature (Supplementary Table 3), to be considered a degenerate neural circuit, we set an inclusion criterion that allowed the time-averaged percentage changes of the neural population firing activities to be within certain acceptable ranges as compared to the activities in Figures 3A, B. Specifically, the time-averaged percentage change in activities of the DA, 5-HT, DRN GABA, VTA GABA, and Glu neural populations for any qualified (degenerate) circuits had to be less than 10, 10, 16, 16, and 10% from those of the template activity profiles, respectively. These mean percentage changes were calculated within a 3 s time duration, from 1 s before cue onset to 1 s after outcome onset, encompassing both baseline and stimulus-evoked activities. Any neural circuit architecture which did not satisfy the inclusion criterion for any neural population activity was discarded. Note that here, we aimed to prove a concept and the specific values of the inclusion criterion were theoretical and not critical—the general results, i.e., overall trends, remained even if the values were altered (not shown).

Various neural circuits were created by systematic addition and modification of connections of our parsimonious, sparsely connected DRN-VTA circuit model (now presented as model architecture “K” in Figure 4). For each model architecture, model parameters were searched, in both reward and punishment tasks using both Types I and II 5-HT neurons, until the activity outputs lie within the inclusion criterion. Both excitatory and inhibitory connections from DRN 5-HT to DRN Glu/GABA and VTA DA neurons were also explored. Once the inclusion criterion was met, the model parameter values were varied to check for robustness (Supplementary Table 2). This was repeated for different model architectures until we reached a high connectivity model structure (architecture “A” in Figure 4).

Based on this extensive search process, we obtained a total of 84 different neural circuit model architectures that recapitulated the activity profiles in Figures 3A, B. Their high-level model architectures were illustrated in Figure 4 (see also Figure 5 and Supplementary Tables 2, 3). For instance, model architecture “A” in Figure 4 (see Supplementary Figure 1, for detailed architecture) actually consisted of 32 distinctive models with different 5-HT neuron types and connectivity signs (Supplementary Table 4). Interestingly, this architecture’s model parameters remained the same with either excitatory or inhibitory connections (Figure 2, diamond connections; Supplementary Table 2).

We observed that all models with connection from 5-HT to DA neurons were relatively weak or not required, consistent with previous work (McDevitt et al., 2014) [but see below, (di Giovanni et al., 2008; Wang et al., 2019)]. Furthermore, a substantial number of other connections in the DRN-VTA model were identified to be redundant or relatively weak, at least in the context of the conditions investigated (Figure 4). In particular, the relatively weak VTA DA-to-GABA connection was consistent with studies that showed either a weak or non-existent direct effect of VTA DA on VTA GABA neurons (Morales and Margolis, 2017).

With these degenerate models, we could also investigate the effects of specific DRN-VTA connectivity in which there are mixed findings or a lack of evidence in the literature. For example, the influence of 5-HT on DRN GABA neurons could be excitatory or inhibitory (Hernández-Vázquez et al., 2019). In the case when this connection was inhibitory, and not excitatory, we found it to involve several degenerate circuits (model architectures “B” to “L”). In another example, with model architecture “B,” when the connection from DRN 5-HT to DRN Glu neurons was inhibitory, it led to more restricted values in its connectivity strength (0–0.1 a.u.) and hence less robust than when it was excitatory (0–1 a.u.) (Supplementary Table 2). It has been reported that there are mixed findings of 5-HT on VTA DA neurons (de Deurwaerdère and di Giovanni, 2017). In our simulations, we found that an inhibitory connection from 5-HT to VTA DA neurons could lead to several degenerate models (architectures “A”-“D,” “I,” and “L”). Note that in model “L” (labeled with an extra asterisk in Figure 4), we also successfully simulated its fast connectivity version, mimicking either 5-HT₃ receptor-mediated transmission or 5-HT-glutamate co-transmission (Wang et al., 2019).

In terms of compensatory network effects, we found that upon the removal of the 5-HT to DRN GABA connection (from architecture “B” to “C”), be it excitatory or inhibitory, the excitatory connection from DRN Glu to 5-HT neurons had to be weakened by ∼16% to satisfy the inclusion criterion. With additional removal of the connection from VTA DA to VTA GABA neurons, we found that the strength of the inhibitory connection from 5-HT to DRN GABA neurons had to be increased by ∼150% (architecture “D”). Interestingly, unlike other models, only the models with architectures “A” and “L” had to include an inhibitory connection from VTA GABA to DRN GABA neurons, in which such a connection was observed (Li et al., 2019).

3.3. Degenerate DRN-VTA circuit models are dynamically stable

After identifying the theoretical existence of degenerate models, we used dynamical systems theory to determine whether they were dynamically stable, i.e., whether (local) perturbation from their steady states would eventually cause a return to their initial steady states (see section “2. Materials and methods”). Specifically, the stability of each neural circuit could be determined by first finding the possible steady state(s) [i.e., fixed point(s)]. This was achieved by setting all the dynamical (differential) equations to zero and finding the algebraic solutions for the dynamical variables (section “2. Materials and methods”). Then the eigenvalues of the system’s Jacobian matrix at the steady states were computed (see section “2. Materials and methods” for the mathematical derivation of the steady states and the Jacobian matrix). For a neural circuit to be dynamically stable, the real part of all the eigenvalues associated with the steady state has to be negative. This was exactly what we found for all the identified degenerate neural circuits in (Figures 4, 5; section “2. Materials and methods”), with no imaginary part in the eigenvalues.

Figure 5A displayed the complete set of the real part of the eigenvalues for model #1 with architecture “A” in Figure 4. This model had inhibitory connections from VTA DA to VTA GABA neurons and from 5-HT to VTA DA/Glu neurons, using Type I 5-HT neurons and under punishment conditions (Supplementary Table 4). It was observed that the eigenvalues with phasic input (blue) were generally larger, magnitude-wise, than those with tonic input (red). This was more pronounced for the eigenvalues with the largest magnitude (maximal eigenvalues) (e.g., see Figure 5A). Importantly, all the eigenvalues were negative, indicating a dynamically stable network model even in the presence of additional phasic stimulus input.

We repeated the analysis for all 84 models, under both phasic and tonic input conditions. This analysis was presented in Figure 5B only for the maximal eigenvalues (red circles and blue crosses). The non-maximal eigenvalues had similar trends across the other models. In general, with phasic activities (blue crosses), the models were more stable than with tonic activities (red circles). However, during phasic activations, there were 18 models with rather small (close to zero) eigenvalues (magnitude wise), albeit still negative. Hence, they are easy to lose the stability. This was not observed for tonic activations, where the (most negative) eigenvalues were found to hover within a small range of values (−0.017 to −0.016), except for models with architecture “L” (∼−0.03). In fact, the latter models, which were the only ones with a fast 5-HT-to-DA connection (Figure 5B, models #77–84), were the most stable under both phasic and tonic conditions. Furthermore, there was no difference identified between the excitatory (models #77–80) and inhibitory (models #81–84) connections. Moreover, most of their eigenvalues under phasic conditions were substantially more negative than their tonic counterparts.

3.4. D2-mediated drugs can distinguish some degenerate DRN-VTA circuits

Given the large number of degenerate and stable DRN-VTA circuits identified, how could one distinguish among at least some of them? To address this, we investigated the neural circuit responses to simulated dopaminergic D₂ receptor agonists due to the extensive D₂ receptor-mediated connectivity within the degenerate DRN-VTA circuits (Figure 4 and Supplementary Figure 1). In particular, in the degenerate models, D₂ receptor-mediated connections involved those from VTA DA neurons to DRN 5-HT, DRN GABA, and VTA GABA neurons, and also the self-inhibitory connection (D₂ auto-receptor-mediated) of DA neurons (see section “2. Materials and methods” for references). To mimic the effects in the model of a D₂ receptor agonist, we gradually increased the strengths equally on the connections mediated by D₂ receptors (Supplementary Figure 1, connections emanating from DA neurons) by some factor (X) and observed how the network activity changed.

As we gradually increased the strengths of these specific sets of connections, subsets of the degenerate models gradually behaved differently from the activity profile template in (Figures 3A, B, 6A, model architectures embedded in non-black regions), i.e., not satisfying the inclusion criterion, thus, possibly allowing us to distinguish them. (see Supplementary Tables 5–9 for details). When at a low dose with a 10-fold (X=10) increase in D₂-mediated connections, all models showed substantial changes to their DA activities, beyond the inclusion criterion, regardless of reward/punishment condition and 5-HT neuronal type composition. Hence, they could not be distinguished at this dosage.

For moderately higher D₂ agonist doses (X = 40 and 70), models with architecture “A” can be distinguished from others, with substantial alteration to the DRN GABA and 5-HT activities (Figure 6A, light blue). Models with architectures “B” and “C” were also beginning to reveal substantial activity changes under reward conditions with Type II 5-HT neurons. Furthermore, with an increased factor of 70 under punishment conditions, a model with architectures “D,” “E,” “H,” and “L” could be distinguished from others with additional changes to their 5-HT activities. With a factor of 100, a maximum of five subsets of model architectures could be distinguished based on an activity enhancement of different combinations of neuronal types. An interesting observation we found was that for the same composition of 5-HT neuronal type (I or II), the subsets of distinguishable neural circuits were different between reward and punishment tasks.

Overall, our drug simulation predicted that a gradual increment of the level of D₂ receptor activation could lead to differential enhancements of firing rate activities that could be used to distinguish subsets of the degenerate DRN-VTA circuits.

4. Discussion

In this work, we have shown that neural circuits which are sources of neuromodulators can degenerate. In particular, the work focuses on computationally modeling and analyzing DRN-VTA circuits, as the DRN and VTA share structural and functional relationships among their constituent neuronal types. In particular, these circuits are involved in the regulation of several cognitive, emotional, and behavioral processes, and are implicated in many common and disabling neuropsychiatric conditions (Muller and Jacobs, 2009). Moreover, their neuronal signaling in reward and punishment tasks are well studied (Hu, 2016).

In this work, we developed biologically based mean-field computational models of the DRN-VTA circuit (Figure 2) with several neuronal types and tested them under classic conditions of (learned) reward and (unexpected) punishment. The modeling was partially constrained by known connectivity within and between the DRN and VTA regions, and their inputs from multiple other brain regions, including mixed combinations of inputs (Watabe-Uchida et al., 2012, 2017; Ogawa et al., 2014; Pollak Dorocic et al., 2014; Beier et al., 2015; Tian et al., 2016; Ogawa and Watabe-Uchida, 2018). Our work demonstrated that degenerate and stable DRN-VTA neural circuits are theoretically plausible.

We found that a parsimonious, sparsely connected version of the DRN-VTA model could reconcile many of the diverse phasic and tonic neural signaling events reported in the DRN and VTA in punishment and reward tasks observed across separate experimental studies (Figures 1, 3A, B). This model was evaluated using Type I and Type II 5-HT neurons in the DRN as defined electrophysiologically in a previous study (Cohen et al., 2015). In the case of Type II 5-HT neurons in reward condition, the model suggested that sustained 5-HT neuron activity between cue and reward outcome (Cohen et al., 2015) would lead to the gradual inhibition of DRN GABA neuron activity and enhancement of VTA GABA neuron activity (Cohen and Uchida, 2012; Li et al., 2016). The sparsely connected model could also reproduce experimental observations (Cohen et al., 2015) of an increase in baseline firing of Type I 5-HT neurons across several trials in the rewarding task, without similar effects on VTA DA neurons, or in the punishment task (Figure 3). Thus, this specific model suggested that slow, across-trial reward-based excitatory inputs could potentially be directly targeted to both DRN 5-HT and VTA DA neurons and that inhibitory 5-HT to GABA to DA connectivity could cancel out the effects of the direct input to DA neurons, rendering only long timescale changes to the baseline activity of 5-HT neurons but not DA neurons (Figure 2). It should, however, be noted that this is just one possible theoretical speculation offered out of many. For instance, another trivial possibility would be simply to have different separate inputs to DRN 5-HT and VTA DA and GABA neurons.

Despite the success of readily recapitulating the main observed phenomena, this relatively simple model was unable to capture some specific activity profiles with Type I 5-HT neurons, namely, the differential baseline activities of VTA GABA neurons and DRN GABA neurons between reward and punishment conditions. We later found that this was the case even for the more complex neural circuit models. Thus, perhaps additional neuronal populations and DRN-VTA circuits could be operating in parallel (Figure 3C) or it could just be these neurons receive different inputs from e.g., different sources. Such parallel circuits could be validated experimentally in the future, for example, using gene-targeting of specific DRN and VTA neuron subtypes and projections. Indeed, there is increasing evidence that DRN 5-HT neurons are more chemically diverse than previously expected, and that there is a high level of functional diversity in output pathways of the DRN and VTA (Watabe-Uchida et al., 2012, 2017; Wong-Lin et al., 2012; Ogawa et al., 2014; Pollak Dorocic et al., 2014; Weissbourd et al., 2014; Beier et al., 2015; Fernandez et al., 2016; Tian et al., 2016; Morales and Margolis, 2017; Zhou et al., 2017; Ogawa and Watabe-Uchida, 2018; Ren et al., 2018; Okaty et al., 2019). Further features will need to be incorporated including additional pathways and mechanisms as they are being uncovered, and this could include VTA glutamatergic neurons (McGovern et al., 2021). Clearly, the number of possible connections and models will be increased.

To demonstrate degeneracy in the DRN-VTA system, we showed that several variants of the DRN-VTA circuit model could readily recapitulate the same neuronal signaling profiles, with only occasional slight changes made to their afferent inputs (Figure 4). Our results also suggest experimental studies on the DRN and/or VTA system, e.g., using tracing methods, to search for more than one neural circuit configuration.

Consistent with the previous work (McDevitt et al., 2014), the degenerate models showed a relatively weaker direct connection from DRN 5-HT to VTA DA neurons as compared to the connection from DRN Glu to VTA DA neurons. However, previous studies had also demonstrated a direct influence of DRN 5-HT on VTA DA neuron activity (e.g., de Deurwaerdère and di Giovanni, 2017). More recent work has shown that DRN 5-HT terminals in the VTA co-release glutamate and 5-HT, eliciting fast excitation (via ionotropic receptors) onto VTA DA neurons and increased DA release in the nucleus accumbens to facilitate reward (Wang et al., 2019). Hence, we also incorporated a model with a fast 5-HT-to-DA connection (model architecture ‘L*’ in Figure 4) and found such circuits to be plausible in terms of capturing the stereotypical reward and punishment signaling.

Next, we applied dynamical systems theory and showed that all the found degenerate circuits were dynamically stable (Figure 5). Interestingly, we found that model architecture ‘L*’ with fast 5-HT-to-DA connection was dynamically more stable than all other architectures investigated (Figure 5B). Future computational modeling work could explore the effects of co-transmission of neurotransmitters on neural circuit degeneracy and functioning and using more biologically realistic spiking neuronal network models across multiple scales (Wong-Lin et al., 2012, 2017).

Finally, we simulated the gradual increase in dopaminergic D₂ receptor activation, due to extensive D₂-mediated connections in the investigated circuits. This was done by increasing the connection strengths emanating from the VTA DA neurons. This allowed us to distinguish subsets of the degenerate DRN-VTA circuits by identifying substantial and differential deviations in specific neural populations’ activities (Figure 6). Interestingly, we found that for the same composition of 5-HT neuronal type (I or II), the subsets of distinguishable DRN-VTA circuits were different between reward and punishment tasks (Figure 6A). Future experimental work could verify this model prediction.

As neuromodulators can selectively regulate degenerate neural circuits (Marder et al., 2014; Cropper et al., 2016), our theoretical work showed that even if the neuromodulator sources vary in physical and functional forms, the effects of neuromodulation on targeted brain areas may remain the same. For example, the co-modulatory effects of DA and 5-HT on prefrontal cortical rhythms (Wang and Wong-Lin, 2013) may be robust to specific interactions between the VTA DA and DRN 5-HT neurons. From a more general perspective, our computational modeling and analytical framework could be applied to the study of degeneracy and stability of neural circuits involving the interactions of other neuromodulators such as norepinephrine/noradrenaline e.g., Joshi et al. (2017). It should be also noted that in the development of the models, we resorted to a minimalist approach by focusing only on sufficiently simple neural circuit architectures that could replicate experimental observations. Future modeling work may investigate the relative relevance of these connections with respect to larger circuits involving cortical and subcortical brain regions across multiple scales, especially during adaptive learning e.g., Wong-Lin et al. (2017) and Zhou et al. (2018).

Taken together, our computational modeling and analytical work suggest the existence of degeneracy and stability in DRN-VTA circuits. Some of these degenerate circuits can be dynamically more stable than others, and subsets of the degenerate circuits can be distinguished through pharmacological means. Importantly, our work opens up a new avenue of investigation on the existence, robustness, and stability of neuromodulatory circuits, and has an important implication for the stable and robust maintenance of neuromodulatory functions.

Data availability statement

The original contributions presented in this study are publicly available. This data can be found here: GitHub, https://github.com/ckbehera/degeneracy.

Author contributions

CB and KW-L conceptualized and designed the study. CB, AJ, D-HW, and KW-L acquired the data and conducted the analyses. KW-L supervised the study. TS provided guidance on the study. CB, AJ, and KW-L wrote the first draft of the manuscript. All authors interpreted the data and revised the manuscript.

Funding

CB was supported by the Ulster University Research Challenge Fund. AJ, TS, and KW-L were supported by the BBSRC (BB/P003427/1). KW-L and D-HW were jointly supported by the Royal Society—NSFC International Exchanges. KW-L was further supported by COST Action Open Multiscale Systems Medicine (Open Multi Med) supported by COST (European Cooperation in Science and Technology), and Northern Ireland Functional Brain Mapping Facility (1303/101154803) funded by the Invest NI and the University of Ulster. D-HW received further support from NSFC under grants 32171094 and 31511130066.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2022.950489/full#supplementary-material

References

Adell, A., and Artigas, F. (2004). The somatodendritic release of dopamine in the ventral tegmental area and its regulation by afferent transmitter systems. Neurosci. Biobehav. Rev. 28, 415–431. doi: 10.1016/j.neubiorev.2004.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Aman, T. K., Shen, R. Y., and Haj-Dahmane, S. (2007). D2-like dopamine receptors depolarize dorsal raphe serotonin neurons through the activation of nonselective cationic conductance. J. Pharmacol. Exp. Ther. 320, 376–385. doi: 10.1124/jpet.106.111690

PubMed Abstract | CrossRef Full Text | Google Scholar

Beier, K. T., Steinberg, E. E., DeLoach, K. E., Xie, S., Miyamichi, K., Schwarz, L., et al. (2015). Circuit architecture of VTA dopamine neurons revealed by systematic input-output mapping. Cell 162, 622–634. doi: 10.1016/j.cell.2015.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Benoit-Marand, M., Borrelli, E., and Gonon, F. (2001). Inhibition of dopamine release via presynaptic D2 receptors: Time course and functional characteristics in vivo. J. Neurosci. 21, 9134–9141. doi: 10.1523/JNEUROSCI.21-23-09134.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Boureau, Y. L., and Dayan, P. (2010). Opponency revisited: Competition and cooperation between dopamine and serotonin. Neuropsychopharmacology 36, 74–97. doi: 10.1038/npp.2010.151

PubMed Abstract | CrossRef Full Text | Google Scholar

Bunin, M. A., Prioleau, C., Mailman, R. B., and Wightman, R. M. (1998). Release and uptake rates of 5-hydroxytryptamine in the dorsal raphe and substantia nigra reticulata of the rat brain. J. Neurochem. 70, 1077–1087. doi: 10.1046/j.1471-4159.1998.70031077.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Challis, C., Boulden, J., Veerakumar, A., Espallergues, J., Vassoler, F. M., Pierce, R. C., et al. (2013). Raphe GABAergic neurons mediate the acquisition of avoidance after social defeat. J. Neurosci. 33, 13978–13988. doi: 10.1523/JNEUROSCI.2383-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, J. Y., Amoroso, M. W., and Uchida, N. (2015). Serotonergic neurons signal reward and punishment on multiple timescales. Elife 4:e06346. doi: 10.7554/eLife.06346

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, J., and Uchida, N. (2012). Neuron-type specific signals for reward and punishment in the ventral tegmental area. Nature 482, 85–88. doi: 10.1038/nature10754

PubMed Abstract | CrossRef Full Text | Google Scholar

Courtney, N. A., Mamaligas, A. A., and Ford, C. P. (2012). Species differences in somatodendritic dopamine transmission determine D2-autoreceptor-mediated inhibition of ventral tegmental area neuron firing. J. Neurosci. 32, 13520–13528. doi: 10.1523/JNEUROSCI.2745-12.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Crawford, L. K., Craige, C. P., and Beck, S. G. (2010). Increased intrinsic excitability of lateral wing serotonin neurons of the dorsal raphe: A mechanism for selective activation in stress circuits. J. Neurophysiol. 103, 2652–2663. doi: 10.1152/jn.01132.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Cropper, E. C., Dacks, A. M., and Weiss, K. R. (2016). Consequences of degeneracy in network function. Curr. Opin. Neurobiol. 41, 62–67. doi: 10.1016/j.conb.2016.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Cullen, M., and Wong-Lin, K. (2015). Integrated dopaminergic neuronal model with reduced intracellular processes and inhibitory autoreceptors. IET Syst. Biol. 9, 245–258. doi: 10.1049/iet-syb.2015.0018

PubMed Abstract | CrossRef Full Text | Google Scholar

de Deurwaerdère, P., and di Giovanni, G. (2017). Serotonergic modulation of the activity of mesencephalic dopaminergic systems: Therapeutic implications. Progress Neurobiol. 151, 175–236. doi: 10.1016/j.pneurobio.2016.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

di Giovanni, G., di Matteo, V., Pierucci, M., and Esposito, E. (2008). Serotonin–dopamine interaction: Electrophysiological evidence. Progress Brain Res. 172, 45–71. doi: 10.1016/S0079-6123(08)00903-5

CrossRef Full Text | Google Scholar

Doya, K. (2002). Metalearning and neuromodulation. Neural Netw. 15, 495–506. doi: 10.1016/S0893-6080(02)00044-8

CrossRef Full Text | Google Scholar

Edelman, G. M., and Gally, J. A. (2001). Degeneracy and complexity in biological systems. Proc. Natl. Acad. Sci. U. S. A. 98, 13763–13768. doi: 10.1073/pnas.231499798

PubMed Abstract | CrossRef Full Text | Google Scholar

Eshel, N., Bukwich, M., Rao, V., Hemmelder, V., Tian, J., and Uchida, N. (2015). Arithmetic and local circuitry underlying dopamine prediction errors. Nature 525, 243–246. doi: 10.1038/nature14855

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernandez, S. P., Cauli, B., Cabezas, C., Muzerelle, A., Poncer, J. C., and Gaspar, P. (2016). Multiscale single-cell analysis reveals unique phenotypes of raphe 5-HT neurons projecting to the forebrain. Brain Struct. Funct. 221, 4007–4025. doi: 10.1007/s00429-015-1142-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Floresco, S. B., West, A. R., Ash, B., Moore, H., and Grace, A. A. (2003). Afferent modulation of dopamine neuron firing differentially regulates tonic and phasic dopamine transmission. Nat. Neurosci. 6:968. doi: 10.1038/nn1103

PubMed Abstract | CrossRef Full Text | Google Scholar

Flower, G., and Wong-lin, K. (2014). Reduced computational models of serotonin synthesis, release, and reuptake. IEEE Trans. Biomed. Eng. 61, 1054–1061. doi: 10.1109/TBME.2013.2293538

PubMed Abstract | CrossRef Full Text | Google Scholar

Ford, C. P. (2014). The role of D2-autoreceptors in regulating dopamine neuron activity and transmission. Neuroscience 282, 13–22. doi: 10.1016/j.neuroscience.2014.01.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Griffiths, D. J., and Schroeter, D. F. (2018). Introduction to quantum mechanics. Cambridge, MA: Cambridge University Press. doi: 10.1017/9781316995433

CrossRef Full Text | Google Scholar

Grossman, C. D., Bari, B. A., and Cohen, J. Y. (2022). Serotonin neurons modulate learning rate through uncertainty. Curr. Biol. 32, 586–599.e7. doi: 10.1016/j.cub.2021.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Haj-Dahmane, S. (2001). D2-like dopamine receptor activation excites rat dorsal raphe 5-HT neurons in vitro. Eur. J. Neurosci. 14, 125–134. doi: 10.1046/j.0953-816x.2001.01616.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hashemi, P., Dankoski, E. C., Wood, K. M., Ambrose, R. E., and Wightman, R. M. (2011). In vivo electrochemical evidence for simultaneous 5-HT and histamine release in the rat substantia nigra pars reticulata following medial forebrain bundle stimulation. J. Neurochem. 118, 749–759. doi: 10.1111/j.1471-4159.2011.07352.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayashi, K., Nakao, K., and Nakamura, K. (2015). Appetitive and aversive information coding in the primate dorsal raphe nucleus. J. Neurosci. 35, 6195–6208. doi: 10.1523/JNEUROSCI.2860-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernández-Vázquez, F., Garduño, J., and Hernández-López, S. (2019). GABAergic modulation of serotonergic neurons in the dorsal raphe nucleus. Rev. Neurosci. 30, 289–303. doi: 10.1515/revneuro-2018-0014

PubMed Abstract | CrossRef Full Text | Google Scholar

Hocking, D. R., Bradshaw, J. L., and Fielding, J. (eds) (2019). Degenerative disorders of the brain. London: Routledge. doi: 10.4324/9781351208918

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, H. (2016). Reward and Aversion. Annu. Rev. Neurosci. 39, 297–324. doi: 10.1146/annurev-neuro-070815-014106

PubMed Abstract | CrossRef Full Text | Google Scholar

Jalewa, J., Joshi, A., McGinnity, T. M., Prasad, G., Wong-Lin, K., and Hölscher, C. (2014). Neural circuit interactions between the dorsal raphe nucleus and the lateral hypothalamus: An experimental and computational study. PLoS One 9:e88003. doi: 10.1371/journal.pone.0088003

PubMed Abstract | CrossRef Full Text | Google Scholar

Joshi, A., Wong-Lin, K., McGinnity, T. M., and Prasad, G. (2011). “A mathematical model to explore the interdependence between the serotonin and orexin/hypocretin systems,” in 2011 annual international conference of the IEEE engineering in medicine and biology society, (Boston, MA: IEEE), 7270–7273. doi: 10.1109/IEMBS.2011.6091837

PubMed Abstract | CrossRef Full Text | Google Scholar

Joshi, A., Youssofzadeh, V., Vemana, V., McGinnity, T. M., Prasad, G., and Wong-Lin, K. (2017). An integrated modelling framework for neural circuits with multiple neuromodulators. J. R. Soc. Interface 14:20160902. doi: 10.1098/rsif.2016.0902

PubMed Abstract | CrossRef Full Text | Google Scholar

Levitan, I. B. (1987). Neuromodulation: The biochemical control of neuronal excitability. Oxford: Oxford University Press.

Google Scholar

Li, Y., Li, C. Y., Xi, W., Jin, S., Wu, Z. H., Jiang, P., et al. (2019). Rostral and caudal ventral tegmental area GABAergic inputs to different dorsal raphe neurons participate in opioid dependence. Neuron 101, 748–761.e5. doi: 10.1016/j.neuron.2018.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Zhong, W., Wang, D., Feng, Q., Liu, Z., Zhou, J., et al. (2016). Serotonin neurons in the dorsal raphe nucleus encode reward signals. Nat. Commun. 7, 1–15. doi: 10.1038/ncomms10503

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Zhou, J., Li, Y., Hu, F., Lu, Y., Ma, M., et al. (2014). Dorsal raphe neurons signal reward through 5-HT and glutamate. Neuron 81, 1360–1374. doi: 10.1016/j.neuron.2014.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ludlow, K. H., Bradley, K. D., Allison, D. W., Taylor, S. R., Yorgason, J. T., Hansen, D. M., et al. (2009). Acute and chronic ethanol modulate dopamine D2-subtype receptor responses in ventral tegmental area GABA neurons. Alcoholism 33, 804–811. doi: 10.1111/j.1530-0277.2009.00899.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Marder, E. (2012). Neuromodulation of neuronal circuits: Back to the future. Neuron 76, 1–11. doi: 10.1016/j.neuron.2012.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Marder, E., O’Leary, T., and Shruti, S. (2014). Neuromodulation of circuits with variable parameters: Single neurons and small circuits reveal principles of state-dependent and robust neuromodulation. Annu. Rev. Neurosci. 37, 329–346. doi: 10.1146/annurev-neuro-071013-013958

PubMed Abstract | CrossRef Full Text | Google Scholar

Matias, S., Lottem, E., Dugué, G. P., and Mainen, Z. F. (2017). Activity patterns of serotonin neurons underlying cognitive flexibility. Elife 6:e20552. doi: 10.7554/eLife.20552

PubMed Abstract | CrossRef Full Text | Google Scholar

May, L. J., Kuhr, W. G., and Wightman, R. M. (1988). Differentiation of dopamine overflow and uptake processes in the extracellular fluid of the rat caudate nucleus with fast-scan in vivo voltammetry. J. Neurochem. 51, 1060–1069. doi: 10.1111/j.1471-4159.1988.tb03069.x

PubMed Abstract | CrossRef Full Text | Google Scholar

McDevitt, R. A., Tiran-Cappello, A., Shen, H., Balderas, I., Britt, J. P., Marino, R. A. M., et al. (2014). Serotonergic versus non-serotonergic dorsal raphe projection neurons: Differential participation in reward circuitry. Cell Rep. 8, 1857–1869. doi: 10.1016/j.celrep.2014.08.037

PubMed Abstract | CrossRef Full Text | Google Scholar

McGovern, D. J., Polter, Abigail, M., and Root David, H. (2021). Neurochemical signaling of reward and aversion to ventral tegmental area glutamate neurons. J. Neurosci. 41, 5471–5486. doi: 10.1523/JNEUROSCI.1419-20.2021

PubMed Abstract | CrossRef Full Text | Google Scholar

Morales, M., and Margolis, E. B. (2017). Ventral tegmental area: Cellular heterogeneity, connectivity and behaviour. Nat. Rev. Neurosci. 18:73. doi: 10.1038/nrn.2016.165

PubMed Abstract | CrossRef Full Text | Google Scholar

Muller, C. P., and Cunningham, K. A. (2020). Handbook of the behavioral neurobiology of serotonin. Cambridge, MA: Academic Press.

Google Scholar

Muller, C. P., and Jacobs, B. (2009). Handbook of the behavioral neurobiology of serotonin, handbook of behavioral neuroscience. Amsterdam: Elsevier.

Google Scholar

Ogawa, S. K., and Watabe-Uchida, M. (2018). Organization of dopamine and serotonin system: Anatomical and functional mapping of monosynaptic inputs using rabies virus. Pharmacol. Biochem. Behav. 174, 9–22. doi: 10.1016/j.pbb.2017.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogawa, S. K., Cohen, J. Y., Hwang, D., Uchida, N., and Watabe-Uchida, M. (2014). Organization of monosynaptic inputs to the serotonin and dopamine neuromodulatory systems. Cell Rep. 8, 1105–1118. doi: 10.1016/j.celrep.2014.06.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Okaty, B. W., Commons, K. G., and Dymecki, S. M. (2019). Embracing diversity in the 5-HT neuronal system. Nat. Rev. Neurosci. 20, 397–424. doi: 10.1038/s41583-019-0151-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Pollak Dorocic, I., Fürth, D., Xuan, Y., Johansson, Y., Pozzi, L., Silberberg, G., et al. (2014). A whole-brain atlas of inputs to serotonergic neurons of the dorsal and median raphe nuclei. Neuron 83, 663–678. doi: 10.1016/j.neuron.2014.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, J., Friedmann, D., Xiong, J., Liu, C. D., Ferguson, B. R., Weerakkody, T., et al. (2018). Anatomically defined and functionally distinct dorsal raphe serotonin sub-systems. Cell 175, 472–487.e20. doi: 10.1016/j.cell.2018.07.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Richards, C. D., Shiroyama, T., and Kitai, S. T. (1997). Electrophysiological and immunocytochemical characterization of GABA and dopamine neurons in the substantia nigra of the rat. Neuroscience 80, 545–557. doi: 10.1016/S0306-4522(97)00093-6

CrossRef Full Text | Google Scholar

Ritter, J., Lewis, L., Mant, T., and Ferro, A. (2008). A textbook of clinical pharmacology and therapeutics. Boca Raton, FL: CRC Press. doi: 10.1201/b13234

CrossRef Full Text | Google Scholar

Schultz, W., Dayan, P., and Montague, P. R. (1997). A neural substrate of prediction and reward. Science 275, 1593–1599. doi: 10.1126/science.275.5306.1593

PubMed Abstract | CrossRef Full Text | Google Scholar

Shepard, P. D., and Bunney, B. S. (1991). Repetitive firing properties of putative dopamine-containing neurons in vitro: Regulation by an apamin-sensitive Ca2+-activated K+ conductance. Exp. Brain Res. 86, 141–150. doi: 10.1007/BF00231048

PubMed Abstract | CrossRef Full Text | Google Scholar

Strogatz, S. H. (2018). Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering. Boca Raton, FL: CRC press. doi: 10.1201/9780429399640

CrossRef Full Text | Google Scholar

Tan, K. R., Yvon, C., Turiault, M., Mirzabekov, J. J., Doehner, J., Labouèbe, G., et al. (2012). GABA neurons of the VTA drive conditioned place aversion. Neuron 73, 1173–1183. doi: 10.1016/j.neuron.2012.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, J., Huang, R., Cohen, J. Y., Osakada, F., Kobak, D., Machens, C. K., et al. (2016). Distributed and mixed information in monosynaptic inputs to dopamine neurons. Neuron 91, 1374–1389. doi: 10.1016/j.neuron.2016.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuckwell, H. C., and Penington, N. J. (2014). Computational modeling of spike generation in serotonergic neurons of the dorsal raphe nucleus. Progress Neurobiol. 118, 59–101. doi: 10.1016/j.pneurobio.2014.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Valencia-Torres, L., Olarte-Sanchez, C. M., Lyons, D. J., Georgescu, T., Greenwald-Yarnell, M., Myers, M. G. J., et al. (2017). Activation of ventral tegmental area 5-HT2C receptors reduces incentive motivation. Neuropsychopharmacology 42, 1511–1521. doi: 10.1038/npp.2016.264

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D. H., and Wong-Lin, K. (2013). Comodulation of dopamine and serotonin on prefrontal cortical rhythms: A theoretical study. Front. Integr. Neurosci. 7:54. doi: 10.3389/fnint.2013.00054

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H. L., Zhang, S., Qi, J., Wang, H., Cachope, R., Mejias-Aponte, C. A., et al. (2019). Dorsal raphe dual serotonin-glutamate neurons drive reward by establishing excitatory synapses on VTA mesoaccumbens dopamine neurons. Cell Rep. 26, 1128–1142.e7. doi: 10.1016/j.celrep.2019.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Watabe-Uchida, M., Eshel, N., and Uchida, N. (2017). Neural circuitry of reward prediction error. Annu. Rev. Neurosci. 40, 373–394. doi: 10.1146/annurev-neuro-072116-031109

PubMed Abstract | CrossRef Full Text | Google Scholar

Watabe-Uchida, M., Zhu, L., Ogawa, S. K., Vamanrao, A., and Uchida, N. (2012). Whole-brain mapping of direct inputs to midbrain dopamine neurons. Neuron 74, 858–873. doi: 10.1016/j.neuron.2012.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Weissbourd, B., Ren, J., DeLoach, K. E., Guenthner, C. J., Miyamichi, K., and Luo, L. (2014). Presynaptic partners of dorsal raphe serotonergic and GABAergic neurons. Neuron 83, 645–662. doi: 10.1016/j.neuron.2014.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitacre, J. M. (2010). Degeneracy: A link between evolvability, robustness and complexity in biological systems. Theor. Biol. Med. Model. 7, 1–17. doi: 10.1186/1742-4682-7-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong-Lin, K., Joshi, A., Prasad, G., and McGinnity, T. M. (2012). Network properties of a computational model of the dorsal raphe nucleus. Neural Netw. 32, 15–25. doi: 10.1016/j.neunet.2012.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong-Lin, K., Wang, D. H., Moustafa, A. A., Cohen, J. Y., and Nakamura, K. (2017). Toward a multiscale modeling framework for understanding serotonergic function. J. Psychopharmacol. 31, 1121–1136. doi: 10.1177/0269881117699612

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, P., He, Y., Cao, X., Valencia-Torres, L., Yan, X., Saito, K., et al. (2017). Activation of serotonin 2C receptors in dopamine neurons inhibits binge-like eating in mice. Biol. Psychiatry 81, 737–747. doi: 10.1016/j.biopsych.2016.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, W., Li, Y., Feng, Q., and Luo, M. (2017). Learning and stress shape the reward response patterns of serotonin neurons. J. Neurosci. 37, 8863–8875. doi: 10.1523/JNEUROSCI.1181-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, H., Wong-Lin, K., and Wang, D. H. (2018). Parallel excitatory and inhibitory neural circuit pathways underlie reward-based phasic neural responses. Complexity 2018:4356767. doi: 10.1155/2018/4356767

CrossRef Full Text | Google Scholar

Zhou, L., Liu, M. Z., Li, Q., Deng, J., Mu, D., and Sun, Y. G. (2017). Organization of functional long-range circuits controlling the activity of serotonergic neurons in the dorsal raphe nucleus. Cell Rep. 18, 3018–3032. doi: 10.1016/j.celrep.2017.02.077

PubMed Abstract | CrossRef Full Text | Google Scholar