Model of the Human Sleep Wake System

A model and analysis of the human sleep/wake system is presented. The model is derived using the known neuronal groups, and their various projections, involved with sleep and wake. Inherent in the derivation is the existence of a slow time scale associated with homeostatic regulation, and a faster time scale associated with the dynamics within the sleep phase. A significant feature of the model is that it does not contain a periodic forcing term, common in other models, reflecting the fact that sleep/wake is not dependent upon a diurnal stimulus. Once derived, the model is analyzed using a linearized stability analysis. We then use experimental data from normal sleep-wake systems and orexin knockout systems to verify the physiological validity of the equations.


Introduction
Most living organisms exhibit daily biological rhythms. One of the most curious and difficult to understand is the sleep/wake cycle. Even though we are very familiar with sleep/wake, the reasons for sleep are not well understood. What is intriguing about the sleep-wake cycle is the complex interplay between the arousal and sleep states. Part of this is due to the cyclic alteration of the rapid eye movement (REM) and non-REM (NREM) stages of sleep. The distribution of REM and NREM varies dramatically across the animal kingdom, to the point that a clear pattern has not emerged. One consequence of this is that an animal model for human sleep has not been found [42]. Another complication is that sleep/wake is sometimes thought to be synonymous with the 24 hour circadian clock. This is not correct. Although sleep/wake is coupled to the circadian clock, it is regulated by a neural system distinct from the circadian timing mechanism [34]. Sleep/wake are also heavily influenced by two other factors. One is homeostatic, where the need for sleep increases with the length of the wake period. The second is allostatic, where the sleep/wake cycle is modified due to external behavioral events [32]. The circadian process is connected to the sleep-wake cycle, while a homeostatic process is inherent to the sleep-wake cycle. To be specific, it is now known that the sleep-wake cycle is independent of the diurnal stimulus, and is instead driven by homeostatic forces [12,17,1]. Therefore, the inherent periodicity of the sleep/wake cycle is not due to a periodic forcing mechanism associated with the diurnal cycle but is an inherent characteristic of the system. The objective of this study is to derive, and then analyze, a model for the sleep/wake cycle that accounts for the circadian and homeostatic influences of sleep/wake, but not necessarily the allostatic affects. This will include both REM and NREM sleep, as well as the waking state and observations about the associated brainwave activity of all states.

Biological Background
A basic assumption is that the brain has two states, sleep and wake, with sharp transitions between them. There are multiple cell groups in the brain contributing to these states, but the monoaminergic nuclei (MN) have been identified as primarily responsible for promoting the awake state [32]. Sleep is primarily promoted by the ventrolateral preoptic nucleus (VLPO), which is a group of cells associated with generating NREM and REM. A consequence of this distinct structural organization is that although the states are coupled, they have different regulatory mechanisms in the brain.

REM/NREM
The REM and NREM phases of sleep have a well defined pattern (refer to Figure 2). Regardless of the conditions of sleep onset, a person always begins in an NREM period of the cycle. After approximately 80-85 minutes the brain switches to the REM stage, and remains in this state for approximately 5-10 minutes before switching back to NREM. This 90 minute pattern continues through the night with the REM-on intervals becoming progressing longer, and the NREM-on intervals shorter. The last REM-on period is about 25 minutes. There is no probability of waking naturally during an NREM interval, as the NREM component occurs during the periods of slow wave sleep. It is postulated that this is due to NREM being the most mentally restorative of the two processes, requiring little brain function and unconsciousness. Natural wake up occurs towards the end of the last NREM period. [39] There are three different ways of triggering sleep onset: exciting sleep active neurons, inhibiting wake active neurons or a combination of the two. The same idea applies to the triggering of wake onset. Consequently, to investigate the source of sleep/wake brain function, it is necessary to examine the neurotransmitter activity of the system.

Location Specific Neurotransmitter Activity
The chemical activity of the brain involves the secretion of hormones and neurotransmitters from specific sections of the brain. The sleep-wake cycle is initiated and maintained by multiple distinct neuronal regions in the brainstem. These regions are located in the hypothalamus and are coupled to each other via electrical and neurochemical signaling. Figures 3 and 4 provide a schematic of location specific neurotransmitter signal flow in terms of excitatory and inhibitory mechanisms. Figure 1 shows the approximate locations of sleep and wake promoting areas in the human brain. Figure 5 provides a simple depiction of the switching between wake and sleep states and well as the switch between REM and NREM states, and the set of neurotransmitters associated with each, as well as their excitatory (+) or inhibitory (-) effects. The major neurotransmitter players are GABA (gamma-aminobutyric acid), orexin (hypocretin), histamine, acetylcholine, noradrenaline (norepinephrine), serotonin and dopamine. Each of these transmitters are located throughout the following neuronal groups: basal forebrain (BF), raphe nuclei (RN), laterodorsal tegmentum (LDT), pedunculopontine tegmentum (PPT), locus coeruleus (LC), lateral hypothalamus (LH), tuberomammillary nucleus (TMN), ventral periaqueductal gray region (vPAG), and ventrolateral preoptic nucleus (VLPO). Adenosine has been proposed to be the homeostatic accumulator of the need for sleep [32]. Acetylcholine effects muscle action, learning and memory. Norepinephrine effects one's control of alertness and wakefulness. Dopamine effects movement, focus and learning. Serotonin regulates mood, appetite, sleep, arousal and pain threshold. GABA contributes to sleep control and appetite control. Also note that the sections of the brain we will be focusing on are the cerebral cortex, hypothalamus and thalamus, which play key roles in homeostatic maintenance. In this way they are the most relevant to the control of the human sleep/wake system. Figure 6 shows one full period of the sleep-wake cycle model in terms of adenosine and GABA concentration levels. The point at which the curves cross indicate wake or sleep initializations, depending on whether one of the variables is increasing or decreasing. Wake Initialization occurs when GABA is decreasing and AD is increasing. Sleep initialization occurs when AD is decreasing and GABA is increasing.  Adenosine, which is a neuromodulator, plays a major role in the sleepwake process. It is a byproduct of adenosine-triphosphate (ATP), the main "energy currency" of the brain, whenever neurotransmitters are fired, adenosine is produced. Adenosine is also continuously converted back to ATP via the metabolic processes that use it as an energy source [37]. The adenosine concentration that we are most interested in is that of the basal forebrain. While inhibiting the expression of all wake active neurotransmitters and disinhibiting the expression of those found in the VLPO [7], adenosine levels in the basal forebrain increase during prolonged wake and fall during recovery sleep, providing an intrinsic periodic force associated with waking [30,29]. The linear relationship between between adenosine accumulation and sleep propensity provides support for using adenosine as representative of sleep propensity. GABA in the VLPO can be viewed as the other major player in the sleep-wake system. GABA A receptors in the brain mediate the tonic (or slow) conductance [40]. This slower type of signaling plays a key role in cell excitability and sleep homeostasis [40]. Neuron populations in the VLPO are inactive during the vigilance state, and have a high firing rate and level during NREM, especially during delta wave (slow wave) sleep [40]. Being the major inhibitory neurotransmitter in the central nervous system, GABA is a key player in many physical systems.
Referring to the Table 1, we can see the detailed activity associated with each neurotransmitter and brain area, allowing for the development of a system of variables. As a simple example, OX LH is the variable associated with the chemical concentration of orexin in the lateral hypothalamus. Since we only discuss orexin in one location, the subscript LH can be dropped. Orexin expression is promoted by acetylcholine in the LDT/PPT and noradrenaline in the LC. Orexin expression is inhibited by adenosine, serotonin, GABA and noradrenaline. Orexin promotes the release of noradrenaline, serotonin, histamine, GABA in the basal forebrain, acetylcholine, and dopamine, as well as promoting its own release (autoexcitation). These dynamics indicate the need for a linear combination of acetylcholine, noradrenaline, serotonin,  dopamine, histamine, GABA, adenosine, and orexin to express the rate of change of OX concentration.

Materials and Methods
In order to test the functionality of the model, we use numerical methods to determine how the system responds to a variety of perturbations. In the calculations to follow, MATLAB 2009a was used to run all the programs written for the Human Sleep/Wake system described here. However, a fourth order Runge-Kutta method was used for solving the system of equations, instead of using one of the MATLAB ode solvers because of the unique perturbations used in the sleep experiments described.

Previous Work
A good survey of the early sleep/wake models can be found in [35]. We will only comment here on those most relevant to our development.
In a recent paper by Phillips and Robinson [28], a model was proposed in order to improve quantitative understanding of the human sleep/wake cycle, which is similar to the intentions of this study. They use the following three neuronal populations to model the human ascending arousal system: monoaminergic (MA) nuclei, acetylcholine related (ACh) nuclei, and populations in the VLPO. The VLPO is then split into circadian and homeostatic drives, which are represented by periodic forcing functions. Each population has its dynamics modeled by a sum of mean cell body potentials and mean neuronal firing rates subjected to the periodic forcing of the VLPO. While their model consistently predicts sleep/wake cycle dynamics observed in previous studies, it is important to note that the homeostatic and circadian drives are represented as periodic forcing rather than dynamics inherent to the physical system.
[10] developed a model of the network of neurons controlling behavioral state transitions in mice, while also identifying special features of human and mouse sleep. It is similar to the Phillips and Robinson model (as well as the model being developed here) in that it defines three main processes being modeled: Sleep activity, Wake activity, and REM activity. Each population is defined by a set of three equations, determined by a wide range in timescales. It incorporates homeostatic REM and NREM sleep drive and uses known sleep active neurotransmitters to determine parameters. It then utilizes extensive experimental data gathered from mouse sleep/wake states to verify the accuracy of the model. A limitation of their approach concerns the separation of vigilance states, which are characterized in a non-quantitative fashion. Similarly, "activity states" are used rather than specific neurotransmitter concentration and/or neuron firing data.
Taking a somewhat different viewpoint, [31] use the "flip-flop" dynamics proposed in [32] to construct a biologically-based mathematical model that accounts for several features of the human sleep/wake cycle. This differs from our model in several ways, but two examples are that they assume, incorrectly, that the SCN projects directly to the VLPO, and an external periodic forcing mechanism is used to drive the model.
In another recent paper by [9], a model of the rat sleep/wake regulatory network is formed utilizing similar techniques as those used to model a mouse sleep/wake regulatory network in her aforementioned work [10]. In this model, neurotransmitter concentrations are coupled to neuronal group firing rates to create a more complete view of the hypothalamic sleep-wake regulatory network. The model is composed of three parts: wake/sleep active neuron population firing rates, wake/sleep active neurotransmitter concentrations and a homeostatic regulator. All physical parameters were determined using extensive experimental data. This paper emphasizes the importance of utilizing neuron firing rates as well as neurochemical concentrations to form a physiologically accurate model of sleep/wake regulation. A downside of this model lies in the use of steady state saturation functions and Heaviside functions. Rather than using inherent neurotransmitter dynamics, the functions of neurotransmitter release are approximated using sigmoidal functions. The homeostatic drive is a physically unquantifiable force composed of Heaviside switching mechanisms.

Model Formulation
In this model for the sleep/wake cycle, we are considering three principal states. They are the wake state (W ), the NREM sleep state (S) and the REM sleep state (R). Table 1 provides a framework for the detailed interactions of all the neurotransmitters involved in the sleep-wake system. Given the number of variables involved, before deriving the equations it is worth describing the dynamics in more qualitative terms.
The sleep-specific and wake specific neurotransmitters operate on a slow time scale, measured in terms of hours, while the REM and NREM specific neurotransmitters operate on a fast time scale system, measured in terms of minutes. The two neurochemicals associated with the slow time scale are adenosine in the BF and GABA in the VLPO. We use these two concentrations as variables indicating promotion of wake and sleep. High adenosine concentrations in the BF promote wakefulness, while high GABA concentrations and fast firing rates in the VLPO promote sleep. Adenosine's slow build up and decay is associated with wakefulness, while GABA in the VLPO is associated with sleep. These observations are the basis for the hypothesis that adenosine is the principal component in determining wake, and GABA V LP O is the principal component in promoting sleep [40,30,29]. We also know from these references that AD disinhibits the VLPO during slow wave sleep, providing a direct connection between the building adenosine sleep pressure during waking and the gradual relief of this pressure during sleep. The other neurotransmitters are assumed to be associated with the fast time scale dynamics. Adenosine and GABA V LP O fluctuate with a period of 24 hours, while the other neurotransmitters fluctuate with a period of several minutes.
An important inherent quality of the sleep-wake cycle is that it should maintain it's state transitions in the absence of external periodic forcing. Mathematically, this indicates the need for the system to have a stable (attracting) limit cycle. This solution should also be stable to certain perturbations, such as those that might occur from genetic mutations causing the production of too much or too little of a neurotransmitter, or lacking the proper neuron groups to produce a neurotransmitter at all. Perturbations may also occur via the removal of specific neuron groups, blockage of neurotransmitter bonding sites, and external forcing. It is well known from experimental data that perturbations to neurotransmitters don't alter the two-state nature of the sleep-wake system [30,29,3,27,16,22,21,20,5,11,23,8,14,4]. As an example, cholinergic agonists applied to the LC in cats leads to increased REM sleep, but the sleep-wake cycle as a whole remains periodic [8]. Therefore one of the central goals will be to locate the area in the phase plane of GABA V LP O and AD where a limit cycle exists and is stable.

Kinetics
For each neurotransmitter involved in the sleep-wake system, we identify a location specific variable representing the effective neurotransmitter concentration. We assume first order kinetics for the neurotransmitters operating on the fast time scale. The fast time system may be formed from the observations compiled about the inhibitory and excitatory nature of the neurotransmitters in Table 1 combined with the dynamics arising in a synapse. If X i is representative of a neurotransmitter in region i, we have Note that inhibitory effects are considered pre-synaptic and are therefore controlling flow from the neuron into the synaptic gap. Despite the fact that inhibitory and excitatory effects can occur pre-or post-synaptically, for simplicity we will assume the blocking or binding process of inhibition as only controlling flow into the gap. The Release function contains rate parameters and neurotransmitter variables that excite and inhibit the expression of X i . The Uptake function depends on X i and its associated vesicular uptake transporter, if existent. The Degradation function depends on X i and an associated rate of decay. The uptake transporters are modeled using known logistic dynamics from neurotransmitter uptake transporter microdialysis measurements [26].
As an example, the release of histamine is promoted by orexin, noradrenaline and acetylcholine in the BF and inhibited by GABA in the VLPO. Histamine is autoinhibitory, like most of the monoamines. The equation for histamine's rate of change, like all of the other neurochemicals involved, also requires an adenosine term. Histamine in the human sleep-wake system has no known uptake mechanism, and so only degrades as a method of removal.

Therefore the equation for histamine is
where the c i 's are rate coefficients. The rest of the equations are formulated in the same fashion.

Uptake Mechanisms
The neurotransmitters involved in the sleep-wake system, as with any neuronal network, have natural rates of removal from the synaptic gap once released. These rates can either be contributed to by degradation, neurotransmitter uptake or a combination of both. All of the neurotransmitters involved degrade over time, but only some have neurotransmitter uptake mechanisms. We will not be considering any of the smaller scale chemical dynamics of the neurotransmitter uptake mechanisms. We care only whether a neurotransmitter is removed by an uptake mechanism, and will not differentiate between plasma membrane-bound transportation or receptor uptake. The key take away from the uptake transportation data is that only GABA, noradrenaline, serotonin and dopamine have significant synaptic uptake transport mechanisms in humans. These mechanisms are characterized by the data given constants GA1, hN ET , hSERT and hDAT for GABA, noradrenaline, serotonin and dopamine, respectively [26]. These neurotransmitter transporters appear in the uptake term of the equations, always multiplied by the neurotransmitter it is dependent on to indicate a linear removal rate. As an example, hSERT appears in equation 9 multiplied by −c 52 S, so as to act as a removal mechanism for serotonin when it is present. However, hSERT does not multiply the natural degradation term −c 50 S.

Nonlinear Dynamics
There are two neurotransmitters, AD and GABA V LP O , operating on a slow time scale that define the sleep-wake cycle. The qualitative portrait in phase space is known from basic sleep-wake cycle activity. We can construct the desired closed orbit and concentric curves from the nullclines of the AD and GABA V LP O equations. We identify the bounding region in the phase space activity by the box contained between (AD, GABA V LP O ) = (AD max , GABA min ) and (AD min , GABA max ). We choose small values of AD min = GABA min = 0.001 rather than exactly zero to prevent inclusion of a potentially unstable steady state. We can scale the size of this box according to the values of AD max and GABA min . These two points are represented by the intersections of the AD and GABA V LP O nullclines. To describe the activity between these two points, we begin with the fact that when AD is at it's maximum value, GABA V LP O is at it's minimum value (0.001). This signifies that at the peak of wake activity, sleep activity is at its minimum. As we move along the AD nullcline, the concentration of AD decreases to 0. We begin by forming the slow time system that will later be shown to be responsible for the limit cycle dynamics of the system. This system is comprised of two equations: the adenosine equation that fluctuates with wake activity, and the GABA in the VLPO equation that fluctuates with NREM sleep activity. We consider NREM sleep to be the ground state of the sleep-wake system, and all other oscillations to be deviations thereof. The adenosine equation must involve a constant source term, due to adenosine's constant creation as a by-product of neuronal activity, as well as a decay term that is dependent on adenosine's natural rate of degradation as well as it's loss from the synaptic gap through the its disinhibition of GABA in the VLPO [7]. It is well known that mechanism for receptor binding is nonlinear [26], and so the disinhibition effect of GABA V LP O via adenosine should be nonlinear. It also follows that this nonlinear term should appear in both the equations for AD and GABA  tion term and a contribution from adenosine as well as the inhibitory and excitatory of other neurotransmitters. A constant loss term is required in the GABA V LP O equation as a representation of its quick breakdown time in the synaptic gap [40]. Though it is being broken down quickly and continuously, GABA V LP O still provides a tonic effect on the sleep-wake system as a whole [40]. Based on these requirements as well as the observations of the glycolysis reaction [33], we obtain the following slow time system: To provide a more physical explanation of the assumed reaction system, we know that adenosine is a bi-product of the constant source of energy in the nervous system, adenosine triphosphate (ATP). We know that adenosine is continuously converted back into ATP [37]. From experiments in the VLPO, we know that activation of VLPO neurons through adenosine-mediated inhibition contributes to the somnegenic actions of adenosine, and only adenosine acts directly on the VLPO in the sleep-wake system [7,3,30,29]. During wake, increased levels of adenosine promote sleep by disinhibiting the release of GABA in VLPO. Adenosine levels in this area increase during wake and decrease during sleep [7]. In this way the physical foundation for the glycolysis system is similar to and provides a good framework for our slow time scale system.

Characterizing REM/NREM
REM and NREM sleep are determined by a nonlinear combination of neurotransmitters. Specifically, GABA V LP O promotes the state of REM sleep, as does acetylcholine from both the LDT/PPT and the BF. While NREM can be seen as an absence of the vigilance state and the contributions from GABA, REM sleep must incorporate contributions from neurotransmitters active in the wake state as well as those in the NREM state. Letting R(t) represent the REM sleep state, we assume that R + ΦR + ΓR = 0. To determine Φ, note that the transition from the waking state to the sleeping state, and vice versa, occurs when GABA V LP O = AD. Given that the REM state has an asymptotically stable steady state during wake, which occurs when AD > GABA V LP O , and has a limit cycle during sleep, which occurs when AD < GABA V LP O , then Φ contains a term of the form (GABA V LP O − AD)(R 2 − R 2 0 ). When entering wake, so GABA V LP O = AD andȦ D > 0, there is an impetus for the impending dormancy of REM and the resulting mathematical requirement is that Φ(Ȧ D, AD) > 0. Likewise, when entering sleep, so GABA V LP O = AD andȦ D < 0, the requirement for the oncoming REM limit cycle is Φ(Ȧ D, AD) < 0. This gives us

Some Special Cases
Given the complexity of (3.13), it is worth considering the solution in particular situations and then using this to help estimate the oscillatory response of equation (3.13). For example, the period for R can be determined from the fact that there are three to four normal REM oscillations over a normal night's sleep. So, when GABA V LP O = AD, Given the slowly varying nature of AD, then R = e rt , where r = 1 2 (−βu ± β 2 u 2 − 4Γ) and u = (AD) . Since u strongly controls the decay, and to have Γ control the number of REM oscillations, we require r = −u ± i √ γ, where √ γ is approximately the number of oscillations. Thus, β = 2 and Γ = γ + u 2 . From this we have that (3.15) where α is determined numerically and 8 < γ < 16.
Based on Lienard's Theorem, we know that the above oscillator has a unique, stable limit cycle [36]. Thus the REM oscillations will continue to function in a stable, periodic fashion even when subjected to reasonable perturbations. Evidence of this result is in the Orexin Knockout section.

Analysis
To find the values of the coefficients of the system and to ensure stability, we present several methods of analyzing the equations. We begin with the limit cycle coming from the Slow Time System. As stated earlier, the adenosine equation drives the wake state dynamics, and the GABA equation drives the sleep state dynamics. The specific equations are It is not possible to determine the rate constants with the available experimental evidence. So, we take a different tack and establish bounds on the constants that will produce the responses known to occur in sleep and wake. To begin, we compute the nullclines of the slow time system in order to create a trapping region for the limit cycle and to determine the coefficients of the system by using an approximate fixed point. We have four coefficients to determine (k 1 , k 2 , k 3 , k 4 ) as well as reasonable values for the parameter . Plugging the approximate fixed point into both nullclines will give us two equations with five unknowns. In order to make the system well determined, we observe the maxima and minima of the equations as well as the linearized slow time system with the approximate fixed point. Knowing the form of the determinant and the trace of this system needed to maintain stability, we can determine the rest of the unknowns. We also know that each of these equations is zero separately when AD = AD max and GABA V LP O = GABA min , as well as when AD = AD min and GABA V LP O = GABA max . We let AD max = GABA max = 2 and AD min = GABA min = 0.01. Taking = 0.3, and using the stated values for the max and min concentrations it follows that k 1 = 0.49, k 2 = 0.1, k 3 = 0.3. k 4 = 0.15. We also note that 0.29 < ≤ 0.32 for the system to remain stable.

Remaining Coefficient Estimation
Now that the coefficients for the slow time system are determined, we can use the requirements that the eigenvalues have either pure imaginary form or negative real parts in the complex form to determine the coefficients of the entire system. We have the solution vector of the linearized system as x = [GABA BF w , GABA BF s , OX, H, ACh BF , ACh LDT /P P T , N A, S, DA, AD, GABA V LP O ] T . The matrix of coefficients and linearized values for the system are given in Table 2.
In order for the system to be stable, we require that the real parts of the eigenvalues to be negative. Restricting the coefficients to all be between 0 and 1, we use an algorithm that generates coefficients that keep the real parts of the eigenvalues negative. We also note that the coefficients in GABA V LP O equation need to be small compared to the rest of the coefficients in order to keep bounded.

Sub-Systems
When analyzing such a large system of equations with undetermined coefficients, it is helpful to break the large system into smaller sub-systems and observe the dynamics on a local scale. We can utilize the neurophysiology of the physical neurotransmitter sub-systems to determine the coefficients as well as the stability profiles of the individual sub-systems. We break the system up into to following categories: NREM Specific, Sleep Specific (NREM Inclusive), Wake Specific, REM + NREM Specific, REM + Wake Specific, Wake Inclusive, Sleep Inclusive, REM Inclusive. Note that there is no set of equations associated with REM only. The categories represent what neuron populations are active during the corresponding states. We observe the steady-state behavior of each sub-system. Since we know the relative magnitude of the neurotransmitter amounts as well as their dynamic behavior during the states of REM, NREM and Wake, the coefficients of these smaller, more manageable systems should be more easily attainable. This also allows us to observe the more manageable systems for potentially extreme behavior, and provides us with bounds to prevent this from happening. This analysis does not give us specific values for the system coefficients, knowing the relative values makes visualizing the system easier as well as simplifying the eigenvalue analysis.

NREM Specific
The NREM specific system involves the neuron population only active during the state of NREM sleep. The neuron population that is only active in NREM sleep is the GABAergic population in the BF. We assume all equations but dGABA BF s dt to be equal to zero, making each neurotransmitter other than GABA BF s constant, and assume the other neurotransmitters in the equation to be at their steady state values when only GABA BF s is active. This is not necessarily a physiologically accurate view of the system for an extended length of time, but we can view a snapshot of this activity. We follow this same criteria for the rest of the small systems.

Wake Specific
The Wake specific system involves neuron populations with the most activity during the wake state. The monoaminergic neuron populations (including OX, H, NA, S and DA) are most active in this state.

Sleep Specific (NREM Inclusive)
The Sleep Specific system involves neurotransmitters that promote sleep states (REMS and NREMS). The populations that are active in this state are GABAergic populations in the VLPO and BF.

REM + Wake Specific
The REM and Wake Specific system involves neurotransmitters that promote REMS and the Wake state. The neuron populations that are active during these states are GABAergic populations in the BF and cholinergic populations in the BF and LDT/PPT.

REM + NREM Specific
The REM and NREM Specific system involves neurotransmitters promoting REM sleep and the NREM sleep. The neuron population that is active during these states is the GABAergic population in the VLPO.

Wake Inclusive
The Wake Inclusive system involves neurotransmitters that promote the Wake state as well as those that promote both Wake and REMS. The neuron populations that are active during these states are the GABAergic population in the BF, the cholinergic population in the BF and LDT/PPT, and the monoaminergic populations.

Sleep Inclusive
The Sleep inclusive system involves neurotransmitters associated with NREM specific, REM and NREM specific, and Wake and REM specific states. The neuron populations that are active during these states are GABAergic populations in the VLPO and BF, and cholinergic populations in the BF and LDT/PPT.

REM Inclusive
The REM Inclusive system involves neurotransmitters associated with the Wake and REM specific and REM and NREM specific states. The neuron populations that are active during these states are the GABAergic population in the VLPO, the wake promoting GABAergic population in the BF as well as the cholinergic populations in the BF and LDT/PPT. This analysis allows for an approximation of the coefficients, but it does not yield exact values for these coefficients. Figures 11 and 12 show the change in concentration levels of Acetylcholine in the BF and GABA in the BF, respectively over a period of 24 hours. The remainder of the neurotransmitters involved show the same fluctuations as that of Acetylcholine, whereas GABA in the sleep active portion of the BF is unique.

Results and Discussion
The results of the implemented normal, unperturbed, system are shown in Figures 11-14. Figure 11  which is a measure of wake pressure. Both exhibit a rapid transition between sleep and wake states as well as ongoing periodic behavior for nine days. The system runs like this for much larger periods, demostrating that unless perturbed, the sleep-wake cycle is a self-sustained, periodic oscillation. The period of each cycle is approximately 24 hours -the wake state occupies roughly two-thirds of this time, while the sleep state makes up the remaining one-third. This is consistent with our knowledge of the structure of the human sleep-wake cycle. Figure 12 depicts one cycle in the GABA − AD phase plane, showing the exact behavior we required from our previous nullcline analysis. Figure 13 is the isolated REM curve for 30 hours. Figure  14 shows the AD, GABA and REM curves simultaneously in a sample 48 hour time slice. The REM/NREM oscillations respond quickly to the transitions between wake and sleep, and cycle through three to four oscillations in the designated eight hour sleep period. From the biological background and neurotransmitter details, this displays physiologically reasonable results.

Orexin Knockout Experiment
Narcolepsy is a neurodegenerative disease where there is a loss of neurons containing orexin (OX) [13]. It has been shown to exist across many mammals, and has been studied extensively in mice. Here, we take information from the sleep-wake cycles of orexin knock-out (which means the targeted destruction of orexin neurons) mice. These mice exhibit the narcolepsy symptoms of sleep fragmentation, sleep onset during REM sleep states and increased spontaneous REM and NREM sleep during the vigilance state [27,41,4,25,15].
A reasonable method to represent the orexin knock-out perturbation to the system is to change the flux of the OX amount equation. Rather than take orexin out all together, we limit the rate at which orexin is produced, but keep the removal factors the same. This is a reasonable way to represent the physiological destruction of orexin neurons, since none of the uptake and removal mechanisms are necessarily effected. It has been shown that levels of orexin decrease in human cerebral spinal fluid by at least 80 percent [13], and so we multiply orexin producing elements of the OX equation (ACh LDT /P P T , N A, OX and AD) by 0.2. Implementing these changes numerically, the periodicity and limit cycle of the system is maintained with orexin decrease. This is illustrated in Figure 15, which clearly shows the same characteristics as seen in Figure 11. However, we also don't see the stereotypical changes to the transition rates between behavioral states characteristic to narcolepsy. We have shown that the system can withstand perturbations [27], but may need to think of implementing another method of making genetically specific perturbations. This also shows that orexin may have to play a larger role in measuring behavioral state instability, and may have to withstand more non-linear changes to its equation when investigated further. [6] describes experiments related to perturbations in an individuals sleep/wake cycle. The trials performed were carried out in a controlled situation with a group of 40 adolescents over a period of one week. The perturbations performed were delayed sleep and wake onset, prolonged sleep and wake periods and forced wake. These perturbations were achieved without the use of external chemical stimuli, and so were essentially achieved in a noninvasive fashion. Salivary cortisol and melatonin levels were recorded daily as well as frequent EEG readings for accurate brainwave activity. The initial time and length of the perturbation period were recorded, and so this data was easy to translate into perturbations in the model. About 50 perturbations were introduced (which translates to a period of about 10 days) and their effect on the system was recorded using graphical and data analysis methods.

Sleep Camp Data
The perturbations to the system were only introduced to one variable -adenosine. Forced sleep and wake were achieved by forcing AD to its maximum at the appropriate forced waking initialization, and forcing AD to its minimum at the appropriate forced sleep time. These sleep and wake times as well as the prolonged sleep and wake periods were taken from a data set of seven adolescent individuals. Once perturbed the amount of times done so in the physical experiment, we allow the system to run without perturbation for several cycles. This entire set of perturbations and the response is depicted in Figure 16. Figure 17 shows what the subject's cycle would look like without perturbation. One can see that the perturbations cause small scale (on the order of a day or two) changes in the length of sleep or wake state. However, once allowed to run unperturbed, the system settles back into a regular 24 hour period. Initial wake and sleep times have drifted, but that is to be expected. When run for a longer time, the system will eventually drift back into its original phase with the same wake and sleep initialization times as the subject started with. We know this is an accurate response to these non-drug related, small perturbations(REF) . Figures 18 and 19 show a quantitative measure of the changes in period in the unperturbed vs. the perturbed cycles. Figure 20 demonstrates the upholding of limit cycle dynamics in the face of perturbations, and also that the AD − GABA phase plane retains its size and shape but drifts in phase.

Conclusions
We have presented a model of the human sleep-wake generating system that identifies specific neuronal regions in the brain with their associated neu- rotransmitters and projections. It accounts for the homeostatic influences of sleep-wake, and has the potential to encompass circadian influences as well. It improves upon previous models of sleep-wake cycles in that the oscillations are self-sustaining without the use of external periodic forcing, and it is formulated specially to look at the specific sleep-wake cycle of humans. We use the known neurochemical dynamics of the human sleep-wake cycle to derive equations and form the stability criteria for mathematical phenomena. Since based on limit cycle dynamics, the system runs with the contributions of all the necessary neurotransmitters and with no need of an external driving force. This system gives reasonable predictions for "natural" sleep bout duration, wake bout duration, REM bout duration and NREM bout duration, as well as estimated natural sleep initiation and wake initiation times. The NREM sleep state is considered the ground state from which all other activity deviates, and REM activity is determined from the difference of the wake and NREM states, as well as the contributions from REM active neurotransmitters. Using experimental data from a mammalian orexin knockout system as well as simply and non-chemically perturbed human sleep-wake cycles, we were able to visualize the perturbed cycle and compare to the experimental results without instability occurring. We can draw a number of conclusions from the resulting model. In the human systems, given forced prolonged waking, sleep duration reacts accordingly to build up appropriate wake pressure. Once the cycle is restored to normal sleep onset and wake onset, the same period for an unperturbed system is maintained. The system exhibits some phase drift with this forcing. Given forced prolonged sleep, wake duration also reacts accordingly to build up appropriate sleep pressure. Once the cycle is restored to normal sleep and wake onset, the periodicity is maintained. The system exhibits some phase drift with this forcing. Given forced wake up during sleep, the subject does not stay awake for long, and falls back into a normal sleep schedule, yielding only some phase drift that is corrected after a few cycles. A subject's REM/NREM cycle goes through three or four oscillations per state (REM or NREM). The first NREM sleep stage is triggered by sleep onset, or when adenosine levels fall below a certain threshold and GABA becomes the dominant neurotransmitter. The last REM cycle is turned off by wake onset, or when adenosine pressure builds back up and GABA decays. Subjects naturally fall asleep and wake up at the intersection points of the GABA and Adenosine thresholds -falling asleep when enough sleep pressure from AD has been accumulated, and waking up when the sleep pressure has decayed once again and when the subject has gone though three to four cycles each of REM and NREM sleep. We can see that given no external forcing (no stimuli, drugs, etc), a "normal" human's sleep-wake schedule unfolds in a periodic, 24 hour fashion that endures in-  definitely unless given some sort of perturbation.
It is important to note the implications of the differences between homeostatic and circadian drive, and the importance of both of these drives' effect on sleep. Unlike most recent neurochemically driven models, ours has no need for external forcing to keep it periodically driven. The periodic, limit cycle activity arises naturally from the neurotransmitter dynamics and the equations themselves. This demonstrates an intrinsic, homeostatic force that keeps the human sleep-wake cycle functioning without a change in periodicity, regardless of whether circadian control is included or not. The perturbations in the experimental section demonstrate that even when subjected to strong forcing, the sleep-wake cycle retains its natural 24 hour period, and shifts in phase are the only real artifacts of outside forcing. The next step would be to develop a more complete predictive model of the human sleep wake cycle for use in the scientific and medical fields.

Acknowledgements
This work was supported, in part, though NSF Grant DMS-0636358. Also, we would like to thank the following for their aid in this project: Dr. Mary Carskadon (Brown University), Dr. Tom Scammell and the Scammell Lab (BIDMC, Harvard Medical School).