Control of sleep-to-wake transitions via fast amino acid and slow neuropeptide transmission

The locus coeruleus (LC) modulates cortical, subcortical, cerebellar, brainstem and spinal cord circuits and it expresses receptors for neuromodulators that operate on a time scale of several seconds. Evidence from anatomical, electrophysiological and optogenetic experiments has shown that LC neurons receive input from a group of neurons called hypocretin neurons that release a neuropeptide called hypocretin. It is less well known how these two groups of neurons can be coregulated using GABAergic (GABA standing for gamma aminobutyric acid) neurons. As the time scale for GABAA inhibition is several orders of magnitude faster than that for the hypocretin neuropeptide effect, we investigate the limits of circuit activity regulation using a realistic model of neurons. Our investigation shows that GABAA inhibition is insufficient to control the activity levels of the LCs. Although slower forms of GABAA can in principle work, there is not much plausibility due to the low probability of the presence of slow GABAA and lack of robust stability at the maximum firing frequencies. The best possible control mechanism predicted by our modeling analysis is the presence of inhibitory neuropeptides, which exert effects on a similar time scale to the hypocretin/orexin. Although the nature of these inhibitory neuropeptides has not been identified yet, it provides the most efficient mechanism in the modeling analysis. Finally, we present a reduced mean-field model that perfectly captures the dynamics and the phenomena generated by this circuit. This investigation shows that brain communication involving multiple time scales can be better controlled by employing orthogonal mechanisms of neural transmission to decrease interference between cognitive processes and hypothalamic functions.

The output of LC neurons is likely to be regulated by GABAergic cells in the sublaterodorsal peri-LC, providing a substantial input to LC cells [26,27]. This introduces a very intriguing interplay between two very discrepant time scales-involving HCRT neuropeptides with a decay time of minutes and GABA A receptors with decay times in the millisecond range. The interplay of these time scales, accommodating several orders of magnitude, is not infrequent in the hypothalamus [28]. These deep neural circuits have been associated with the control of REM sleep atonia [29,30] and bursts of LC overexcitation lead to (reversible) behavioral scenarios associated with neuropsychiatric disorders [21]. Therefore, imbalances in excitation/inhibition in the LC may underlie sleep disorders, indicating LC regulation to be an important issue that needs to be understood in depth. Other forms of corticothalamic feedback control have been studied before-as key gears beneath brain oscillations (e.g., spindles) during slow-wave sleep [31][32][33][34][35].
Here we investigate through realistic and reduced mean-field models the interplay between a fast amino acid (GABA) and a slow neuropeptide (HCRT) in the sleep-to-wake transition regulation. We have modeled the system as three groups: LC, HCRT and GABA populations. The main hypothesis is that the GABA inhibition controls the LC activity. LC and HCRT model cells from previous publications [11,16] were used, and the GABA model is introduced in section 2. We then show that, unexpectedly, the presence of a GABA population actually increases the activity of the LC cells. We use phase-resetting curves to provide an understanding of this actually being due to the active integration inherently present in LC cells. In order to suggest a possible regulation mechanism, we propose a new group which releases a neuropeptide acting on the same time scale as the HCRT peptide, but with an inhibitory effect on LC cells. This new population can balance the excitation that comes from HCRT cells and accurately regulates the activity of LC cells, not only decreasing its firing rate but also allowing a nimble rise of activity on the onset of HCRT excitation as well as a swift decay at the end of excitation. We also show a very simple mean-field model, based on simplified Wilfon-Cowan equations, that captures the dynamics of this circuits. We then discuss how our hypothesis can be tested experimentally.

The network model
In previous studies [11,16], we have used two conductance-based compartment models to describe LC and HCRT cells, as well as chemical rate equations for the release of neurotransmitters to simulate AMPA synapses. These models were fitted to agree with electrophysiological recordings of neural populations present in the hypothalamus and are described elsewhere [16]. Here we add a third group, probably located at the peri-LC, which inhibits both HCRT and LC cells through GABA A receptors (see figure 1). For the jth GABA cell, we use two compartments: the soma,   (1) and (2) are described in appendix A.
We model the GABAergic synapses similarly to the AMPA and HCRT connections used previously [16,36]. The release of GABA neurotransmitter by neuron j, which should bind at GABA A receptors in LC and HCRT cells, evolves according to the kinetic equation . Then the current going into a cell labeled k, in either the LC or HCRT populations, is defined as . To draw the connections, we set a probability of connection between groups and then choose pairs of neurons uniformly. Following observations of the existence of indirect GABAergic inhibition from the HCRT population [37], the GABA neurons receive excitatory inputs from HCRT and inhibit the LC neurons (without any feedback). The connection probabilities are shown in figure 1. In the simulations we used the backward differentiation formula and Newton iteration implemented within CVODE, with relative and absolute tolerances of 10 −9 [38]. We used 20 neurons per group, unless stated otherwise, resulting in a total of 845 coupled differential equations.
Carter and collaborators [21] have shown that the stimulation of LC cells during an interval of the order of some seconds results in the wake transition. Thus we have set a stimulation protocol: we inject into the HCRT neurons (soma) a current of magnitude I DC within a time interval of 10 s. Then we can assess how the LC responds to a sudden increase of activity in the HCRT cells, which is observed in experiments whenever a sleep-to-wake transition takes place. We show in figure 2 the activity generated by this protocol. As the HCRT neurons are excited (starting precisely at 70 s), we see an immediate increase in activity of LC and GABA neurons. As the stimulation stops acting, the LC neurons eventually stop firing, which is expected. We also show some activities of individual neurons from each of the three groups, where a considerable change is seen when this external current takes effect.

The role of the GABA population
It is likely that this GABA population regulates the LC activity; it is especially so because an activity overload may lead to behavioral attacks similar to episodes of neuropsychiatric disorders [21]. For this reason we wanted to understand exactly how the GABA population is affecting the LC (and HCRT) populations. We show in figure 3(a) a completely unexpected result: as the GABAergic synapses get stronger, the firing rate of LC cells increases instead of decreasing. We show HCRT and GABA inputs to an LC cell, which have completely different time scales, in figure 3(b). In figure 3(c) we show that, on indefinitely exciting the HCRT cells, the firing rates of all three populations eventually reach a steady level, which mainly depends on g GABA . The steady firing rates may also depend on other parameters (strength of the synapses between HCRT and LC groups, amplitude of the HCRT excitation, connectivity probabilities, etc), but we focus on the dependence on g GABA in order to investigate the control of the LC activity. Then, averaging these steady firing rates over time, we can see the LC firing rate dependence on g GABA (see figure 3(c), right panel). It is clear that both GABA and LC neurons increase their activities, showing that this effect occurs for a wide range of g GABA . Since HCRT receives feedback connections from GABA (see figure 1), we see a meager increase in the firing rate (≈1.0 Hz) of the GABA neurons as g GABA becomes stronger, due to a mild synchronization of HCRT neurons.
To help make sense of the g GABA values that we present in figure 3, we have checked, by direct inspection of the postsynaptic membrane potential, whether the amplitude of the inhibitory postsynaptic potential (IPSP) happened to be in an acceptable range (visible in figure 4(b)). For all g GABA values considered, the IPSP amplitudes are always less than 4 mV, which may be considered a reasonably strong inhibition in most systems.
This result is puzzling and opens the way to a few questions. For instance, if we simulate exactly the same system, but with a GABA population slightly larger, will the LC firing rate get smaller? We have tested this with a GABA population consisting of 30 and 40 neurons (i.e., The firing rate of all populations as g GABA is raised. The particular case without the activity of GABAergic neurons is precisely the same as that considered in [11]. (b) The activity and input current of a randomly chosen LC neuron. The time scales of HCRT and GABA are completely different. During HCRT stimulation, GABA current peaks have a smaller amplitude due to the relative refractory time (originated by the delayed K + current). Nevertheless, it is clear that the frequency of these pulses increases significantly, driven by the HCRT stimulation. (c) We show that if the current that excited HCRT cells remains for enough time, then the firing rates of all three populations achieve an equilibrium (on the left a simulation with = g 500 GABA nS is shown). Then we use these steady firing rates to show (right panel) the decrease of the LC firing rate as g GABA grows.
twice the LC and HCRT populations) and the result is that the LC firing rate can get even a little higher. Also, in general, spikes induced by inhibition are usually due to low-threshold currents, and LC neurons present both calcium I h and I T (T-type channel) low-threshold currents. On shutting both currents down (independently), this increase in firing rate with g GABA persists; thus it is not caused by the presence of such currents. Finally, this result does not substantially change if we modify the probability of connection from the LC group onto itself. The only restriction is that this probability must be higher than ≈0.2, below which exciting the LC group becomes much slower and does not comply with experimental observations. To some extent, these loopback connections help the LC group to engage some activity level quickly.
One last question is that of whether this whole effect is robust against changes in the decay time constant k 1 GABA r , defined in section 2. This time constant actually defines the time scale of the GABA A action on LC cells, and reports show that this constant can be as large as 8 ms. In fact, there are observations of two kinds of GABA A inhibitory postsynaptic currents [39]: GABA A,fast , which accounts for this range of the decay constant, and GABA A,slow , which manifests a decay constant larger than 30 ms and is present in all of the brain, though less frequently. We started investigating the GABA A,fast using several values in the range 1.5-8 ms, Figure 4. How LC neurons, whose parameters were fitted to behave as in previous experiments [16], behave when receiving GABA neurotransmitter. (a) Two cells, one LC and one GABA, are simulated. First we insert a weak current into the LC soma to make it spike intermittently. After a while (20 s), the GABA neuron receives a similar current and starts firing, inhibiting the LC cell. Although pulses of membrane hyperpolarization are observed in the LC voltage time series, the firing rate of the LC cell actually increases significantly. (b) Membrane potentials of both perturbed and unperturbed (reference) model LC neurons. Although an inhibitory pulse is usedwhich generates the hyperpolarization around 9 s-the next spikes get advanced with respect to the unperturbed case. (c) Phase-resetting curves (PRCs) of the LC cells for an inhibitory bi-exponential pulse, fitted from the output GABA currents, with several intensities I p . Notice that prompting such inhibitory perturbation over a considerable range of phases (for most cases, over 80% of the whole range) actually results in an advance of the phase (the same effect as is seen in (b)). This is the root cause for the increase in firing rate. and the same effect persists. Also, for a time constant of about 5-6 ms, the increase in firing rate becomes even more pronounced.
However, when we used a decay time scale consistent with GABA A,slow (namely, 30-100 ms), the GABA population was able to decrease the firing rate of the LC cells. Although at first sight this might appear to be a possible mechanism for LC activity regulation, there are features that suggest otherwise, such as the fact that most connections into LC neurons are gap junctions [40], ruling out the GABA A,slow known underlying mechanisms (diffusion and affinity of GABA A receptors) [39]. This makes this hypothesis very unlikely, but still possible. Also the observed inhibition in our simulations is extremely noisy, introducing non-observed oscillations in the LC activity (not shown), and thus it does not seem to be as good a control mechanism as the hypothesis that we shall present in section 5.

Inhibition-induced advances in the spike timing
GABA is the main inhibitory neurotransmitter present in mammalian central nervous systems, yet it is increasing the firing rate of LC neurons somehow. In figure 4(a), we isolate this effect by simulating only two neurons: one LC cell and one GABA cell. The GABA connects to the LC via equation (4) and they both receive a steady input current, starting at different times. First a 26 pA current is injected at the LC cell, making it fire intermittently at a constant rate of nearly 2.7 Hz (which is comparable to the behavior of these cells during awake states [21]); then, after 20 s, the GABA neuron receives a similar current and starts perturbing the LC neuron. One can see that the firing rate of the LC doubles at this moment, even though the voltage time series of the LC shows IPSPs (see 4(b), and more details below). To maximize this effect, we have used = × g 3 10 GABA 3 nS. We next use phase-resetting curves [41][42][43] to achieve an understanding of what happens when an LC cell receives an inhibitory presynaptic potential.
Phase-resetting curves (PRCs) are widely used to investigate synchronization and phaselocking issues, both theoretically [44] and experimentally [45,46]. The main idea is to assess, once an oscillator is trapped on a limit cycle, how it responds to an external stimulus prompted at different phases ϕ of its cycle. If after receiving this stimulus the oscillator remains in the basin of attraction of its limit cycle, then eventually it will return the original oscillation, but with a different phase, given the trajectory deviation. A PRC is the graph of this phase difference Δϕ against the phase ϕ at which the stimulus is prompted to the oscillator. Following our definition of a PRC [41], a positive Δϕ means an advance in spike timing-the next spike happens sooner; a negative Δϕ means a delay.
In figure 4(b) we show our protocol for assessing the PRCs. We injected a current into an LC cell to make it spike intermittently, with a steady frequency of approximately three spikes per second. The inhibitory pulse is prompted near 9.0 s. We fitted the GABA A receptor current from our model and used it as stimulation, with a scaling factor I p to make the pulse as strong as one wants. There is a direct correspondence between g GABA and I p , but they are not exactly the same, since I p is just a factor in the fitted bi-exponential formula 4 . Then, after a small transient (usually taking no more than one spike), the phase difference Δϕ is measured. This leads us to the PRCs that we show in figure 4(c), which are characteristics of the so-called type II oscillators (having positive and negative Δϕ). In the majority of cases, more than 80% of the 4  phase range has a positive Δϕ, which means that, although the inhibitory pulse causes a hyperpolarization in the cellʼs membrane voltage, the end result is an advance in the forthcoming spikes. Also, in figure 4(b) we can already see this phase advance effect, as after an IPSP the spikes of the perturbed time series arrive sooner than the reference time series spikes. And as this phase advance is often manifested in these cells, an inhibitory train of pulses may increase the firing rate, explaining figures 4(a) and 3. We believe that this is the effect that is driving our observations in section 3.

Regulation of LC activity: a hypothesis
LC activity cannot be disrupted into long bursts of activity, as this leads to sleep disorders, and it is believed that some neural population should exist to deter this undesired scenario. Through a GABA neurotransmitter alone, this is not accomplished. We have tested, for instance, what happens if we change the strength of the synapse of the HCRT group with the LC if either there is or there is not a GABA population present. In all cases tested, whenever the GABA population was present, the LC firing rate was boosted, and we see no regulation at all.
This follows from the inner characteristics of our LC neurons: they are actively integrating their inputs, with a fast rise after a hyperpolarization. Other neural dynamics might actually show a decrease in firing rate with intermittent GABA excitation, especially passive models (such as leaky integrate-and-fire models). Unfortunately, this is not the case with LC neurons, as shown in recordings in previous experiments [11]. In short, the PRC should present a considerable portion of phases with negative Δϕ, instead of positive. Except for an exotic dynamical system capable of both reproducing observed population features of LC recordings and also showing an adequate PRC, GABA neurotransmitter is not the answer for LC activity regulation.
We then have formulated the following hypothesis: what if the regulation is given by the competition between two neuropeptides with large time scales (comparable to that of HCRT), one excitatory and another inhibitory? This would generate a balanced input current in the time scale of tens of seconds, regulating the activity of LC cells. This may not be common, but in specific cases localized peptides can exert opposite effects on postsynaptic neurons, and one example is the case of the HCRT population itself, which receives both hypocretin as an excitatory neuromodulator and dynorphin as an inhibitory neuromodulator [47]. There are plenty of neuropeptides with synaptic currents that are very similar to those of the HCRT model we use here (see figure 3(b)) [48].
To answer this question, we add another population into our model which inhibits the LC population using an inhibitory neuropeptide, with time scales comparable to those for HCRT. We shall refer to this population as inhibitory neuropeptides (INPs for short). For simplicity, we shall use for the INPs exactly the same model as for the HCRT cells, but with an inverted sign of the synaptic strength. In this way, the net contribution to an LC neuron becomes a competition between HCRT excitation and this unknown peptide inhibition (shown in figure 5(c)).
Since the idea is that this neuropeptide acts whenever the LC population becomes overloaded with activity, a connection from the LC to this inhibitory neuropeptide makes sense. We have also connected into this population the HCRT, so that it can engage in some residual activity and respond to LCs more quickly. We show our main result in figure 5(a). This inhibitory neuropeptide is actually capable of regulating LC activity.
We also show for completeness, in figure 5(b), the activity of a candidate neuron per group, showing a significant increase in activity right after the LC cell starts firing. We can finally compare, in figure 5(d), how the firing rate of each population depends on the inhibitory synaptic strength g of the INP (as opposed to the case shown in figure 3(c)). We use the same protocol as was proposed before: we stimulate the HCRT population indefinitely and wait for all populations to reach their stationary firing rates. With this INP population, the LC firing rate decreases with g, thus suggesting that this could be a mechanism underlying the LC activity regulation. The scale of g (kS in figure 5(d)) matches the strength of the synapse from the HCRT to the LC used. There are no significant differences, from this point of view, between having or not having the GABA neurons, although they may be important in other ways. The decrease in the INP firing rate is due to a lower excitation coming from the LC cells, and thus a lower recruitment.
It is important to note that the INP population with a right choice of parameters-without a need for fine tuning-combines a lower firing rate with a swift rise of LC activity and a faster end of activity: comparing figures 5(b) and 3(a) we see that the activity of the LC cells dies out almost 10 s before, and that the highest firing rate achieved during the stimulation of HCRT cells is only slightly smaller. Although the onset of the LC activity is delayed with the INP population (by less than 3 s), the rising time scale is considerably faster. These features suggest the INP as providing a nimble regulation mechanism for LC activity. This observation may be the most important result and the most compelling signal that the existence of such a neurotransmitter could be of great importance to this hypothalamic circuit.
Also, this activity regulation does not only hold when the inhibitory and excitatory neuropeptides have precisely the same parameters, which would be a fairly unrealistic situation. We have found that the key feature beneath this effect is the use of an inhibitory neurotransmitter with a comparable time scale. To test this, we have changed the decay time constant for INP release over a wide range (from 10 −5 to 10 −3 m s −1 ) and checked that this population can still regulate LC activity.

The simplified network model
We now propose a mean-field-like model, based on simplified Wilson-Cowan equations, that can reproduce the firing rates of all three populations with their corresponding mechanisms of information transmission. We intend to provide a simplified tool that can be related to realistic simulations and be used for further analysis.
First, we propose the following system of coupled nonlinear differential equations for all populations: where F j (t) is the firing rate of population j (with 0 → HCRT, 1 → GABA, 2 → LC and 3 → INP), θ j being a population-specific gain function, A j a rate of excitation, γ j a damping term and F j0 the residual activity when there is no excitation. I(t) is the external current applied only to the HCRT group (i.e., δ = b j j ,0 ). We have used as the gain function To find the appropriate coefficients for each population separately, we use a least-squares minimization and neglect some small coefficients (see appendix B), ending up with the following equations:   . Additional values of the fitted coefficients are described in the appendix B.
For further testing our model, we compare the solutions of these mean-field equations with the findings from simulations, where the parameters differ from those used to fit them. For instance, on changing the amplitude of the current used in the protocol shown by figure 1, our simplified model can reproduce simulations within an error of 20%. We show in figure 6 a case in which we have fitted the model with a protocol close to that of figure 3(c), and then used it to predict the behavior of the system with a current square pulse as in figure 3(a). The errors may vary in the range 5% to 30%.
These simplified equations capture the trends of all populations with a reasonable agreement, giving a more intuitive insight into how these populations interact with each other.

Discussion
In this paper we have modeled a neural system composed of three populations of twocompartment conductance-based model neurons: the HCRT cells, composed of neurons expressing the HCRT neuropeptides, the LC cells, responsible for triggering the sleep-wake transition, and the GABA interneurons, which inhibit both HCRT and LC neurons through GABA A receptors. Although direct optogenetic stimulation of LC cells leads to wakefulness, vigilance state transitions are not exclusively managed by this neural group. LC neurons integrate two completely different time scales: one of the order of milliseconds, used by AMPA and GABA neurotransmitters, and another of the order of seconds, used by the HCRT neuropeptides. We have thoroughly analyzed the mechanism of LC activity regulation, which should be capable of suppressing overloads of activity that lead to behavioral states similar to episodes of neuropsychiatric disorder [21].
We have shown that due to GABA A neuroreceptors, the activity of LC cells actually increases. This unexpected effect is a result solely of how these cells respond to inhibitory pulses. An analysis of their phase-resetting curves (PRCs) shows that 80% of the phase range leads to an advance in the upcoming spike times, which thus leads to an increase in firing rate. This is due to the fact that the LC cell model integrates its inputs actively, in contrast to passive models such as the leaky integrate-and-fire model. Thus a fast hyperpolarization pulse consequently provokes a fast rising in the membrane potential and an advance in spike timing.
This effect of increased activity induced by GABAergic neurons has been previously observed, although through different mechanisms [49]. Also an inhibitory synapse that causes Figure 6. The mean-field equations (6)-(9) after parameters were fitted using a leastsquares minimization (see appendix B for details). Bright lines are simulation results with the same parameters for comparison. In this particular case, we have fitted the mean-field model using a steady current applied from t = 60 s onwards (like for figure 3(c)). This is shown in the left panel. For the right panel, we have used the parameters from this fitting and used the model with a square-pulse current (the same as that used in figure 3(a)). This shows that the mean-field equations capture the trends of all populations and can be used to give a simplified view of this hypothalamic network. spike timing advances is not excessively rare, and we can cite for instance bursting neurons in central pattern generators [50].
Since a fast inhibition would not consistently regulate LC activity, we propose the existence of a population that releases a neuropeptide similar to hypocretin/orexin but with an inhibitory effect-i.e., working on a slower time scale. A simplified view of the new network configuration is depicted in figure 7. This population, which we refer to as INPs, could use for instance any known inhibitory neuropeptide with a long enough time scale. We have shown that with this new model the LC activity can be regulated, especially if there is a feedback connectivity: not only does the INP inhibit the LC cells, but also the LC excites the INP. We have proposed this connectivity since a large activity coming from the LC population would immediately evoke a larger inhibition by INP neurons. Carefully choosing the parameters, we have shown that this can actually lead to a nimble rise of LC activity and, when the stimulation of the HCRT is over, a similarly swift end of activity. This is probably the most compelling argument for the existence of a neuropeptide with such characteristics perhaps being of great importance to this hypothalamic circuit. Among the possible INPs acting on LC neurons are MCH [51][52][53] and opioid [54,55] peptides. Finally we proposed a mean-field model that captures the main elements in the interaction among those groups and can be used as a simplified view of this hypothalamic network for further analysis.
We have also reported that a slower GABAergic inhibition, similar to what is usually called GABA A,slow currents, may also decrease the LC activity. However, this kind of current is unlikely to arise in LC cells and the regulation is less stable than the results achieved by this hypothetical INP population.
The interplay of two neuropeptides balancing each other has been studied before as well, though from different perspectives [47,48,56]. Experimental tests can be performed on this hypothalamic circuit to find out whether our predictions are borne out-in particular, starting with active stimulation (by current clamping, for instance) of these GABAergic cells and studying the effect on the LC population; this should induce an increase of the firing rate. This empirical observation, regardless of the end result, may provide new ingredients helping to provide an understanding of the role of each of these three populations and how LC activity is regulated. Figure 7. The new configuration proposed for active control of the LC activity. We have introduced a population which releases a slow neuropeptide inhibiting LC group. We propose a feedback connection from the LC to the INP for a more robust control mechanism. Among possible INPs acting on LC neurons, we can cite MCH (melaninconcentrating hormone) and opioid peptides.

Appendix A. The GABAergic interneuron model details
We provide here details additional to those provided in section 2 of the currents used in equations (1) and (2). The time dependence and units are intentionally dropped to simplify the notation whenever they are clear from the context.
We start with the currents involved in the axon compartment. First we consider the sodium current, modeled as   where we have set the calcium dynamic dissipation as μ = 1.5, and I Ca is the calcium current, which takes effect in the soma compartment, described next. The calcium current was modeled using the Goldman-Hodgkin-Katz formalism [58], as = − ( ) Equation (5) defines a class of models from which we want to find one that reproduces our numerical results. The idea then is to find a fit from the data A j , a jk , β j , σ j and γ j for all ∈ j k , {0, 1, 2, 3}. Since the F j0 are easily accessible from simulations, and in order to decrease the complexity of the fitting process, we have neglected them and assumed the values given below. We use a simple procedure of minimization of a least-squares objective function. First, for clarity we shall use the following notation: β σ γ ( ) 3 , to explicitly express the dependence of the parameters through the mean-field equation (5). Suppose  j is the firing rate time series of population j, as for instance shown in figure 2. In order to make the convergence of the parameters during the minimization easier, we have at this point used a Gaussian smoothing algorithm on these time series. Then we can estimate the derivative by calculating , which then is found by using a Broyden-Fletcher-Goldfarb-Shanno algorithm.
However, to actually use β σ γ F t A a ( ; , , , , ) j j j j j j in equation (B.2) we need the firing rates of other populations. This is an omitted dependence and we have used, only during the minimization process, the time series  j to meet this need.
Given the number of parameters, we have removed some of them in order to make the analogy to figure 1 as clear as possible. For instance, a 02 was set to zero because there is no direct contribution from the LC group to HCRT cells. Also other parameters were removed from the final mean-field equations (6)-(9) because they were found to be too small, thus only adding more noise and iterations to the algorithm. This is the case for a 01 , representing the feedback from GABA to HCRT.
Finally, we provide the parameters fitted for . Notice that F j0 agrees with the residual firing rate for each population.