Diversity and Noise Effects in a Model of Homeostatic Regulation of the Sleep-Wake Cycle

Recent advances in sleep neurobiology have allowed development of physiologically based mathematical models of sleep regulation that account for the neuronal dynamics responsible for the regulation of sleep-wake cycles and allow detailed examination of the underlying mechanisms. Neuronal systems in general, and those involved in sleep regulation in particular, are noisy and heterogeneous by their nature. It has been shown in various systems that certain levels of noise and diversity can significantly improve signal encoding. However, these phenomena, especially the effects of diversity, are rarely considered in the models of sleep regulation. The present paper is focused on a neuron-based physiologically motivated model of sleep-wake cycles that proposes a novel mechanism of the homeostatic regulation of sleep based on the dynamics of a wake-promoting neuropeptide orexin. Here this model is generalized by the introduction of intrinsic diversity and noise in the orexin-producing neurons, in order to study the effect of their presence on the sleep-wake cycle. A simple quantitative measure of the quality of a sleep-wake cycle is introduced and used to systematically study the generalized model for different levels of noise and diversity. The model is shown to exhibit a clear diversity-induced resonance: that is, the best wake-sleep cycle turns out to correspond to an intermediate level of diversity at the synapses of the orexin-producing neurons. On the other hand, only a mild evidence of stochastic resonance is found, when the level of noise is varied. These results show that disorder, especially in the form of quenched diversity, can be a key-element for an efficient or optimal functioning of the homeostatic regulation of the sleep-wake cycle. Furthermore, this study provides an example of a constructive role of diversity in a neuronal system that can be extended beyond the system studied here.


Introduction
Disorder, which originates from both noise and diversity, is naturally present in all biological systems. In neuronal systems some examples are the random opening and closing of ion channels, the multitude of stochastic input currents in the neurons, and the diversity of shapes, sizes, and electrophysiological properties of the neurons [1,2]. Disorder is often considered to be harmful to the systems' functioning and to information encoding. However, it was likewise repeatedly demonstrated that a certain level of disorder can facilitate signal encoding by enhancing system's response to an external stimuli. For instance, quenched diversity clearly shows its constructive role in the phenomenon of diversity-induced resonance, in which an assembly of heterogeneous excitable units presents an optimal response to an external forcing for a suitable intermediate degree of heterogeneity [3,4,5]. Similar constructive effects can be observed in the presence of noise. For example, interplay of noise and nonlinear forces produces the directed motion of motor proteins [6], order-disorder transitions, oscillations, and synchronization in assemblies of excitable units [7,8,9], and an optimized system response in the ubiquitous phenomenon of stochastic resonance [10,11], e.g. in ion-channels and neurons [12,13,14,15,16,17].
In the present study we examine the effects of noise and diversity (heterogeneity) in a physiologically based neuronal model of sleep-wake cycles [18]. This model introduces a novel mechanism of the homeostatic regulation of sleep based on the dynamics of a wake-promoting neuropeptide orexin (also called hypocretin), assuming depression of orexinergic synapses during wakefulness and their recovery during sleep. This mechanism is based on the experimental findings of the essential role of orexin system in maintaining wakefulness and its ability to integrate the sleep-wake relevant information coming from many brain areas [19,20] and respond to changes in the body external and internal environments by encoding the body activity state, energy balance, sensory and emotional stimuli [21,22].
In the original model interaction between only two representative neurons is simulated: the orexin neuron and the local glutamate neuron that are reciprocally connected to each other according to the experimentally established physiological connec-tions [23]. Both orexin and glutamate neurons are firing during wakefulness and are silent during sleep. The transitions between firing and silence are governed by the interplay between the circadian input and homeostatic mechanisms as initially proposed by Borbely [24]. For simplicity, in this model only a single type of orexin neurotransmitter (instead of the two types actually known) is considered, and it is assumed that the system can be either in the wake state or in a generic non-Rapid Eye Movement sleep state, without specifying ultradian structure of sleep. Also this model did not consider noise effects, and diversity could not be included since there are only two neurons present.
In the present paper we extend the above described two-neuron model to a more realistic multi-unit model with heterogeneous neurons. The aim of the study is to first of all investigate how the presence of diversity in the neuronal population affects sleep-wake transitions, since it is well-known that neurons are highly heterogeneous by their nature. In particular, within the orexin neurons population significant intrinsic diversity can be found: different electrophysiological properties, sizes in the diameter range 15-40 mm, and various shapes such as a spherical, fusiform, or multipolar [19,25,22]. Secondly, also stochastic fluctuations, representing current noise, are added to the model and the response of the system is studied for different levels of noise. The question naturally arises, to what extent noise and diversity are essential ingredients for the functioning of assemblies of neurons and other complex systems, and what is the optimal level of noise and diversity required for the emergence of an optimal response to external stimuli. It is shown below that the model under study presents both diversity-induced resonance and stochastic resonance, but the former appears more clear and robust, since it is always associated with a regular almost-periodic spiking-silence activity, rather than to the irregular random transitions characterizing the stochastic resonance regime.

Materials and Methods
In this section the two-neuron model of sleep-wake cycles [18] is described and some examples of dynamics in the presence of an external periodic signal are illustrated. Further, this model is extended to account for multiple neurons dynamics and heterogeneity, and a simple quantitative criterion to estimate the quality of a sleep-wake cycle is introduced. This criterion will be used in the Results section to compare sleep-wake cycles dynamics obtained at different parameter sets.
The two-neuron model The original model of the homeostatic regulation of sleep has a minimal structure consisting of two representative interacting neurons A and B, as depicted in Fig. 1. The neuron A simulates a representative neuron from the orexinergic neuronal population, while the neuron B represents a local glutamate interneuron (for details see [18]). The state of wakefulness or sleep is determined by the firing regime of neurons A and B, since these neurons are known to fire during wakefulness and be almost silent during sleep (see e.g. [22]).
Interaction between the neurons A and B takes place through glutamate and orexin neurotransmitters, as detailed below. The neuron A is acted upon by a stimulus in pace with the circadian rhythm, here treated as a periodic external signal -a simplification justified by its independence from the homeostatic process [26]. The homeostatic process itself is described by an additional macroscopic variable M(t) simulating availability of orexin.
Dynamics of the neurons A and B are based on a Hodgkin-Huxley-type model [27]. The membrane potentials of the neurons A (V A (t)) and B (V B (t)) are thus calculated as: where C p (p~A,B) are the membrane capacitances per unit area of the respective neurons, I a (a~L,Na,K,gl,ox) are the ionic currents, g a are the maximum conductances, and E a are the equilibrium potentials. The capacitance values are taken as C A~CB~1 mF=cm 2 . The values of all the other model parameters are listed in Table 1.
In the following we give a detailed explanation of different parts of the model. N External forces. The current I ext acting on the neuron A and the noise currents j p (t), p~A,B, can be considered as external forces, in the sense that they do not depend on the system variables.
The external current I ext (t) is assumed to simulate a stimulus associated with the circadian rhythm. For simplicity in the present study a periodic pulse input is used to introduce circadian activation of the system: t, I ext (t)~I ext (tzt). Such current can be interpreted as an awakening effect of an alarm clock or some other disturbance coming with a period of 24 hours. In the following we employ a train of rectangular pulses with length t 0 (t 0 vt) and height I 0 , as depicted in Fig. 2 where n is an integer. This simple form is chosen because it is convenient for carrying out a systematic study of the neuron response at different parameters sets. However, it should be kept in mind that it represents a drastic simplification, and more realistic shapes of circadian currents can also be used [18].

Author Summary
All biological systems are inherently noisy and heterogeneous. Disorder is mostly expected to disturb proper functioning of a system, like it can be the case with noise in a radio signal. However, it has been demonstrated by numerous studies that noise can actually improve signal encoding -the so-called stochastic resonance phenomenon.
Recently, it was discovered that quenched diversity (heterogeneity) can also enhance the response of a system to an external perturbation (diversity-induced resonance). In this study we investigate the role of noise and diversity in a neuronal model of sleep-wake cycles based on the dynamics of the wake-promoting orexin neurons that is crucial for stability of wake and sleep states. We demonstrate that suitable levels of diversity introduced in the orexin neurons can significantly improve the quality of the sleep-wake cycle, and may be essential for proper sleep-wake periodicity. Noise, on the other hand, provides only a mild improvement.
The noise term j p (t) represents fluctuating currents that are known to be always present in neurons. For simplicity, we assume zero-average Gaussian white-noise processes: with D p being the noise intensity.
N Internal dynamics. The leakage, sodium, and potassium currents I p,a (p~A,B; a~L,Na,K) in the equation of the neuron p depend only on the variables of the same neuron p and, thus, describe the neuronal internal dynamics.
The leakage currents I p,L~gL (V p {E L ) represent a flow of ions with a small conductance g L &0:1 mS=cm 2 driving the membrane potential toward the negative value E L &{60 mV.
The depolarizing Na-currents I p,Na~gNa (V p {E Na )a p,Na have a maximum conductance g Na~3 mS=cm 2 and a large positive equilibrium potential E Na~5 0 mV. The activation variables a p,Na (t), with 0ƒa p,Na ƒ1, represent the fraction of open ion- Figure 1. Scheme of the two-neuron model of the sleep-wake cycle [18]. The A?B red arrow from the orexin-producing neuron A (red circle) to the neuron B (blue circle) represents the glutamate projection as well as the orexin projection regulating the homeostatic process. The blue arrow represents the B?A glutamate projection. The neuron A is also acted upon by a periodic signal representing the effect of the circadian clock. doi:10.1371/journal.pcbi.1002650.g001 Table 1. Parameters of the two-neuron model [18]. channels contributing to the Na current. Because of their fast activation relative to the other time scales, the Na-current is assumed to be activated instantaneously, according to its voltagedependency: where W(x) is the sigmoid function S Na is the steepness of the sigmoid function and W Na is the halfactivation potential. The repolarizing K-currents I p,K~gK (V p {E K )a p,K are characterized by a maximum conductance g K~4 mS=cm 2 , a large negative equilibrium potential E K~{ 90 mV, and a longer activation time than the depolarizing Na-current, namely t K~2 ms. Consequently, the dynamics of the K-currents activation variables are modelled as where W(x) is defined in Eq. (6).
Couplings. The neurons A and B are mutually coupled by chemical synapses through the glutamate-induced (I p,gl ) and the orexin-induced (I ox ) currents. Unlike the Na and K currents, I p,gl and I ox depend on the activity of both presynaptic and postsynaptic neurons. The activation variables a p,gl and a ox depend on the appearance of a spike in the presynaptic neuron, i.e. on the presynaptic voltage. Additionally these currents depend on the voltage of the postsynaptic neuron, similarly to other ionic currents. Both glutamate and orexin are excitatory neurotransmitters, so they are assumed to open depolarizing ion channels, such as Na-channels.
The activations of the glutamate-induced currents are modeled as: This equation is similar to Eq. (7) but has the important difference that the equilibrium value W(S gl (V p p {W gl )) for the activation variable a p,gl depends on the membrane potential V p p of the other neuron p p ( p p~B if p~A, p p~A if p~B). The time constant t gl~3 0 ms accounts for the delay coming from the activation of glutamate receptors, and the following activation of ion channels.
The orexin-induced current represents the effect of orexin produced by the neuron A and acting on the neuron B. It is modeled in a form similar to the glutamate-induced current. This current provides a simplified description of the effects of orexin on the neuron B which appear after a complex series of processes, involving production of orexin in the soma of the neurons, its release in the synaptic cleft, and activation of G-protein coupled metabotropic receptors. The dynamics of the activation variable a ox (t) depend not only on the membrane potential V A (t), but are also related to the availability of orexin at time t, described by the additional variable M(t) (0ƒMƒ1). The dynamics of the variables a ox (t) and M(t) are defined by the equations: dM dt~{ The term M|W(S ox (V A {W ox )) in the Eq. (9) reflects activation of the synaptic current due to appearance of a spike in the presynaptic neuron A. At the same time it determines the rate of orexin availability reduction in Eq. (10) due to spiking of the neuron A with a time constant t { ox . The first term in Eq. (10) determines recovery rate of the orexin availability with time constant t z ox .
The meaning of the product M(t)|W(S ox (V A {W ox )) is that there is orexin-induced activity in neuron B if (1) there is enough orexin available above a critical threshold [M(t)wM critical ], and (2) the neuron A is in the firing state [W(t)&1].
The time constants t + ox accounting for the orexin dynamics are much longer than the time constants associated with ionic current terms. The time constant t ox of the homeostatic regulation process is even longer, being of the order of magnitude of the daily period t.
For numerical convenience, simulations are made over rescaled daily and orexin time scales: the daily period was assumed to be t~24 s, instead of t~24 h, achieved through a suitable rescaling, which was applied to the orexin time scale t ox and the production and reduction times t + ox . The other time parameters are left unchanged. Since such rescaled t ox and t + ox are still much larger than any other time scale of the microscopic dynamics, the rescaling does not change the main results of the simulations. See [18] for a detailed validation of such rescaling procedure.
All the parameter values for the currents are listed in Table 1. It is assumed that the neurons A and B share the same parameter values, unless specified otherwise. Such an assumption is justified, because the major properties of these neurons required for the model are the tonic firing (periodic single spike activity) and silent states. Without any external input both neurons should be in a silent state, while they are brought to firing activity in response to depolarization. Therefore, change of parameters in a physiologically allowed range would primarily lead to the different amount of depolarization needed to excite neurons, and would not affect the major outcomes of the simulations.
The system defined above is essentially an excitable feedback system, i.e. both the external input of sufficient strength and the AB coupling are essential elements for maintaining firing activity of the neurons. Orexin-related dynamics, with the associated long time scales, are expected to direct the homeostatic sleep process, which regulates the sleep-wake transitions. The healthy sleep-wake cycles in this system are realized as follows: N Initiating wakefulness. A sufficiently strong or long external signal or a stimulus associated to the circadian rhythm, e.g. the idealized rectangular pulse considered here, activates the system and induces firing activity in the neuron A. Due to the excitatory synaptic connection from A to B, the neuron B is also activated. N Maintaining wakefulness. Once the pulse is finished and the external current is zero, the system remains in the wake state (i.e. both neurons A and B are firing) due to reciprocal excitation between the neurons. The firing activity lasts for a fraction q of the daily period t. Ideally one can assume q~2=3, corresponding to 16 hours for a day of 24 hours, i.e. 16 seconds for the daily period t~24 s with the time scales of the model considered here.
N Initiating sleep. The firing stops in both neurons due to decreased availability of orexin according to the dynamics of M(t). This is associated witch the transition from wake to sleep.
Two examples of the two-neuron model dynamics without noise are illustrated in Fig. 2. The left part of the figure represents the response obtained for a pulse length t 0~5 00 ms and height I 0~0 :895 mA=cm 2 . In each period orexin is depleted during the neuronal activity and recovered while the neurons are silent. The stimulus parameters used in this example have been intentionally chosen close to the critical firing threshold, so that by slightly reducing the pulse height or length, the periodic appearance of a continuous time interval of spiking regime is lost. Such case is demonstrated in the right hand side of the figure, where the current pulse height is slightly lower, I 0~0 :893 mA, while all the other parameters are kept the same. There the prolonged wake state is induced only every other day, because the input is insufficient to induce sustained spiking at the same levels of orexin availability. By reducing the pulse amplitude or duration even further it is possible to observe different behaviors such as triple or higher-order periodicities.

The heterogeneous model
As a step toward a more realistic model we generalize the twoneuron model into a heterogeneous multi-neuron model. For simplicity we first increase the number of orexin neurons only. To do this we replace the single neuron A by a set of N A neurons fA i g (i~1, . . . ,N A ), while still maintaining only one neuron B. Also, in this paper we assume that the diversity is constant in time in order to consider the simplest case possible.
In reality a certain level of heterogeneity is observed in all neuronal parameters. However, given that our model neurons are simple pacemaking neurons such diversification of different model parameters (in a physiologically allowed range) would simply lead to slightly different firing rates of the neurons. This, in turn, will result in diversity in activations of synaptic currents, which can be mimicked by simply diversifying their activation thresholds. Thus, in the following we can limit ourselves to studying the effects of diversity in activation thresholds of synaptic currents without loss of generality. Furthermore, as a first step, the heterogeneity is only introduced in the glutamate-induced currents to avoid having a too complicated system, which would become difficult to understand.
With regard to the coupling topology among the orexin neurons, so far there is no detailed experimental data. Therefore, for simplicity, we chose an all-to-all coupling via gap junctions, but other variations can be tested in the future. The intensity of the coupling has been chosen large enough to ensure that the neurons A i respond in pace to the external current. The equations of the two-neuron model are modified accordingly.
N Dynamics of the neurons A i . The membrane potentials V (i) A (t) of the neurons A i , are described by equations analogous to Eq. (1): The current terms are similar to those in the two-neuron model, apart from the additional coupling currents between two generic neurons A i and A j , , with i,j~1, . . . ,N A , where k int is the gap junctions conductance that can be treated as coupling strength. The currents' activation variables a (i) A,Na and a (i) A,K are modeled in accord with the equations of the two-neuron model. Note that the specific values of the activation variables will be different for different neurons since they depend on voltages of each particular neuron A i .
For simplicity the same external current I ext (t) given by Eq. (3) is assumed to act on all neurons A i (see Fig. 3). The noise terms j (i) A (t) as well as the noise j B (t) acting on the neuron B (see below) are also defined similarly and assumed to be statistically independent from each other. For convenience the properties of all stochastic forces are written together (i,j~1, . . . ,N A ): N Connections from the neuron B to the neurons A i . The neuron B has glutamatergic synaptic inputs to each of the neurons A i as depicted in Fig. 3. Diversity is introduced in the activation thresholds of the glutamate-induced currents according to the following equation for the activation variables: The thresholds W (i) B,gl adopt different values for each neuron A i that are independently extracted from a probability distribution defined later in the text.
N Connections from the neurons A i to the neuron B. Each of the neurons A i has synaptic projections to the neuron B. This is translated in the model by replacing the single glutamate-and orexin-induced currents with their averages such that Eqs. (8) and (9) for the activation variables become: Note that diversity is again introduced in the activation thresholds of the glutamate-induced currents W (i) A,gl corresponding to heterogeneous (A i ?B) synapses located at the neuron B. Due to the differences in the A i neurons, the orexin availability function M (i) (t) is different for different neurons, although still following Eq. (10).
The above described set of equations constitutes the multineuron heterogeneous model of the homeostatic regulation of sleep. Numerical results were obtained using a variation of the Runge-Kutta 2nd-order method, which is suitable for equations with stochastic terms, namely the Heun method [28]. Identical initial conditions were assumed for all neurons, corresponding to a silent state.

Quantifying the quality of the sleep-wake cycle
In this section a heuristic criterion is introduced in order to evaluate and compare the quality of the system responses obtained for different external signals or internal parameter values.
For this purpose, the period t is divided into a ''day'' wakefulness sub-period of length t 1~q t and a ''night'' sleep sub-period of length t 2~( 1{q)t, with t~t 1 zt 2 . The quantity q is defined as a wake fraction. A typical sleep-wake cycle with an eight-hour sleep sub-period has q~2=3. For the day corresponding to the n-th period (nt,(nz1)t), the ''day'' is represented by the sub-interval (nt,ntzt 1 )~(nt,(nzq)t), which covers the first fraction q of the period, while the ''night'' extends in the complementary fraction (1{q) of the period in the time interval (ntzt 2 ,(nz1)t)~((nzq)t,(nz1)t).
For each period n~0,1, . . ., we compute wakefulness time intervals Dt (1) n and Dt (2) n spent by the system in the wake state during the day, Dt (1) n , and night, Dt (2) n . The wake/sleep state is identified with the spiking/silent regime. A simple quantitative estimate of the quality of the sleep-wake cycle can, thus, be done through the following linear function of the wakefulness time intervals, Figure 3. Scheme of the heterogeneous model. Example of model system with N A~5 orexin-producing neurons fA i g, i~1, . . . N A (red spheres) and one neuron B (blue sphere). The neurons A interact with each other through an all-to-all coupling (red lines). Blue and red projections have a meaning similar to those of Fig. 1: the neuron B is coupled to the neurons A i through parallel glutamate projections, while each neuron A i is coupled to neuron B through a glutamate and an orexin projection. The neurons A are also acted upon by a stimulus representing the effect of the circadian clock (gray arrows). doi:10.1371/journal.pcbi.1002650.g003 r(Dt (1) ,Dt (2) where Dt (a)~P n Dt (a) n =N sp , a~1,2, represent the average of the wakefulness time intervals during the day (a~1) and during the night (a~2), with N sp being the total number of periods of the simulation.
The fractions Dt (a) =t a (a~1,2) can vary in the interval (0,1); then the coefficient r in Eq. (16) is limited in the interval ({1,1). The maximum value r~1 corresponds to an optimal cycle with Dt (1)~t 1 (wakefulness during the entire day) and Dt (2)~0 (sleep during the entire night); any deviation from the optimal state (r~1) comes either from values Dt (1) vt 1 (implying some sleep during the day) or values Dt (2) w0 (meaning at least some wakefulness in the night). See Text S1 for further details on the definition of the time intervals Dt (1) , Dt (2) and the coefficient r.

Results
In this section we study how the presence of disorder affects the system response and discuss the main differences compared to the two-neuron model. The term ''disorder'' is used here to refer to either noise, i.e. disorder in time (stochastic terms in the external current), or diversity, i.e a quenched heterogeneity in the neuronal parameters. These two aspects are studied separately. For the sake of simplicity, we examine the response of the system to a periodic stimulus represented by a train of short rectangular pulses as defined in Materials and Methods.
In each of the examples considered, the initial configuration in the absence of noise and diversity is the same as the sub-threshold state illustrated in Fig. 2-right with a double-periodic response. It is obtained for a reduced height of the current pulse I 0~0 :893 mA, while the other parameters are unchanged as given in Table 1. The reason for starting from such an under-threshold non-optimal configuration is that it is most sensitive and, thus, best illustrates the effects of added noise or heterogeneity. While a response with a double-periodicity may seem unrealistic, this starting configuration is intended to be an example of non-optimal response rather than a standard reference state. In fact, in realistic situations noise and heterogeneity are always present so that such a state without noise or diversity represents a hypothetical system that would be obtained if one could switch off noise or replace heterogeneous synapses with perfectly identical ones. The results presented below suggest that a multi-periodic sleep-wake cycle can be turned into a regular (singleperiodic) one by adding a suitable degree of disorder.

Effect of noise
Here we investigate the effects of the noise currents in the equations for the membrane potentials. For clarity only the cases in which noise currents are present either in the neurons A i or in the neuron B are considered.
Noise in the neurons A i . To study the effects of the noise currents j (i) A (t), i~1, . . . ,N A acting only on the neurons A i (as per Eq. (11)) we set D B~0 . Also, no diversity in the characteristic parameters of neurons A i is introduced. We have simulated a system with N A~2 0 identical neurons and a single B neuron on a time interval t[(0,100t). A raster plot for the activity of the neuron  Fig. 4-A. The plot shows that N for small D A &0 the system's configuration corresponds to the assumed non-optimal double-periodic solution; N the system's response becomes slightly more regular and periodic as D A is increased, despite the fact that the neuron cannot initiate a firing event at the beginning of each period; N as D A becomes even larger the neuron B keeps firing tonically for a longer and longer time interval (even longer than a single period) thus deteriorating the general quality of the response.
A sample of time dependence of the main variables in the interval (0,4t) for D A~1 mV and D A~2 mV is illustrated in Fig. 5. In general, the type of variability induced by noise currents acting on the neurons A i affects both the firing initiation and, especially, its duration. However, it is difficult to establish an actual improvement of the quality of such a response as a function of the noise intensity D A , as even the coefficient r, shown in Fig. 4  N at larger values of D B the state of sleep is frequently interrupted by almost isolated spikes at random times.
A representative example of time dependence of selected variables of the neurons A 1 and B are shown in Fig. 7. Note the different type of behavior induced by a high levels of noise acting on the neuron B, compared to the case in which noise acts on the neurons A i . In the first case irregular switching between the firing and silent states is observed more often, especially considering the transient firings in the otherwise silent sleep state. Furthermore, this random firing appears only in the neuron B, but is insufficient to also induce spiking in the A 1 neuron. This activity may represents intermittent awakenings, which are likely due to the ability of noise to favor the ignition of spiking events. Such random spikes are not observed when noise acts on the neurons A i only, even at much larger noise intensities. This may be related to the coupling between the neurons A i , which constrains them in the same (spiking or silent) state. In order to excite all neurons A i together one would need an input signal affecting all of them in the same way, which is highly improbable in a realistic system.
The dependence of the coefficient r on D B is shown in Fig. 6-B. Again, only a mild stochastic resonance behavior is suggested by the data when varying the noise intensity. It should be noted that in this particular configuration with noise acting only on the neuron B, the response of the neuron B does not depend on the number N A of homogeneous neurons fA i g, due to the equivalence to the configuration of the two-neuron model, as we have checked numerically. Thus, the plots of neuron A 1 in Fig. 7 are representative of all other neurons A i . In fact, the external current I ext (t) as well as the coupling currents are the same for each neuron A i , which produces the same response. According to the equations of the heterogeneous model, the effective current acting on the neuron B is the arithmetic average of the currents coming from the various neurons fA i g and, therefore, coincides with that of any single neuron A i . We use here a homogeneous multineuron generalized model only for a better comparison and consistency with the rest of this study.

Effects of heterogeneity
The effects introduced by a heterogeneity in the neurons are dramatic compared to the effects of noise. The corresponding improvement of the system response for suitable intermediate amounts of diversity can be detected very clearly. This is the main result of this paper and it is illustrated in this section. Noiseless neurons are assumed for easier estimation of the heterogeneity effects (D A~DB~0 ).
As in the study of noise described above, we carry out the study of diversity starting from the same configuration with a nonoptimal double-periodic response to the external periodic stimulus, corresponding to a zero diversity (homogeneous system). Heterogeneity is then introduced in the glutamate-induced currents, either in the thresholds W (i) A,gl regulating the response of the A i ?B synapses at the neuron B or in the thresholds W (i) B,gl of the B?A i synapses at the neurons A. This is done by randomly extracting values W from a probability density f p (W ) and assigning them to the threshold parameters W (i) p,gl (p~A,B). The probability density used here has a bell-shape f p (W )~P((W {W p,gl )=dW p,gl ), where P(x)!1=cosh(x) 2 , the quantity W p,gl~S W T represents the average value, while dW p,gl measures the dispersion of the distribution f p (W ) around the average value and is related to the standard deviation s p by dW p,gl~p s p = ffiffiffiffiffi 12 p . For further details see Text S1. The width dW p,gl is assumed in the following as the measure of neuronal diversity. In order to carry out meaningful comparisons with the homogeneous (two-neuron) model, the average values are set equal to the corresponding parameters of the homogeneous two-neuron model, The other parameters are unchanged compared to the two-neuron model, see Table 1. Diversity in the B?A i synapses (neurons A i ). Diversifying the potential thresholds W (i) B,gl implies heterogeneous glutamate synapses located at the neurons A i , see Eq. (13) and Fig. 3. That is, each neuron A i responds in a different way to the stimulation from the neuron B. Notice that this is a truly heterogeneous system which cannot be reduced to an effective two-neuron model-as in the case of heterogeneous synapses at neuron B considered in the next section. We studied a system with N A~2 0 neurons A i with diversified threshold parameters W (i) B,gl , i~1, . . . ,N A . The system dynamics were examined for different sets of thresholds fW (i) B,gl g extracted from distributions f B (W ) with different widths dW B,gl but always the same average value W B,gl~WB,gl .
The resulting raster plots of the activity of the neuron B are shown in Fig. 8-A, and a sample of time dependencies for the neurons A 1 and B is shown in Fig. 9. The existence of an optimal degree of diversity, corresponding to a value dW B,gl approximately Diversity and Noise in a Sleep-Wake Cycle Model PLOS Computational Biology | www.ploscompbiol.org between 1 and 1:5mV, can be clearly seen both from Fig. 8-A and from the dependence of the coefficient r on the diversity degree dW B,gl , in Fig. 8-B.
The underlying mechanism leading the system from the doubleto the single-periodic response as diversity is increased can be interpreted following the prototype mechanical model of diversityinduced resonance introduced in Ref. [3]. In this model a set of interacting oscillators moving in a bistable potential is subjected to an external periodic force, which pushes the system toward the left and the right barrier alternately. If the oscillators are identical, i.e. they have the same parameter values corresponding to an underthreshold regime, then the system of oscillators cannot perform jumps on the other site of the barrier under the action of the applied periodic force. However, when the parameters are diversified (keeping constant the corresponding average value) some oscillators respond more promptly to the force and jump to the other side of the barrier, gradually pulling the rest of the system. In the present case, each neuron A i corresponds to a nonlinear oscillator of the example, while the parameter which is diversified is the activation thresholds W (i) B,gl of the glutamateinduced currents.
To show that this is the actual mechanism in action, Fig. 10 (left) illustrates the response of the heterogeneous system by depicting the time dependence of the glutamate activation variables a (i) A,gl (t) of the neurons A i , i~1,5,10,15,20, with different values of the thresholds W (i) B,gl , at the beginning of a new period in the presence of the periodic current pulse. In Fig. 10 also the raster plots for all neurons in the same time interval are shown. One can notice that the activation variables a (i) A,gl (t) behave differently from each other. Those associated to the lowest values of the activation threshold (indicated by small i values) respond stronger to the current pulse than those with the highest values of the threshold (largest values of i). The system is observed to reach the spiking regime faster than in the homogeneous case, which is shown in the right part of Fig. 10 through the comparison between the glutamate activation variable of the homogeneous system, a A,gl (t), and the average activation variable Sa A,gl (t)T~N {1 A P i a (i) A,gl (t) of the heterogeneous system. Eventually, a A,gl (t)?0 and the homogeneous system goes back to the silent state, while the average activation variables of the heterogeneous system (and their average Sa (i) A,gl (t)T) continue to oscillate around positive values, signaling the stability of the reached firing state.
Diversity in the A i ?B synapses (neuron B). In order to study the effects of added heterogeneity in the glutamate synapses located at the neuron B, one has to diversify the potential threshold parameters W (i) A,gl , see Eq. (14). For this particular case, it is possible to simplify the model into a two-neuron model with a single effective AB coupling. This is possible because heterogeneity only enters Eq. (14), while other model equations reduce to the same equations in the case of identical neurons A i , so that all neurons A i behave in the same way. Thus, the effective glutamateinduced current to the neuron B is where V A is the common value of the membrane potentials of the neurons A i . Here f (W ) is the probability density of the corresponding thresholds W~W (i) A,gl , which is assumed to be the same bell-shaped probability density f (W )~P((W {W )=dW ) as discussed above, with W~W A,gl and dW~dW A,gl .
We now consider two limiting cases of the effective current given by Eq. (18). In the limit dW %S {1 gl , when diversity is very small on the scale S {1 gl , it can be assumed that the following approximation holds, f (W )~P((W {W )=dW )&d(W {W ) and the integral (18) can be reduced to the homogeneous result, In the complementary limit of high diversity level, dW &S {1 gl , the smooth function W(S gl (V A {W )) can be approximated with Heaviside step functions H(W {V A ), and the effective current becomes gl , so to a very good approximation the effective current can be written as I eff (V )&W(l(V A {W )), where the parameter l depends on the ratio dW =S {1 gl and varies monotonously between S gl and dW {1 , as dW =S {1 gl varies between 0 and ?. The system's response at different levels of heterogeneity, dW A,gl , is presented in Fig. 11-A through the raster plots for the neuron B. A sample of time dependence is shown in Fig. 12, while Fig. 11-B shows the dependence of the coefficient r on the diversity level. In Fig. 11, it can be seen that for small values of the diversity dW A,gl the response of the neuron B presents the double periodicity of the reference configuration. Single periodicity is recovered for higher levels of diversity. At even higher values of dW A,gl , the coefficient r begins to decrease. The points of the curve corresponding to the highest values of r suggest an optimal degree of diversity dW A,gl &1 mV.
The main difference compared to the case in which noise intensity is varied, is that the response of the system remains more regular also at the highest levels of diversity considered, i.e. without random spikes appearing during the silent state and with a typical cycle well shared between a day and a night sub-period.

Discussion
In the present work we have introduced a heterogeneous multineuron version of the previously developed physiologically motivated model of the homeostatic regulation of sleep. The multi-neuron model is composed of a population of conductance-based orexinproducing neurons and a single representative glutamatergic neuron. In this model the glutamatergic and orexinergic neurons are undergoing transitions between firing and silence depending on the external circadian input and internal homeostatic mechanisms. These transitions correspond to the transitions between wake (firing) and sleep (silence), with the homeostatic mechanism being dependent on the availability of orexin. The specific aim of this study was to explore the effects of noise and diversity in the regulation of sleep-wake cycles in such a model. It is clear that diversity and noise are integral parts of all biological systems, including the orexinergic neuronal population in the lateral hypothalamus. However, the role of disorder, and especially diversity, is rarely considered in the physiologically based mathematical models of sleep-related systems [29,30,31,32,33,34,35]. To our knowledge, diversity had so far been included only in one such model, i.e. the model of interacting circadian oscillators [5], and here we present another example of the constructive role of diversity in regulation of sleep.
We have demonstrated the existence of a diversity-induced resonance, leading to a clear and strong improvement of the quality of the sleep-wake cycles, at a physiologically justified intermediate level of diversity of the orexin-producing neurons. However, only a mild improvement was found with varying noise intensity (stochastic resonance phenomena).
We have considered the simplest system with only 20 heterogeneous orexin neurons and one local glutamate neuron. Also we have used a very simple all-to-all network topology for the connections among orexinergic neurons. However, it can be expected that constructive effects of diversity will be found also in other model configurations. In the future, more realistic modifications of the model with a larger population of glutamatergic neurons and more sophisticated inter-populations connections should be considered. Furthermore, in the future studies interplay between noise and diversity should likewise be investigated, since in nature both types of disorder are normally present.
The validity of the result obtained within this model may be more general, since diversity-induced resonance is known to take place for suitable values of the parameters in general networks of interacting (non-linear) oscillators. A question then naturally arises: whether the phenomena encountered here could also characterize other systems where there is a coupling between two very different time scales or, in other words, if homeostatically regulated biological systems may take advantage from a suitable level of heterogeneity of their components.

Supporting Information
Text S1 Some details on the numerical simulation. The supporting Information file Text S1 contains some more detailed information about: A: the definition of the wakefulness time intervals during the day and the night, Dt (1) and Dt (2) , respectively, used for evaluating the quality of the sleep-wake cycle; B: the extraction procedure of the diversified threshold potentials W i and the form of the corresponding probability distribution f (W ). (PDF)