Neural Circuit Interactions between the Dorsal Raphe Nucleus and the Lateral Hypothalamus: An Experimental and Computational Study

Orexinergic/hypocretinergic (Ox) neurotransmission plays an important role in regulating sleep, as well as in anxiety and depression, for which the serotonergic (5-HT) system is also involved in. However, little is known regarding the direct and indirect interactions between 5-HT in the dorsal raphe nucleus (DRN) and Ox neurons in the lateral hypothalamus (LHA). In this study, we report the additional presence of 5-HT1BR, 5-HT2AR, 5-HT2CR and fast ligand-gated 5-HT3AR subtypes on the Ox neurons of transgenic Ox-enhanced green fluorescent protein (Ox-EGFP) and wild type C57Bl/6 mice using single and double immunofluorescence (IF) staining, respectively, and quantify the colocalization for each 5-HT receptor subtype. We further reveal the presence of 5-HT3AR and 5-HT1AR on GABAergic neurons in LHA. We also identify NMDAR1, OX1R and OX2R on Ox neurons, but none on adjacent GABAergic neurons. This suggests a one-way relationship between LHA’s GABAergic and Ox neurons, wherein GABAergic neurons exerts an inhibitory effect on Ox neurons under partial DRN’s 5-HT control. We also show that Ox axonal projections receive glutamatergic (PSD-95 immunopositive) and GABAergic (Gephyrin immunopositive) inputs in the DRN. We consider these and other available findings into our computational model to explore possible effects of neural circuit connection types and timescales on the DRN-LHA system’s dynamics. We find that if the connections from 5-HT to LHA’s GABAergic neurons are weakly excitatory or inhibitory, the network exhibits slow oscillations; not observed when the connection is strongly excitatory. Furthermore, if Ox directly excites 5-HT neurons at a fast timescale, phasic Ox activation can lead to an increase in 5-HT activity; no significant effect with slower timescale. Overall, our experimental and computational approaches provide insights towards a more complete understanding of the complex relationship between 5-HT in the DRN and Ox in the LHA.


Introduction
Mood and neuropsychiatric disorders such as depression have a close relationship with sleep disturbances, and are instantiated by the overlap of emotional processing and the sleep-wake regulation neuronal circuitries [1]. The neuropeptide hormone orexin/ hypocretin (Ox) has been known to regulate sleep and its deregulation is related to narcolepsy, and novel drugs to facilitate sleep induction by activating Ox receptors are currently under development [2,3,4]. Recent studies also suggest a role for Ox in depression, emotional processing, reward seeking behaviour and in the regulation of endocrine functions [5,6,7,8,9,10]. Ox neurons comprising of neuropeptides Ox A and Ox B are found predominantly in the lateral hypothalamus (LHA) [11,12] and are known to function through OX 1 R and OX 2 R G-protein coupled receptors, respectively [13,14,15,16].
The neurotransmitter/neuromodulator 5-hydroxytryptamine (5-HT) released by serotonergic neurons, substantially located in the midbrain's dorsal raphe nucleus (DRN) is often associated with mood and emotional processing, and its dysfunction is related to mood and neuropsychiatric disorders [17,18,19]. In addition, perturbations of 5-HT have also been found to influence sleep [18,20,21,22]. Numerous drugs to treat depression are already on the market that target 5-HT neurotransmission [23,24].
It is important to understand how drugs that target Ox and/or 5-HT systems alter neuronal activity and signal transmission in order to understand the underlying mechanisms of antidepressive and sleep inducing effects. Therefore, we set out to map what subtypes of 5-HT receptors are expressed by Ox neurons, and how neuronal transmission and signal transduction in neuronal circuits may be controlled by these receptors.
Up till now, only the (inhibitory) 5-HT 1A R has so far been found in LHA's Ox neurons [25,26,27]. In addition to their inhibitory (5-HT 1A ) autoreceptors [28], 5-HT also excites the GABAergic inhibitory neurons within the DRN for self-regulation [29,30,31,32,33]. Within the LHA, in addition to their selfexcitatory Ox autoreceptors [4], Ox neurons can send direct longrange excitation to 5-HT neurons, and the GABAergic neurons in the DRN [31], mediated by both OX 2 R and OX 1 R [30,32]. Thus, Ox can have both direct and indirect influences on the DRN's 5-HT neurons. However, it remains unknown whether 5-HT can reciprocally indirectly influence LHA's Ox neurons by influencing the LHA's GABAergic neurons, and whether this connection is effectively excitatory or inhibitory. It is also not known whether Ox can innervate its local GABAergic neurons similar to how 5-HT neurons excite their local GABAergic neurons [34].
Knowledge of direct and indirect circuit connections is important to provide a more complete understanding of diversified neural circuit dynamics and regulations within the DRN-LHA system [35,36]. Furthermore, it is not known how the interplay between fast synaptic transmission and slow currents induced by 5-HT and Ox can affect the relay of information in these circuits. In this work, primarily through immunofluorescence (IF) staining, we attempt to map out a more complete neural circuit between the principal (5-HT and Ox) and non-principal (GABAergic and glutamatergic) neurons. Based on some of these findings, we present a computational model to investigate possible neural circuit dynamics between the DRN and LHA. As our goal in the modelling is to understand how the various neural circuit architecture and connectivity timescale affect the DRN-LHA activity, we keep the model as simple as possible. Hence we make use of population-based or ''mean-field'' firing-rate type approach, which compromises between previous biophysical models [37,38,39,40] and more abstract mathematical models [41,42,43].
Our experimental results reveal various other 5-HT receptor subtypes expressed in Ox and GABAergic neurons in the LHA, provide more evidence to support a unidirectional relationship between these LHA's neurons, and suggest that Ox can project to DRN's 5-HT neurons indirectly through local non-5-HT neurons. Our computational modelling results show that if 5-HT is weakly excitatory or inhibitory on LHA's GABAergic neurons, the network can exhibit slow oscillation. This is not observed if the connection is strongly excitatory. Furthermore, we show the importance of the timescale for the Ox-to-DRN connection during transient behaviour.

Animals
Twelve-week-old C57BL/6 male mice (n = 4) were used for each qualitative experiment and n = 6 mice for each quantification experiment described in this study. An orexin/enhanced green fluorescent protein (Ox-EGFP) breeder pair was a kind gift from Prof. Takeshi Sakurai (Kanazawa University, Japan). Brain sections from these mice show green fluorescence in the Ox neurons when excited at 488 nm wavelength. Breeding was set up in-house and the male pups were aged to ten-twelve weeks before the start of the experiments. Animals were maintained on a 12/ 12 h light/dark cycle (lights on at 8:00 A.M., off at 8:00 P.M.), in a temperature-controlled room (21.561uC). Animals received food and water ad libitum. All animal experiments were licenced by the UK Home Office in accordance with the Animals (Scientific Procedures) Act of 1986 and in agreement with UK and EU laws.

Perfusion, Fixation and Sectioning
Mice were anaesthetised with pentobarbitone (0.3 ml; Euthanal, Bayer AG, Leverkusen, Germany) and perfused transcardially with 0.1 M PBS (pH 7.4) buffer followed by ice-cold 4% paraformaldehyde in PBS. The brains were removed and fixed in 4% paraformaldehyde for at least 24 hr and cryoprotected in a 30% sucrose solution in PBS overnight at 4uC. Brains were then snap frozen using Envirofreez, and coronal sections of 45 mm thickness were cut using a Leica cryostat. According to the mouse brain atlas by Paxinos and Franklin (2004), LHA and DRN sections were cut at a depth of 20.34 mm to 22.80 mm bregma and 24.04 mm to 25.20 mm bregma respectively. Sections were chosen according to the stereological rules, with the first section taken at random and every sixth section afterward. In the case of LHA, by taking every 6 th section (in total 54 sections per LHA per brain half), at least 10-11 sections were taken per immunostaining experiment (n = 4). In the case of DRN, by taking every 3 rd section (in total 25 sections per DRN per full brain), at least 11-12 sections were taken per immunostaining experiment (n = 4). 10-11 sections from 4 mice brain halves (in case of LHA) and full brain (in case of DRN) were processed so, at least 40 sections were considered in total for every immunostaining experiment.
All primary antibodies used in this study were previously characterized and their specificity was verified according to the respective manufacturer. After blocking in 1% BSA and 5% donkey normal serum in TBS buffer (pH 8, Sigma Aldrich) to avoid nonspecific antibody binding, sections were incubated in the primary antibody overnight at 4uC. The following day, sections were incubated in the secondary antibody for an hour at room temperature and mounted using Vectashield mounting medium (Vector Laboratories) on the slides coated with 3-aminopropyl triethoxy silane (Sigma Aldrich). For double IF stainings, a simultaneous method was used where sections were incubated with two primary antibodies together for 48 hrs at 4uC. In the case of triple IF stainings, three primary antibodies were added together and sections were incubated for 72 hrs at 4uC. Labeled donkey IgG (H+L) anti-goat Alexa Fluor 488 (1:800 dilution, Cat. # A11055, Molecular Probes), anti-rabbit Alexa Fluor 546 (1:1000 dilution, Cat. # A11056, Molecular Probes), anti-chicken CF 594 (1:1000 dilution, Cat. # BTIU20167, Biotium) and antiguinea pig CF 633 (1:1000 dilution, Cat. # BTIU20171, Biotium) secondary antibodies were used in the study. Negative controls were performed for single IF staining by omitting the primary antibody and for double/triple IF staining by omitting the primary and secondary antibody.

Antigen Retrieval
Antigen retrieval was done while performing NMDAR1 IF staining. Sections were incubated in 10 mM sodium citrate (pH 6) at 80uC for 30 min, before blocking the sections.

Microscopy
Fluorescence microscopy. All single IF stainings in EGFP brain sections for 5-HT receptors on the Ox neurons in the LHA were visualized using Qimaging (Chromaphor). Microscopy was performed using an Olympus BX51 (Surveyor version 5.5.5.30, automated specimen scanning for the OASIS automation control system).
Confocal microscopy. Imaging was performed using a confocal microscope (Leica Microsystems; SP5 LAS IF Software). For quantification experiments, three sections of similar density of Ox neurons in the LHA were analyzed per brain (n = 6). 4-5 images were obtained from each section thus, 70-90 images were analyzed for each quantification experiment.

Co-localization Quantification
For quantification of the 5HT receptor subtype on the Ox neurons in LHA, images were acquired using 636 objective in a Leica SP5 confocal microscope. Once the conditions such as photomultiplier gain for each channel and pinhole settings were adjusted to minimize background noise and saturated pixels, parameters were kept constant for all acquisitions. Triple-stained images were obtained by sequential scanning for each channel to eliminate the cross talk of chromophores and to ensure reliable quantification of co-localization. Ambiguity and inconsistency are the two major issues affecting colocalization analysis. In the context of digital imaging, colocalization means the colours emitted by the fluorescent molecules occupy the same pixel in the image [53,54]. Therefore, we have used the JACoP (Just Another Colocalization Plug-in) tool of Image J for colocalization analysis. The degree of Ox neuron (Alexa488, green) signal colocalizing with 5HT receptor (Alexa546, red) signal was quantified on single-plane 8-bit color images using the JACoP plugin [55]. A simple way of measuring the dependency of pixels in dual-channel images is to plot the pixel grey values of two images against each other. The intensity of a given pixel in the green image is used as the x-coordinate of the scatter plot and the intensity of the corresponding pixel in the red image as the ycoordinate. Results in JACoP are displayed in a pixel distribution diagram called a scatter plot or fluorogram in addition to the calculated co-localization coefficients such as Pearson's and Overlap coefficient.
Pearson's correlation coefficient (Rr) is the most quantitative estimate of colocalization that depends on the amount of colocalized signals in both channels in a nonlinear manner and is a well-defined and commonly accepted means for describing the extent of overlap between image pairs [56]. It is used for describing the correlation of the intensity distributions between channels. It takes into consideration only similarity between shapes, while ignoring the intensities of signals. The values of Pearson's coefficient range from 21 to 1, with values from 0.5 to 1.0 indicating colocalization and 21.0 to 0.5 indicating nocolocalization. As Pearson's Correlation does some averaging of pixel information and can return negative values another method, the Overlap Coefficient, is simultaneously used to describe overlap. Manders' overlap coefficient (R) is based on the Pearson's correlation coefficient with average intensity values being taken out of the mathematical expression. This new coefficient varies from 0 to 1, with values from 0.6 to 1.0 indicating colocalization and 0 to 0.6 indicating no colocalization. Overlap coefficient according to Manders indicates an overlap of the signals and thus represents the true degree of Colocalization. Costes randomization (number of randomization rounds = 1000) was used to exclude any co-localization of pixels that might have occurred due to chance [57]. P-value for each image pair was 100.0% (calculated from the fitted data).

Statistics
Statistical analyses were performed using Prism 5 (GraphPad software Inc. USA) with the level of probability set at 95% and the results are expressed as means6SEM. Data for 5-HT receptor quantification was analysed by two-tailed unpaired t-test.

Computational Model
To investigate the consequences of the DRN-LHA circuit architecture on systems dynamics, we implement neural network model that is an extension and modification of our previous model [58]. The aim of the model is to understand how the circuit connectivity and timescale affect the DRN-LHA activity.
Neural units. Our neural network model consists of 4 populations of neurons, namely, Ox neurons in the LHA, local LHA inhibitory GABAergic neurons, 5-HT neurons in the DRN, and local DRN inhibitory GABAergic neurons. Glutamatergic neurons will be ignored in this work primarily due to the evidence showing that glutamatergic effects in DRN are locally weaker when compared to local GABAergic influence [59], and that we can implicitly encompass the effects of the LHA's glutamatergic neurons on Ox neurons [4] with self-excitatory Ox connections. Furthermore, incorporating two additional (glutamatergic) neural populations can lead to more free parameters in the model.
The chosen model is of the population-averaged or ''meanfield'' firing rate type model [60,61,62,63]. This simplifies a population of neurons into its representative unit. The 4 neural populations to be considered are the LHA's Ox and GABAergic neural populations, and the DRN's 5-HT and GABAergic neural populations. With support from electrophysiological data, the input-output or current-frequency relationship (f-I curve) for each neural population can be described by threshold-linear functions [64,65]: where f j and (I Local,j -I Background,j ) denote the population-averaged firing rate activity and total afferent input of the j th neural population, respectively. I Background,j is the background current, consisting of inputs from the rest of the brain areas. g j and I 0,j determine the input-output slope and the current threshold, respectively. The threshold-linear function [z] + = z if z.0, and 0 otherwise. Based on their known neuronal electrophysiological properties [4,27,66]), we determine and constrain the values of the g j 's (i.e. g 5-HT and g GABA for DRN; g GABA and g Ox for LHA) and the I 0,j 's, while consider I Background,j as a free parameter. Note that we have assumed the dynamics of the neural population (neuronal membrane time constant , 10 ms) to be relatively instantaneous and slaved to the much slower timescale of the connections (, s) [63]. This simplification reduces 4 model parameters (neuronal membrane time constants) and 4 dynamical (differential) equations for the 4 neuronal types [60].
2.8.2. Inputs and connections: Since each neural unit receives inputs from self-feedback and from 2 other neural units (e.g. 5-HT neuronal population receives both afferent inputs from DRN's GABAergic neurons and longer range connection from Ox neurons), the local afferent input (minus the background input) can be described by: where the subscript ''self'' denotes a self-feedback connection (e.g. due to autoreceptors in Ox and 5-HT neural populations, or GABAergic synapses in the two inhibitory neural populations), and the subscript ''10 or ''20 denotes the afferent inputs due to the other 2 neural populations. For local GABAergic neurons, they receive self-inhibition, and projections from their local principal neural population (5-HT or Ox if from DRN or LHA, respectively), and long-range projection from the other brain region (LHA or DRN). If the net effect of the i th neural population on the j th neural population is inhibitory or excitatory, the coefficient in front of I j,i will be 21 or +1, respectively. For example, suppose the 5-HT neural population receives inhibitory autoreceptor influence, direct projection from Ox neurons, and local GABAergic influence. Then, To simulate the phasic response of the circuit, an extra term I stim is added to the 5-HT or/and Ox inputs The dynamics of each input are filtered by its (synaptic/ effective) time constant (t syn, j, i ) as follows: In principle, there are 10 such similar dynamical equations describing the currents caused by ionotropic (4 equations) and metabotropic (6 equations) receptors. The synaptic time constants associated with the ionotropic receptors are obtained from electrophysiological data while the effective metabotropic time constants are deduced from the associated G protein-coupled inwardly-rectifying potassium (GIRK) current, or if unavailable, the temporal change in firing rate due to the injection of 5-HT or Ox [67]. For simplicity, a simple linear relationship is assumed between the input current I j,i and activity level f i . The other important parameter, J syn,j,i is the connection strength within or between the neuronal populations. Under steady state condition (dI j,i /dt = 0), J syn,j,i is defined as the ratio of the current I j,i and the associated (presynaptic) activity f i . These currents are obtained from various experiments [4,26,27,31,68,69]. The in vivo baseline firing rates (f i ) for the Ox and GABAergic populations in LHA are ,5 Hz (3-8 Hz in experiments) [70,71,72]. In DRN, the baseline neuronal firing rate of 5-HT and GABAergic neuronal population is ,5 and ,15 Hz, respectively [73]. The relationship among these baseline activities will be used to constrain our model parameters, namely the J i , g i and I j,0 ( Table 1). Consistent with our assumption on ignoring the relatively much faster neuronal membrane dynamics (,10 s ms), we shall also ignore the dynamics for the relatively fast GABAergic synapses (,4 ms), assuming they attain instantaneous steady states. This further reduces 4 dynamical equations in describing the associated currents.
Thus Further details of the model parameter values, justifications, and their related references are summarized in Table 1.  Simulations and analysis. We use XPPAUT [74] for neural circuit dynamics simulations and stability analysis. The Runge-Kutta 2 numerical integration algorithm with a time step of 10 ms is used. Smaller time steps do not affect our results.     (Figure 2). This suggests that the DRN-to-LHA connection may transmit signals fast. Figure 3 shows representative z-stack max (30 mm thickness) confocal images showing colocalization of Orexin A with 5HT 1A R ( Figure 3A) and 5HT 3A R ( Figure 3B), respectively. As we demonstrated the additional 5HT receptor subtypes on the Ox neurons in LHA for the first time, the next obvious step was to quantify the double-labeled Ox neurons for each 5HT receptor subtype. Using the JACoP plug-in tool, we quantified the degree of co-localization by calculating Pearson's correlation coefficient ( Figure 4B) and overlap coefficient by Manders ( Figure 4C). It is important to note, that the two coefficients, namely Pearson's Correlation Coefficient and Overlap Coefficient, showed similar pattern of changes while revealing the different aspects of the colocalization process, proving the applicability of the calculations to investigate the degree of colocalization of serotonin receptor subtypes and orexin A (Ox neurons).

5-HT 3A R and 5-HT 1A R Receptors on Ox as well as GABAergic Interneurons in the LHA
To examine the serotonergic long-range connections from DRN in the midbrain to the Ox and GABAergic neurons in the LHA, triple-label IF staining was performed. We searched for 5-HT 3A R

OX 1 R, OX 2 R and NMDAR1 Receptors on Ox Neurons but not on GABAergic Neurons in the LHA
To identify the role of Ox neurons in regulating the activity of local GABAergic neurons, we carried out a triple IF labeling in the LHA. To study this inter-relationship, we searched for OX 1 R and OX 2 R receptors on the Ox and GABAergic neurons. We identified OX 1 R ( Figure 6A) and OX 2 R ( Figure 6B) on the Ox neurons but found these to be absent on the GABAergic interneurons. This sheds light on the self-regulatory mechanism of Ox neurons via OX 1 R and OX 2 R auto-receptors. Previous studies have demonstrated the presence of inhibitory GABA B receptors on the Ox neurons [33,75,76,77] and there is evidence that GABA A receptors are also present on Ox neurons (Backberg et al., 2004, Kokare et al., 2006. This could suggest a one-way  Ox neurons form a distinct group of a hypothalamic neuronal population that project to multiple brain regions and coordinate many physiological functions [9,78]. Also, there is a convergence of signals that regulate the activity of Ox neurons by neurotransmitters, hormones, etc. Glutamate is an important neurotransmitter for Ox neurons and excitatory AMPA receptors have been shown to mediate the miniature EPSC in Ox neurons [79]. Also, NMDA receptors activate Ox neurons in the perifornical region of LHA [80]. To verify whether NMDA receptors are present on the overall Ox neurons and not confined solely to the perifornical area, we did triple-label IF staining. Here, we show the presence of N-methyl D-aspartate receptor 1 (NMDAR1) on the Ox neurons ( Figure 7) and their absence on the GABAergic interneurons ( Figure 7A, overlay image).

Ox Axonal Projections Receive Glutamatergic and GABAergic Post-synaptic Inputs in the DRN
To study the long-range connections of Ox neurons in the DRN, we performed IF labeling in the DRN region for PSD-95 (a marker for glutamatergic synapses) and gephyrin (a marker for GABAergic synapses) post-synaptic proteins [81]. Firstly, we carried out double-label IF staining for PSD-95 and Ox in the DRN and observed PSD-95 immunopositive Ox axonal terminals ( Figure 8) indicating glutamatergic input to the Ox axonal projections. In another set of experiments, we examined GABAergic post-synaptic input to the Ox terminals as it has been shown that Ox fibres project to both 5-HT and GABAergic neurons in the DRN [31]. A double-label IF analysis for gephyrin and Ox in the DRN (Figure 9) showed gephyrin immunopositive

Map of the DRN-LHA Circuit
Based on the above results and in other previous studies (see below), a tentative map of the DRN and LHA is illustrated in Figure 10A. Connection (i) in the figure involves 5-HT 1A R [26], 5-HT 1B R, 5-HT 2A R, 5-HT 2C R, 5-HT 3A R (found in our current study). Although connection (ii), found in this study, consists of 5-HT 1A R and 5-HT 3A R, the specific types of connection (e.g. effectively excitatory or inhibitory) and their strengths are not known. Connection (iii) and its receptor types are not known yet. Similarly, the receptor types for connection (iv) (Ox 1 R and Ox 2 R found in this study are consistent with previous findings [30,32]), and connections (v) and (vi) (found in this study) are unknown. Connection (vii) receptors are not known and specifically Ox 1 R and Ox 2 R are not found in our study.
Other connections in the circuit based on other previous work: (viii) excitatory connection and receptor types are not known [82]; (ix) GABA A/B [33,75,76,77]; (x) AMPAR [79] and NMDAR1 (found in the current study); (xi) GABA B [77]; (xiii) Ox 2 R [4], Ox 1 R and Ox 2 R (found in the current studies); (xv) 5-HT 1B/1D [34]; 5-HT 1A/2A/2C [69]; (xvi) 5-HT 7 [34]; (xvii) GABA A/B [34]; (xviii) AMPAR and NMDAR [34]; (xix) AMPAR and NMDAR [34]; (xxi) 5-HT 1A R [83]. For the connections (xii), (xiv), (xx) and (xxii) (shown by black dotted circle/arrow in Figure 10A), we implicitly assume that non-principal neurons (GABAergic and glutamatergic neurons in LHA and DRN) are self-coupled.  The experimental findings in this study and in other previous work now provide sufficient information to build a computational model to investigate how the DRN-LHA neural architecture can influence the systems dynamics. We purposely simplified the implementation of the neural population units and connections in order to better illustrate the effects of network topology on dynamics (see Materials and Methods section). We also omit modelling the glutamatergic neurons in the DRN and LHA because there is evidence that glutamate has a weaker effect than GABA in the DRN [59]; and that indirect excitation from LHA's glutamatergic neurons on Ox neurons may be represented implicitly by the Ox autoreceptors. This significantly reduces the number of model parameters and dynamical equations involved. In our model, we find that the various receptor subtypes will not affect the steady-state activities of the system. We shall henceforth not elaborate on the specific receptor subtype till we investigate the transient activation part of the results (section 3.6.3) when we study the influence of various timescales induced by the various receptors.
Next we investigate how the unknown model parameters and different timescales affect the system behaviour. Specifically, the focus is to understand the effects of: (i) 5-HT on LHA's GABAergic neurons; (ii) local Ox and GABAergic interactions; and (iii) how connection timescales affect phasic 5-HT or Ox activations.

Oscillatory DRN-LHA Behaviour if 5-HT Weakly Excites or Inhibits LHA's GABAergic Neurons
It is not yet known whether the connection from 5-HT neurons to LHA's GABAergic neurons (J 5-HT-to-GABA(LHA) ) is effectively excitatory or inhibitory. Here, we shall explore these possibilities. If this connection is excitatory, we find that as its connection strength J 5-HT-to-GABA(LHA) is increased, the steady-state firing rate activities for most of the neural populations decrease (Figure 11). Note that LHA's GABAergic neural population barely increases. This slight change is due to the self-inhibition within these GABAergic neurons. In contrast, the relatively larger changes for the other neural populations (especially when the connection strength is weaker) are due to the strong inhibitory projection from the LHA's GABAergic neuron onto the Ox neurons (see Table 1), which subsequently affect the neurons in the DRN.
When the connection from 5-HT to LHA's GABAergic neurons becomes very weak (J 5-HT-to-GABA (LHA) close to 0), the DRN-LHA circuit will begin to exhibit slow oscillatory behaviour. Figure 12 shows such transition towards oscillatory behaviour for one sample neural population (the rest of the neural populations look similar). The oscillatory period can be as slow as a few minutes (Figure 12, inset). It can also be observed that as the connection strength decreases, the amplitude of oscillation (bounded by the top and bottom lines in grey region) increases. When this particular connection becomes inhibitory (J 5-HT-to-GABA(LHA) ,0), the oscillation amplitude can become very large (not shown). It is well-known that excitatory-inhibitory network can easily create oscillation [60,61]. Similarly, the observed oscillation phenomenon can be explained by the interplay between strong inhibition and selfexcitatory auto-regulation in the Ox neural population; the oscillation disappears when Ox autoreceptors are removed from the model e.g. when J Ox-auto = 0 (not shown). As far as we know, there has yet to be any observable slow oscillation of Ox in the timescale of minutes. Hence, our model suggests that strong excitatory connection from 5-HT to LHA's GABAergic neurons is more plausible, and we shall proceed with this assumption from here on.

DRN-LHA System is Resilient to Changes in the Ox-to-GABAergic Neurons in the LHA
One of our experimental findings in this study shows that the LHA's GABAergic neurons do not have Ox receptors, consistent with indirect results from previous works [4]. Using our model, we check for the significance of such Ox receptors' absence. In our model, we gradually increase the connection strength of Ox to LHA GABAergic neurons (J Ox-to-GABA(LHA) ). This only marginally decreases the DRN-LHA steady-state activities ( Figure 13). A similar explanation as that for the results in Figure 11 can be used to account for this phenomenon.

Transient Ox and 5-HT Activities can be Affected by Different Receptor Timescales
Having studied the tonic activities (stable steady states) of the LHA-DRN system, we shall now proceed to investigate how transient or phasic activations of 5-HT or Ox can affect the system. The motivation for this is that it has been known that 5-HT and Ox can phasically activate in the presence of behaviourally relevant stimuli [84,85,86]. Furthermore, our current experimental finding suggests that 5-HT can influence LHA neurons over multiple timescales, through both the slow G-protein-coupled (non-5-HT 3A ) and fast ligand-gated (5-HT 3A ) receptors. To simulate this, the (more plausible) excitatory connection J 5-HT-to-GABA(LHA) is considered, and a pulse stimulus current of the duration of 0.5 s with an amplitude of 150 pA is applied to either the 5-HT or Ox neural populations. To understand the individual role of the slow and fast timescales, and minimize any confounding effect, we simulate the two timescales separately. Figure 14A (left) shows that a 0.5 s stimulation of 5-HT neurons rapidly increase its activity, followed by a slower decay back towards baseline. The 5-HT activity decay is due to feedback inhibition from its autoreceptors and the local GABAergic neurons (hence the undershoot below baseline). The Ox activity responses are generally suppressed by the strong 5-HT-mediated inhibition ( Figure 14A, right). When the 5-HT-to-Ox and 5-HT-to-GABA (LHA) connections are fast (mimicking 5-HT 3A R), the rebound upon removal of stimulus (, after the 2 s mark in Figure 14A, right) is higher, suggesting the faster disinhibition than that for the slower connections. Varying these timescales do not affect the transient 5-HT activity in DRN (overlapping curves, Figure 14A, left).
When a similar stimulus is applied to the Ox neurons, the activity of the Ox neurons increases considerably throughout the stimulus duration ( Figure 14B, right), which is due to the selfamplifying effect of Ox autoreceptors dominating over the local GABAergic inhibition (feedback inhibition is not strong due to the weak excitatory connection from Ox neurons to their local GABAergic neurons). These Ox neurons, affect 5-HT neurons either directly or indirectly (through the DRN's GABAergic neurons). From our simulations, we find that when the Ox-to-5-HT connection acts on a fast timescale (5 s), they can transiently excite the 5-HT neurons ( Figure 14B, left). If the direct timescale is  much slower (e.g. 60 s), then 5-HT neurons can hardly be activated. Varying Ox-to-GABA(DRN) connection does not affect the system significantly. Note the post-stimulus undershoot due to recurrent inhibition, suppressing the Ox and 5-HT transiently below baseline. Taken together, when Ox-to-5-HT acts on a faster timescale, the phasic influence of Ox on 5-HT activity is much larger.

Discussion
We have mapped out the direct and indirect connections between and within the DRN and LHA brain regions using IF staining to identify the receptors and trace long-range connections. To consider indirect connections, non-principal neurons have to be involved. In our study, the non-principal neurons are the inhibitory GABAergic and excitatory glutamatergic neurons, although there exist other neuronal types (e.g. neurons containing neuropeptide Y, etc) [87].
We have confirmed a previously identified 5-HT 1A receptor on Ox neurons, and have also identified multiple major 5-HT receptors in the LHA. They include, on Ox neurons, 5-HT 1B R, 5-HT 2A R, and 5-HT 2C R and also fast ligand-gated 5-HT 3A R. It is interesting to note that Muraki et al. (2004) has shown that the 5-HT 1A receptor antagonist WAY-100635 can completely block the 5-HT hyperpolarizing effect on Ox neurons. It could perhaps be that the 5-HT 1A receptors have the highest affinity as compared to the other 5-HT receptor subtypes found in this study. Further  experimental work would help to clarify this issue. In addition, we found both fast and slow 5-HT receptors (5-HT 3A R and 5-HT 1A R) on LHA's GABAergic neurons. Hence, these 5-HT (especially 1A) receptors could directly and indirectly influence Ox neurons, and further studies would be required to compare their relative affinities, and hence relative influences on Ox neurons. Within the LHA, we have also found the presence of OX 1 R, OX 2 R and NMDAR1 receptors on Ox neurons but not on GABAergic neurons. For LHA-to-DRN projections, we have identified Ox neurons receiving glutamatergic and GABAergic post-synaptic inputs in the DRN. The results from our experimental work and from previous work are summarized in Figure 10A. Our results also suggest that LHA's GABAergic neurons could be isolated from direct (excitatory) afferent influences from local glutamatergic and Ox neurons, but could be influenced directly by 5-HT neurons, or through a unidirectional closed loop that involves the Ox and 5-HT neurons ( Figure 10A).
From Figure 10B, the local LHA architecture looks similar in some respects to the well-studied cortical column architecture, with the Ox neurons acting as pyramidal neurons, with excitatory feedback among themselves [88]. However, our present work and in a previous study [4] have shown Ox neurons to have very weak influence on their local GABAergic neurons. This differs from the stereotypical excitatory-inhibitory feedback loop in a cortical column model. The DRN also looks similar to that of a cortical column if we include glutamatergic neurons (or glutamatecontaining 5-HT neurons) to provide excitatory feedback. But excitation in the DRN is known to be weaker than GABAmediated inhibition [59]. To provide a systemic understanding of the identified DRN-LHA architecture, we incorporated some of our current findings and previous published data into a computational model, which is an extension of our previous model [58].
We purposely simplified the neural unit and neuronal interaction implementations to focus primarily on this unique neural architecture of the DRN-LHA system. As a first step, glutamatergic neurons were justifiably omitted in our model simulations and analyses. We have included various constraints to the model, specifically on the relative values of the input-output slope, total input currents, and time constants and relative strengths of the connections. More importantly, we intentionally constructed a model to demonstrate the complex consequences of the neural circuit architecture we have established from our current experimental study.
We first use our model to explore the consequences of 5-HT's effect on LHA's GABAergic neurons. We found that the system becomes oscillatory when the connection strength is weak or inhibitory ( Figure 12). This oscillatory behaviour has so far not been observed in experiments. Thus, based on these results, we hypothesize that this connection is excitatory. It would be interesting to test the strength of this connection using 5-HT 3A R agonists/antagonists on LHA's GABAergic neurons. This hypothesis would also imply that DRN may send inhibition to Ox neurons both directly and indirectly (through the GABAergic neurons), i.e. no balanced projections. This is in contrast with the long-range projection of Ox to the DRN, which excites both 5-HT neurons and GABAergic neurons in the DRN ( Figure 10B).
Another interesting experimental finding in our work is the indication of an absence of Ox receptors on LHA's GABAergic neurons. Our model simulations show that this particular connection does not affect the coupled DRN-LHA system significantly ( Figure 13). This means that Ox receptor on LHA's GABAergic neurons will have little influence on the circuit's dynamics. It is interesting to speculate that the absence of these Ox receptors could be due to a lack of significant functional roles at the circuit level.
From our experimental work, we have identified both slow 5-HT G-protein-coupled and fast ligand-based receptors on both Ox and GABAergic neurons in the LHA. Our model attempted to mimic these different 5-HT receptor mediated timescales separately to investigate how the DRN-LHA circuit as a whole can be affected. We found that the 5-HT timescales do not change the tonic (steady-state) activities of the system, but can greatly affect the transient activations ( Figure 14). In general, a faster transient 5-HT influence on Ox neurons does not affect the suppression much but can result in a faster disinhibition of Ox neurons. This could mean that the faster 5-HT 3A R could be useful for quickly resetting the Ox neurons back to baseline after phasic 5-HT activation. More interestingly, a faster transient Ox influence can excite phasically 5-HT activity while slower timescale does not.
In summary, we have established aspects of the neurobiological circuitry function between the levels of 5-HT and Ox through direct and indirect pathways between the DRN and LHA. This work could have important implications in clinical neuroscience and neuropsychopharmacology as this DRN-LHA loop has been interpreted in two ways. It has been hypothesized that lower levels of 5-HT (common in depression) provide weaker inhibition to Ox neurons. Taking into account the effects of the circadian regulation of Ox and its influence on other neurotransmitters or neuromodulators, the Ox level may increase by this change in 5-HT. This increase in Ox levels may then generate a state similar to insomnia and other mood alterations. Similarly, in the reverse direction, it has been argued that lower levels of 5-HT are due to either a weaker excitatory connection from the Ox neurons (LHA) or because of lower levels of Ox (LHA), a situation commonly seen in hypersomnia or in narcolepsy [3,89]. To explore the potential effects of Ox knock-out mice on 5-HT activity, our model can simulate such an effect by removing all the Ox effects in the circuit, and we observed that 5-HT neurons can still fire at ,3.5 Hz (not shown). Furthermore, drug studies often do not consider integrating multiple targeted and non-targeted but connected brain areas. For example, in the DRN-LHA circuit considered in this study, an administration of 5-HT 1A R agonist can directly affect not only 5-HT autoreceptors, but also all the 5-HT 1A R in the DRN's GABAergic neurons, and LHA's Ox and GABAergic neurons. Other brain regions without 5-HT or Ox receptors, but connected to the affected DRN and LHA, will also be indirectly affected. Thus, the overall effect is complex, and this could be one important reason underlying serious side effects of various neuropharmacological drugs. A promising approach to gain a holistic understanding of such complex neurobiological systems is to perform more intensive computational modelling, simulations and analyses.