Mathematical Biology - Semantic Scholar

J. Math. Biol. (2010) 60:615–644 DOI 10.1007/s00285-009-0276-5

Mathematical Biology

A mathematical model of the sleep/wake cycle Michael J. Rempe · Janet Best · David Terman

Received: 1 August 2008 / Revised: 1 May 2009 / Published online: 26 June 2009 © Springer-Verlag 2009

Abstract We present a biologically-based mathematical model that accounts for several features of the human sleep/wake cycle. These features include the timing of sleep and wakefulness under normal and sleep-deprived conditions, ultradian rhythms, more frequent switching between sleep and wakefulness due to the loss of orexin and the circadian dependence of several sleep measures. The model demonstrates how these features depend on interactions between a circadian pacemaker and a sleep homeostat and provides a biological basis for the two-process model for sleep regulation. The model is based on previous “flip–flop” conceptual models for sleep/wake and REM/NREM and we explore whether the neuronal components in these flip–flop models, with the inclusion of a sleep-homeostatic process and the circadian pacemaker, are sufficient to account for the features of the sleep/wake cycle listed above. The model is minimal in the sense that, besides the sleep homeostat and constant cortical drives, the model includes only those nuclei described in the flip–flop models. Each of the cell groups is modeled by at most two differential equations for the evolution of the total population activity, and the synaptic connections are consistent with those described in the flip–flop models. A detailed analysis of the model leads to an understanding

This work was partially funded by the NSF under agreement 0112050 and by the AFOSR grant FA9550-06-1-0033 to D. Terman and J. Best. M. J. Rempe (B) Mathematical Biosciences Institute, Ohio State University, Columbus, OH 43210, USA e-mail: [email protected] J. Best · D. Terman Department of Mathematics, Ohio State University, Columbus, OH 43210, USA e-mail: [email protected] D. Terman e-mail: [email protected]

123

616

M. J. Rempe et al.

of the mathematical mechanisms, as well as insights into the biological mechanisms, underlying sleep/wake dynamics. Keywords

Sleep · REM · NREM · Sleep homeostat · Circadian pacemaker

Mathematics Subject Classification (2000)

92B99 · 37N25

1 Introduction Despite nearly a century of study, sleep and its underlying processes hold many mysteries. It remains unclear how identified brain regions interact to bring about the different stages of sleep and wakefulness, how the timing of sleep depends on the length of time spent awake and work load, and how pathologies associated with sleep, such as narcolepsy, arise. Mathematical modeling has been used extensively to address these issues (Borbély 1982; Daan et al. 1984; Jewett and Kronauer 1999; McCarley and Hobson 1975; Strogatz 1987); however, the success of these models in predicting sleep propensity under normal conditions, forced desynchrony and sleep deprivation comes despite a lack of a detailed understanding of biological mechanisms that underlie the observed sleep and EEG dynamics. In this paper, we develop a mathematical model for the human sleep/wake cycle that is based upon physiological processes and interactions of neuronal cell groups. Thus, the model can be applied to test hypotheses concerning the neural circuitry and how identified brain regions interact to bring about the different stages of sleep and wakefulness. One goal of this paper is to further strengthen the biological underpinnings of the two-process model for sleep regulation (Achermann 2004; Achermann and Borbély 1990; Borbély 1982; Borbély and Achermann 1999). This model addresses the timing of sleep and wakefulness, positing that sleep regulation arises via the interactions of a roughly 24-h circadian process (C) with a wake-time-dependent homeostatic process (S). The model conceptualizes a reference level of sleep that is homeostatically maintained between two thresholds that are in turn modulated by the oscillations of C. Process S builds up during waking hours, so that large values of S correspond to a strong sleep drive. When S reaches an upper threshold, the system “falls asleep” and S begins to decay. S continues to decay during sleep until it reaches a lower threshold and the system “wakes up.” The dynamics of S are chosen based on the observation that EEG activity in the delta range of 0.5–4.5 Hz —slow wave activity (SWA)—reflects sleep intensity (Daan et al. 1984). Spectral analysis of EEG activity not only showed that SWA is highest in the early part of the sleep episode and, on average, decreases over the course of the night as sleep need is satiated, but also provided a quantitative measure of its time course. The increase is fit with a saturating exponential function, and, in the simplest version of the model, the decrease of S during sleep is also fit to an exponential function. The two-process model has been quite successful in predicting the timing of wake and sleep under normal conditions, acute sleep deprivation and naps despite a lack of a detailed understanding of neural mechanisms that underlie the observed sleep and EEG dynamics. We note that while mechanisms underlying the circadian pacemaker are well known, the nature of the homeostatic process remains poorly understood. It is also mysterious why the timing of waking up and falling asleep occur when

123

Mathematical model of sleep

617

some combination of the homeostat and circadian drive reach two somewhat arbitrarily constructed thresholds, especially as it is now clear that these processes are not independent (see for example Colwell and Michel 2003). The past decade has brought significant progress in identifying the neural substrates of sleep regulation. Our model builds on the work of Saper et al. (2001) who introduced the idea of a hypothalamic sleep/wake switch. Their conceptual model relies on evidence that sleep-promoting neurons within the ventral lateral preoptic nucleus (VLPO) have mutually inhibitory connections with wake-promoting monoaminergic cells in nuclei including the tuberomammillary nucleus (TMN), the dorsal raphe nucleus, and the locus coeruleus (LC). They proposed that these mutually-inhibitory connections function as a flip–flop switch: any increase in activity on one side generates positive feedback via disinhibition, avoiding dangerously slow transitions between sleep and wake. Saper and colleagues also proposed that neurons containing the neuropeptide orexin (also called hypocretin) stabilize the switch by reinforcing the monoaminergic cell groups (Saper et al. 2005). Orexinergic cells of the lateral hypothalamus project to many of the monoaminergic nuclei, providing circadian-modulated excitation. Saper et al. hypothesized that, beyond simply promoting wakefulness, orexin functions to stabilize the switch against over-frequent transitions. In support of this hypothesis, experiments show that the loss of orexin causes many symptoms found in narcolepsy, including frequent transitions into sleep and fragmented sleep episodes, indicating decreased stability of both wake and sleep states (Chemelli et al. 1999; Mochizuki et al. 2004; Sakurai 2007). Saper and colleagues have also recently extended the flip–flop model to incorporate REM and NREM sleep stages (Lu et al. 2006), positing an additional flip–flop switch between brain stem areas involved in controlling REM and NREM sleep. The ventrolateral part of the pariaqueductal grey matter (vlPAG) and the lateral pontine tegmentum (LPT) were identified as REM-off regions and the sublaterodorsal nucleus (SLD) as well as the precoeruleus (PC) were identified as REM-on regions. Experiments suggest that the VLPO contains two distinct populations of cells: the cluster portion (clVLPO) and the extended portion (eVLPO). The NREM-promoting nuclei receive inhibition from eVLPO, while the REM-on populations are inhibited by monoaminergic nuclei; these connections link the two flip–flop switches into a verbal model for the neural circuit regulating sleep including REM and NREM activity. Recently, several physiological models have been developed that are based in part on the flip–flop concept (Diniz Behn et al. 2007, 2008; Phillips and Robinson 2008, 2007; Tamakawa et al. 2006), but several of these studies investigate sleep in rodents, not humans (Diniz Behn et al. 2007, 2008; Tamakawa et al. 2006) while others model human sleep, but do not include REM dynamics (Chou 2003; Phillips and Robinson 2008, 2007). One model combines aspects of both the two-process model and the sleep flip–flop conceptual model in human sleep (Phillips and Robinson 2008), but is focused primarily on sleep deprivation, and does not include a REM flip–flop or an analysis of the effects of orexin. Our model is based on the flip–flop switches, but also reproduces the important features of the two-process model. The model is consistent with results of experimental studies used to investigate how observed features of sleep arise from the circadian

123

618

M. J. Rempe et al.

and homeostatic processes and from their interactions. Our results demonstrate that the neuronal components in the flip–flop model, together with a sleep homeostat and a circadian pacemaker, are sufficient to account for several features of the sleep/wake cycle, including a loss of orexin leading to more frequent transitions as in narcolepsy. However, some features of REM/NREM dynamics require additional hypotheses. We provide a detailed analysis of the model that leads to an understanding of the mathematical mechanisms, as well as insights into the possible biological mechanisms, underlying the sleep/wake dynamics. 2 Materials and methods We describe our model in two steps. First we describe the model for the sleep/wake flip–flop switch and then we add cell groups to this model for the REM-NREM switch. 2.1 The sleep/wake flip–flop model A schemata of our model for the sleep/wake and REM-NREM switches is shown in Fig. 1. The model includes the sleep-promoting neurons in the ventrolateral preoptic region of the hypothalamus (VLPO), the wake-promoting monoaminergic cell groups (AMIN), orexin neurons (ORX), a circadian pacemaker corresponding to activity within the suprachiasmatic nucleus (SCN), and input from cortical areas. We also assume that there is a sleep homeostat, (HOM), that increases while awake and decreases during sleep. Note that the sleep-promoting and wake-promoting cells inhibit each other. We assume that output of SCN inhibits VLPO. While experiments have demonstrated that direct connections from SCN onto VLPO are sparse, SCN does project onto other cells groups including DMH (Chou et al. 2002) and this has an inhibitory effect on VLPO (Chou et al. 2003). In our model, AMIN receives excitatory input from ORX. The orexin neurons receive excitatory input from SCN and inhibitory input from VLPO. Fig. 1 A schemata of the model

SCN

IoV

ORX

VLPO

AMIN HOM

eVLPO

REM-off

123

REM-on

inhibition excitation

IAo

Mathematical model of sleep

619

Experiments have demonstrated that these cells are typically inactive while asleep (Sakurai 2007). The connections used in the sleep–wake portion of the model are the same as those used in Phillips and Robinson (2007), and are similar to those used by Tamakawa et al. (2006). The VLPO and AMIN are each modeled as a system of two ordinary differential equations. The model used for each population is a Morris–Lecar system (Morris and Lecar 1981). This is a common model for excitable systems and is a reduced version of the Hodgkin–Huxley model (Hodgkin and Huxley 1952). Note that using a Morris– Lecar system to model the activity of populations of cells rather than individual cells is not a typical application of this mathematical system. This model was originally developed to model the voltage of an individual cell, but we are using it here to model the activity of distinct populations of cells. The equations for AMIN can be written as δ A x !A = f A (x A , y A ) − I VLPO + I ORX + I0A − I HOM + InA

(1)

y !A = g A (x A , y A ), while the equations for VLPO can be written as δV x V! = f V (x V , yV ) − I AMIN − I SCN + I0V + I HOM + InV

(2)

yV! = gV (x V , yV ).

Here, x A and x V represent the overall population activity of the wake-promoting cells and sleep-promoting cells, respectively, and δ A and δV are constants. The variables y A and yV are recovery variables. The nonlinear functions f and g are of the form f (x, y) = 3x − x 3 + 2 − y and g(x, y) = ε(γ H∞ (x) − y)/τ (x) where H∞ (x) = 21 tanh(x/0.01) is a smooth approximation of the Heaviside stepfunction; that is, H∞ (x) ≈ 0 if x < 0 and H∞ (x) ≈ 1 if x > 0. Moreover, τ (x) is a step function that is of the form τ (x) = τ1 + (τ2 − τ1 )H∞ (x). All parameter values can be found in Table 1. The constants I0A and I0V represent background cortical drives and InA and InV are noise terms included to test the robustness of the model. They are described in more detail in the appendix. The terms I VLPO , I AMIN , I SCN and I ORX correspond to inputs from VLPO, AMIN, SCN and ORX, respectively. We assume that I VLPO = gvlpo H∞ (x V ),

I AMIN = gamin H∞ (x A ),

I SCN = gscn C(t)

(3)

where gvlpo , gamin and gscn are constants and C(t) is the circadian pacemaker which we model as in Achermann and Borbély (1994). The details of the circadian input to the model can be found in the appendix. To simplify the model, we do not include a separate differential equation for ORX. Instead, we assume that I ORX = I SCN (1 − H∞ (x V )). This implies that ORX is silent while VLPO is active; on the other hand, if VLPO is silent, then ORX follows SCN.

123

620

M. J. Rempe et al.

Finally, I HOM represents an input corresponding to the sleep/wake homeostat. We assume that this is of the form I HOM = ghom h(t)

(4)

where ghom is a constant and h(t) decays exponentially while the system is asleep and increases while awake. We scale the variables so that the system is defined to be asleep if x A < 0 and awake if x A > 0. Then h(t) satisfies h ! = αh (h max − h) while awake (when x A > 0) and h ! = −βh h while asleep (when x A < 0). Here, αh , βh and h max are constants. Note that while awake the homeostatic inhibition of AMIN grows stronger and the homeostatic drive to VLPO increases. That is, while awake there is an increasing pressure to fall asleep. During sleep the opposite is true. We note that our simple model for each cell group, without inputs, shares many properties of well-known neuronal models including the Morris–Lecar model and the FitzHugh–Nagumo equations (FitzHugh 1961; Morris and Lecar 1981) and (Nagumo et al. 1962). We chose a simple form for the nonlinear functions f and g to simplify the analysis and our discussion of how parameters should be chosen to generate the desired behavior. 2.2 The REM-NREM flip–flop model We now discuss the REM-on and REM-off cell groups of the model introduced in the previous section. Notice in Fig. 1 that the REM-on and REM-off cell groups mutually inhibit each other. We have now divided the VLPO into two cell groups: the clustered VLPO (clVLPO) and the extended VLPO (eVLPO). The eVLPO inhibits REM-off neurons and AMIN inhibits the REM-on neurons. This is consistent with experimental results (Lu et al. 2002). We further assume that clVLPO inhibits eVLPO. This assumption will be discussed in detail later. The REM-on and REM-off cell groups are modeled in a manner very similar to before. The equations for the REM-on cells are δ R x !R = σ

!

f R (x R , y R ) − I AMIN − INREM + I0R + InR

y !R = σ (g R (x R , y R )) ,

"

(5)

"

(6)

while the equations for the REM-off populations are δ N x N! = σ

!

f N (x N , y N ) − IeVLPO − IREM + I0N + InN

y N! = σ (g N (x N , y N )) .

Here, δ R , δ N , I0R and I0N are constants, InR and InN are noise terms, I AMIN is as defined in (3), except the constant gamin is different, and IeVLPO represents input from eVLPO. The parameter σ determines the frequency of the REM/NREM cycling. The

123

Mathematical model of sleep

621

clVLPO is modeled precisely as in Eq. 2; in fact, we assume that (x V , yV ) now represents the overall population activity and recovery of clVLPO. We model eVLPO by the scalar equation xe! = −xe + (ce − ae h e (x V + cv , k xv ) − be h e (x A , k x ))

(7)

where xe represents the overall activity of eVLPO neurons, h e (v, k) =

1 tanh(v/k) 2

and ae , be , ce , cv , k xv and k x are constants. The variable h e depends on the activity of the wake-promoting cells in a step-like manner; however, its dependence on clVLPO is more graded. Finally, we let IeVLPO = gevlpo xe where gevlpo is a constant. Note that xe receives inhibition from clVLPO. This is an assumption that we have made that was necessary to faithfully reproduce correct REM/NREM dynamics. This is discussed in greater detail in Sect. 5. Note that neither the laterodorsal tegemental nucleus (LDT) or the pedunculopontine tegmental nucleus (PPT) are included in our model of the REM flip–flop as they are mentioned, but not included, in the original mutually inhibitory flip–flop switch (Lu et al. 2006).

3 Results 3.1 Timing of sleep, wake, and ultradian dynamics A solution of the model under normal conditions is shown in Fig. 2a. Note that the system sleeps for approximately 8 h and is awake for approximately 16 h with rapid transitions between sleep and wakefulness. Sleep occurs during the trough of the circadian cycle. Figure 2b shows the activity of the REM-on and REM-off cell groups during one sleep episode. The model reproduces several known properties of the ultradian rhythm (Swick 2005): sleep begins with an episode of NREM activity and there are then several transitions between REM and NREM sleep. Moreover, the NREM episodes become increasingly shorter and REM episodes increasingly longer while asleep so that late sleep is dominated by REM activity and the system wakes up from REM sleep. The model parameters have been selected so as to match the desired dynamics. The parameter values are shown in Table 1. Experimental data have shown that aminergic nuclei are most active during wakefulness, are essentially silent during REM sleep, but do show some activity during NREM sleep (Saper et al. 2001). Our modeling approach does not readily achieve intermediate levels of activity, so the aminergic group is considered inactive during both NREM and REM sleep. The transitions between sleep and wake and between NREM sleep and REM sleep take 1–2 min. To switch from sleep to wake, the VLPO shuts down 1–2 min before

123

622

M. J. Rempe et al.

(A)

2 1

x,xV

xV x

0 −1 −2

Waking

Sleep

Waking

Sleep

C(t) 12AM

8AM

4PM

12AM

8AM

4PM

12AM

time (hours)

(B) 2

x xN xR

1

x,x R,xN 0 −1 −2

8PM

12AM

4AM

8AM

12PM

(C)

2

x,xV,xN,xR

time (hours)

1

x xV xR xN

0 −1 −2 8AM

8PM

8AM

8PM

8AM

time (hours) Fig. 2 a Normal sleeping and waking. The solid trace represents the wake-active population (AMIN), and the dashed trace represents the sleep-active population (VLPO). b REM/NREM cycling during one sleep episode. The thin dashed trace is the REM-off population and the solid trace is the REM-active population. c Without orexin, there are more frequent transitions into and out of sleep, and more REM episodes than NREM episodes, including several REM episodes during waking

the AMIN turns on and the switch from NREM to REM occurs when the NREM population shuts down 1–2 min before the REM population activates. When the system transitions from wake to sleep, the VLPO activates first and shuts down the aminergic group after 1–2 min. Similarly, the transition from REM to NREM occurs when the NREM population activates and shuts down the REM population. The in between time when one population inactivates and the other activates is considered as a transition zone, which is relatively short on the time-scale of the overall sleep structure. This is consistent with the flip–flop hypothesis and is discussed more in Sects. 4.3 and 4.6.

123

Mathematical model of sleep

623

3.2 Removal of orexin leads to unstable switching With the input from ORX to AMIN removed the model reproduces a few key characteristics of the narcolepsy phenotype (Fig. 2c). Note that AMIN may cycle several times during which time VLPO is silent. We interpret this as several short wake episodes interrupted by very brief sleep episodes. This behavior is consistent with the narcoleptic phenotype which is characterized by excessive daytime sleepiness often leading to several short daytime sleep episodes (Dauvilliers et al. 2007). This is also consistent with the hypothesis proposed by Saper and colleagues: that orexin neurons act to stabilize the sleep/wake switch. Also note that overall there are more REM sleep episodes than NREM sleep episodes. This is consistent with the REM behavior of narcoleptics, who have inappropriate intrusions of REM sleep, presumably because without orexin the REM/ NREM balance has been tipped in favor of REM-on firing (Boeve et al. 2007). The model makes a suggestion about how that balance is tipped: if VLPO does not activate, it does not inhibit eVLPO so eVLPO is free to inhibit the REM-off group which has the effect of disinhibiting the REM-on group. The model also demonstrates REM episodes that intrude on waking and more sleep-onset REM (SOREM) than in the normal case, which is also consistent with the narcolepsy phenotype (Krahn et al. 2001). It is not our intention to model cataplexy, which frequently accompanies narcolepsy. It has been reported that monoaminergic nuclei become discoordinated during cataplexy (John et al. 2004; Wu et al. 1999) so modeling this behavior is beyond the scope of this model since it treats the monoaminergic nuclei as a single population. Though only reproducing a few characteristics of narcolepsy, it is important to note that although considerable inter- and intra-individual variation have been observed in the degree and nature of symptoms, these two features have been observed repeatedly in narcoleptic subjects and are considered indicative of the disorder.

3.3 The two-process model and sleep deprivation To compare the output of our model with that of the two-process model we plotted I HOM , the input corresponding to the sleep/wake homeostat, along with two other curves. These curves are defined as L ≡ I SC N + K L and U ≡ I SC N + K U where K L and K U are constants. We will explain how these constants are chosen in the analysis section. The system wakes up when I HOM reaches the lower threshold L and falls asleep when I HOM reaches the upper threshold U , consistent with the two-process model of sleep regulation (Daan et al. 1984). We used the model to simulate sleep deprivation and the results were identical to those of the two-process model. We simulated keeping the system awake in two-hour increments for up to 24 h past the usual bedtime. This was accomplished by increasing the cortical drive I0A to the wake-promoting cells AMIN. We then returned I0A to its default value and this allows the system to fall asleep. The network then wakes up when I HOM crosses L and then falls asleep when I HOM crosses U .

123

624

M. J. Rempe et al.

15 0 540 480 420 360

37.0

0

90 180 270

0

90

180 270

CIRCADIAN PHASE (degrees)

0

30 (%)

SOREMS

30

36.5

model

(B)

experiment

CIRCADIAN INPUT TOTAL SLEEP TIME (normalized) (min)

BODY TEMPERATURE TOTAL SLEEP TIME SOREMS (degree) (%) (min

(A)

15 0 600 400 200 0 1 0.5 0

0

90 180 270

0

90

180 270

0

CIRCADIAN PHASE (degrees)

Fig. 3 The model exhibits a circadian dependence of total sleep time and sleep-onset REM sleep (SOREMS). a Data from human sleep experiments. Upper panel total sleep time plotted as a function of circadian phase at midsleep. Middle panel the circadian distribution of SOREMS. Lower panel circadian wave form of core body temperature. b The same sleep measures as computed using our model. In the lower panel, the trough of the circadian input was assigned phase 0. a Copyright 1995 by Society for Neuroscience. Reproduced with permission of the Society for Neuroscience via Copyright Clearance Center

3.4 Circadian modulation of total sleep time and sleep-onset REM One purpose of the model and of many sleep experiments is to investigate how observed features of sleep arise from the circadian and homeostatic processes and from their interactions. In the laboratory, circadian and homeostatic modulatory effects can be teased apart with a “forced desynchrony protocol” (Kleitman 1987) that involves isolating individuals from time cues while imposing a rest-activity schedule so different from the near-24-h intrinsic circadian period (Campbell et al. 2000; Czeisler et al. 1999; Wright et al. 2001) that subjects are unable to entrain. As rhythms driven by the circadian pacemaker proceed with a period differing from that of the sleep–wake cycle, phase-resetting stimuli are distributed uniformly across the circadian cycle, avoiding feedback-resetting effects (Campbell et al. 2000). Dijk and Czeisler (1995) observed eight men subjected to a 28-h rest-activity cycle, studying the extent to which circadian and homeostatic processes separately and together influence several aspects of sleep. Here we describe how output from the model compares with their results. In order to replicate this experiment with the model, we shifted the circadian curve by 15-min increments and started each simulation with the homeostat at the value it takes on after a normal night of sleep. Noting that Dijk and Czeisler’s protocol results in 18 h 40 min of wakefulness prior to each sleep episode, we temporarily increased the cortical drive I0A to maintain wakefulness for a similar period before withdrawing the cortical drive and allowing an opportunity for sleep. Experimental measurements of total sleep time revealed a broad peak encompassing circadian phase 0 (the trough), having a rapid rise and more gradual decline (Dijk and Czeisler 1995). As seen in Fig. 3, these features are also reproduced in the model.

123

Mathematical model of sleep

625

Total sleep time has a gradual decline and a rapid rise (top panels) with a peak corresponding to the trough of the circadian cycle (bottom panels). Dijk and Czeisler found that sleep tended to be undisturbed in the early part of sleep episodes regardless of the circadian phase at sleep initiation (Dijk and Czeisler 1994); however, there was a circadian dependence of sleep efficiency (proportion of time in bed spent asleep) during the latter half of the scheduled sleep episodes. Accordingly they plotted total sleep time versus circadian phase at midsleep, and we have done the same in order to facilitate comparison. The results are consistent with the relation found by Strogatz (1986). Our model output differs notably from the data of Dijk and Czeisler in that there are phases at which the model does not go to sleep. The circadian phases at which the model failed to sleep coincide with phases previously described as a wake-maintenance zone (Strogatz 1986) or forbidden zone for sleep (Lavie 1986); unentrained, free-running subjects rarely choose to initiate sleep during this circadian phase (Strogatz 1986). The subjects in Dijk and Czeisler’s study, instructed to try to sleep during this zone, were able to do so albeit with a pronounced increase in sleep latency. Both model and data share the property that total sleep time versus circadian phase results in a left-skewed curve. Dijk and Czeisler’s study also revealed a pronounced circadian modulation of REM sleep, with a peak REM sleep propensity shortly after the minimum of the circadian cycle (Dijk and Czeisler 1995). During entrained conditions, this peak falls nearly two hours before habitual waking and is reflected in the observation that the amount of REM sleep increases across the sleep episode. An additional expression of this modulation arises in the forced desynchrony protocol in terms of the likelihood, with respect to circadian phase, of experiencing sleep-onset REM sleep (SOREMS). Dijk and Czeisler (1995) found that the percentage of sleep episodes with SOREMS (defined as REM sleep latency of at most 10 min instead of the more usual 70–80 min or longer (McCarley 2007; Plazzi et al. 2008) showed a broad peak beginning shortly before circadian phase 0 and preceded by an expanse of circadian phases during which sleep rarely if ever began with REM sleep. There is disagreement between labs as to what constitutes “sleep onset” and therefore what counts as SOREMS. We have used 25 min as the cutoff for SOREMS, which is consistent with more recent sleep studies (Sasaki et al. 2000; Takeuchi et al. 2003) which chose a 25-min cutoff for scoring SOREMS because of the observed bi-modal distribution of sleep-onset REM periods (SOREMP) (Bes et al. 1996; Sasaki et al. 2000) and a 25-min cut-off was able to capture both modes. Computing the frequency of SOREMS in our model, we find a pattern similar to that in the Dijk and Czeisler study (Fig. 3b). Note that in this figure the data are double-plotted to illustrate the circadian dependence and to be consistent with the presentation of the experimental data. In our model SOREMS occurs at roughly the same frequency and at the same circadian phase as seen experimentally. In the model, the occurrence of SOREMS is dependent on the circadian phase at sleep onset because of the following mechanism: if a sleep episode begins when the circadian input is low, sleep is initiated by VLPO activating, and not by AMIN inactivating. However, when circadian input is high, the AMIN cubic is low, so the fixed point in the active state is close to the slow nullcline. With the addition of noise the cubic nullcline may shift downward, causing this fixed point to lose stability and jump down. When sleep begins with this mechanism AMIN releases REM from

123

626

M. J. Rempe et al.

Fig. 4 Phase plane diagrams of a relaxation oscillator for three different values of a stimulating input

inhibition allowing it to fire right away causing SOREMS. As a result, noise is crucial for SOREMS to occur in this model. When sleep begins with VLPO jumping up, eVLPO releases the REM-off population before the REM-on population has time to activate. 4 Analysis We analyze solutions of the model using phase plane methods. The analysis helps to explain the mathematical mechanisms at work in the model and leads to insights and new hypotheses for the biological mechanisms underlying these rhythms. Moreover, the analysis is extremely useful in demonstrating how parameters in the model should be chosen to obtain the desired behaviors. The analysis is described in a series of model equations of increasing complexity. We first introduce the basic phase-plane approach by considering one population of cells, with constant, and then time-varying, input. We then demonstrate how these methods apply to multiple populations by considering two mutually inhibitory nuclei with constant inputs. Finally, we consider the full models for the sleep/wake dynamics. 4.1 One nucleus without coupling We analyze solutions using phase plane methods. We begin by considering one population of cells with a constant input, ITOT , without noise. This is modeled as x ! = f (x, y) + ITOT y ! = g(x, y).

(8)

Figure 4 shows the phase plane for three different values of ITOT . Note that there exist I1 < I2 so that if I < I1 then there exist three fixed points as shown in Fig. 4a. It is not hard to show that the fixed point along the left branch of the cubic x-nullcline is stable, while the other fixed points are unstable. If I1 < ITOT < I2 , then there is a unique fixed point that lies along the middle branch of the x-nullcline and it is unstable. In this case, (8) exhibits a stable limit cycle. Finally, if ITOT > I2 , then there

123

Mathematical model of sleep

627

Fig. 5 Phase planes when there is time-varying input. A stable fixed point may disappear by a saddle-node bifurcation as shown in b

are again three fixed points as shown in Fig. 4c. The fixed point that lies along the right branch of the x-nullcline is stable. In what follows, we will say that a cell-group (x, y) is silent if x < 0 and is active if x > 0. We next consider one population of cells with time-varying input, ITOT = ITOT (t). Suppose that ITOT (0) > I2 and ITOT (t) is a decreasing function. When t = 0, the cell-group is active and lies at the stable fixed point along the right branch of the x-nullcline. This is shown in Fig. 5a. With increasing t, ITOT decreases and this lowers the x-nullcline. We assume that ITOT (t) changes slowly so that the solution tracks near the point where the y-nullcline intersects the right branch of the x-nullcline, as long as this point exists; that is, as long as the “right knee” of the x-nullcline lies above the y-nullcline. When ITOT (t) = I2 , the right knee of the x-nullcline touches the y-nullcline as shown in Fig. 5b. At this time, the stable fixed point along the right branch of the x-nullcline comes together with the unstable fixed point along the middle branch; this corresponds to a saddle-node bifurcation for the system in which ITOT is constant. As t increases further, these two points disappear and the solution of (8) is forced to “jump down” to the silent phase, as shown in Fig. 5c. As we shall see, this describes the basic mechanism for how the full system “falls asleep” and “wakes up.” 4.2 The network without circadian and homeostatic input Now consider a network in which VLPO and AMIN mutually inhibit each other, but there is no circadian or homeostatic input. The equations are then: x !A = f A (x A , y A ) − gvlpo H∞ (x V ) + I0A y !A = g A (x A , y A )

x V! = f V (x V , yV ) − gamin H∞ (x A ) + I0V yV! = gV (x V , yV ).

(9)

We consider the separate phase planes associated with each cell-group. In fact, there are two phase planes associated with each cell-group depending on whether that cellgroup receives inhibitory input from the other cell group. Figure 6a shows the phase plane associated with VLPO. There are two x V -nullclines; these depend on whether

123

628

M. J. Rempe et al.

Fig. 6 a The nullclines for the VLPO system while asleep and while awake. b The nullclines for the AMIN system while asleep and while awake

AMIN is active or not; that is, whether x A > 0 or x A < 0. Note that if AMIN is active, then VLPO is silent; this is because the stable fixed point corresponding to VLPO lies on the left branch of its cubic nullcline. On the other hand, if AMIN is silent, then VLPO is active, because the corresponding stable fixed point now lies on the right branch of the cubic nullcline. Figure 6b shows the phase plane associated with AMIN. As above, there are two x A -nullclines depending on whether VLPO is active or not. If VLPO is active, then AMIN is silent, while if VLPO is silent, then AMIN is active. Note that there are two stable fixed points of the full (four-dimensional) system (9). At one of these fixed points, VLPO is silent and AMIN is active; that is, the system is “awake.” At the second stable fixed point, VLPO is active and AMIN is silent, so the system is “asleep.” 4.3 The sleep/wake switch We now consider the full system (1), (2) and step through the complete sleep/wake cycle. We begin with the system “asleep” so that VLPO lies at a stable fixed point along the right branch of the x V -nullcline and AMIN lies at a stable fixed point along the left branch of the x A -nullcline. Now the circadian input C(t) and the homeostat h(t) will change and, while asleep, the total input to VLPO will decrease. As described earlier, this has the effect of lowering the x V -nullcline. Eventually there is a saddle-node bifurcation when the right knee of the x V -nullcline just touches the yV -nullcline. At this time, the system “wakes up,” by which we mean the following: Immediately after the saddle-node bifurcation, (x V , yV ) jumps down to the silent phase as shown in Fig. 7a. This then “releases” AMIN from inhibition. That is, IVLPO = gvlpo H∞ (x V ) switches from gvlpo to 0. Therefore, the x A -nullcline also changes as shown in Fig. 7b. When IVLPO = 0, the x A and y A nullclines do not intersect along the left branch of the x A -nullcline; however, they do intersect along the right branch of the x A -nullcline. Hence, (x A , y A ) must jump up to the active phase and approach the stable fixed point along the right branch of the x A -nullcline. Once (x V , yV ) jumps down the process

123

Mathematical model of sleep

629

Fig. 7 The mechanism for waking up. a Waking up is initiated when there is a saddle-node bifurcation in the active phase of the VLPO system. When VLPO jumps down, this releases AMIN from inhibition and it jumps as shown in b

(A)

(B)

MIN LPO

Fig. 8 The mechanism for falling asleep. a The system falls asleep when there is a saddle-node bifurcation of the VLPO system and it jumps up. The resulting inhibition onto AMIN lowers the AMIN cubic nullcline, forcing AMIN to jump down as shown in b

of AMIN being released from inhibition and transitioning to its active phase takes 1–2 min, so the transition from sleep to wake is rapid on the timescale of the overall sleep structure. These 1–2 min of transition are a gradual evolution between the states, which is consistent with the idea that sleep/wake transitions are on the order of minutes (Merica and Fortune 2004). We continue to step through the sleep/wake cycle, now assuming that the system is awake. We note that while awake, the total input to VLPO will increase. Similar to above, this has the effect of raising the x V -nullcline. Eventually there is a saddle-node bifurcation when the left knee of the x V -nullcline just touches the yV -nullcline as shown in Fig. 8a. Immediately after this, (x V , yV ) jumps up to the active phase, this inhibits AMIN and (x A , y A ) jumps down as shown in Fig. 8b. That is, the system falls asleep. In this transition, once VLPO activates, it takes about 1–2 min before AMIN jumps down to its inactive state, giving a transition time of 1–2 min for falling

123

630

M. J. Rempe et al.

asleep. Again, this is rapid on the timescale of the sleep/wake pattern, but the process itself is somewhat gradual on the timescale of minutes, but rapid in comparison to the time spent in each state. This transition time is also consistent with human sleep experiments (Merica and Fortune 2004). Note that in our model, VLPO is, in some sense, responsible for initiating both sleep and wake. That is, sleep begins when VLPO jumps up to the active phase and this forces AMIN to jump down to the silent phase; the system wakes up when VLPO jumps back to the silent phase and this releases AMIN from inhibition. This mechanism is consistent with results that suggest that the VLPO drives the transition from wake to sleep (McGinty and Szymusiak 2003; Szymusiak et al. 2007). 4.4 Relationship to the two-process model Note that the system wakes up and falls asleep when there are saddle-node bifurcations of the reduced system (8) in which the total input is taken to be a constant bifurcation parameter. Since the total inputs to VLPO and AMIN involve both a circadian and a homeostatic component, the position of the saddle-node bifurcations defines a relationship between the circadian and homeostatic inputs when the system either wakes up or falls asleep. It is this relationship that leads to a direct correspondence between the model presented here and the two-process model. Recall that the system falls asleep when VLPO exhibits a saddle-node bifurcation at its left knee and the system wakes up when VLPO exhibits a saddle-node bifurcation at its right knee. Suppose that these saddle-node bifurcations occur when the total input ITOT V to VLPO are ITOT V = I1V and ITOT V = I2V , respectively. Now, ITOT V = I0V + I HOM − I SCN . Hence, the system falls asleep when I HOM (t) = I SCN (t) − I0V + I1V and the system wakes up when I HOM (t) = I SCN (t) − I0V + I2V . To plot the homeostat and two circadian thresholds in the two-process style, the lower threshold L is the curve defined by the function I SCN (t) − I0V + I2V ≡ gscn C(t) + K L and the upper threshold U is the curve defined by the function I SCN (t) − I0V + I1V ≡ gscn C(t) + K U .

123

Mathematical model of sleep

(A)

631

(B)

VLPO

6 4

4

IA=0

yV

yA

IA=gA

2 0

IV=0 with ORX

AMIN

6

sleep

IV=0 no ORX

2

IV=gV 0

wake -2

0

-1

xV

1

2

sleep -2

-1

0

xA

1

2

(C) 5

Ihom

4 3

*

2 1 8AM

** 4PM

12AM

8AM

time

4PM

12AM

8AM

Fig. 9 The loss of orexin causes more frequent sleep/wake transitions. a As in Fig. 8 the system may fall asleep when there is a saddle-node bifurcation of the VLPO system. g A = gamin . b Without excitatory ORX input, there is no longer a stable fixed point of AMIN corresponding to the wake state and the system falls asleep when AMIN reaches the right knee of its corresponding cubic. gV = gvlpo . c The graph of I HOM (t) for the solution shown in Fig. 2c

By choosing L and U in this manner, the homeostat changes direction when it reaches L or U . Note that this approach is similar in spirit to the approach used by Phillips and Robinson (2008) to show that their model is consistent with the two-process model. However, they do not explicitly show how to compute the appropriate offsets for the circadian curves.

4.5 Orexin neurons help to stabilize the flip–flop switch Figure 2c demonstrates that if we remove the excitatory input from the orexin neurons to AMIN, then the model exhibits more frequent transitions between the sleep and wake states. Phase plane analysis is very useful for understanding the mechanism underlying this behavior. Consider the phase plane corresponding to AMIN while the system is awake. It is useful to compare the phase planes of AMIN with and without input from ORX. Recall that, while awake, I ORX = I SCN ; hence, removing I ORX is equivalent to removing excitatory input from SCN. Figure 9b shows the phase plane of AMIN with and without excitatory input from ORX. Note that with ORX input, the x A -nullcline intersects the y A nullcline at a stable fixed point along the right branch of the x A -nullcline. This corresponds to the wake state. Removing excitatory ORX input lowers the x A -nullcline. If parameters are

123

632

M. J. Rempe et al.

chosen appropriately, then there will no longer be a stable fixed point representing the wake state. Now consider what happens when the system “wakes up.” Recall that the waking up process begins when there is a saddle-node bifurcation of the VLPO system and (x V , yV ) jumps down to the silent phase as shown in Fig. 9a. This then releases AMIN from inhibition and (x A , y A ) jumps up to the active phase. The system is now awake and (x A , y A ) evolves up the right branch of the x A -nullcline. If there were excitatory input from ORX, then (x A , y A ) would approach a stable fixed point in the wake state (see Fig. 7). However, without ORX input there is no longer such a stable fixed point. In this case, (x A , y A ) is forced to jump down to the silent phase when (x A , y A ) reaches the right knee of the x A -nullcline, at which time the system falls asleep. If (x A , y A ) evolves quickly along the right branch of the x A -nullcline, then the system remains awake for a short time. This explains the rapid transitions from the wake to the sleep state. The mechanism for the rapid transitions from the sleep to the wake state is more subtle. Note that VLPO does not receive input from ORX, so we cannot explain these rapid transitions by considering how removing ORX changes the VLPO nullclines. There are actually two reasons why removal of ORX leads to short sleep episodes. The first is because the homeostat does not have enough time to build up while awake. Perhaps the easiest way to understand this is to consider Fig. 9c which shows I HOM (t) for the solution shown in Fig. 2c. Consider the location marked with ∗. Note that, as before, the system wakes up when I HOM touches the lower threshold L. This is when there is a saddle-node bifurcation in the active phase of the VLPO (x V , yV ) system. While awake, I HOM builds up. However, because the system is awake for only a short time, I HOM builds up by only a small amount. Then the system falls back asleep, I HOM decreases and the system remains asleep until I HOM returns to L. Because I HOM increased only a small amount while awake, the duration of this sleep episode must be small. Note that, in our model, the duration of the sleep episode depends on whether it occurs along the rising or falling phase of the circadian rhythm. As illustrated in Fig. 9c, the sleep episode will be shorter if it occurs along the rising phase of the circadian rhythm and will be longer if it occurs along the falling phase of the circadian rhythm. The second reason why removal of ORX leads to short sleep episodes can be understood by examining Fig. 2c. Note that around 4 PM the second day the AMIN population activates several times while VLPO remains silent. Comparing Figs. 2c and 9c, we see that this AMIN cycling occurs during sustained, elevated levels of the circadian input I SCN . This location is marked in Fig. 9c with ∗∗. Since I SCN inhibits VLPO, when I SCN is elevated, there is increased inhibition to VLPO. This increased inhibition shifts the (x V , yV ) cubic nullcline downward. Because the (x V , yV ) cubic nullcline is lower than normal, when (x A , y A ) jumps down from its active state and releases VLPO, the left knee of the VLPO cubic is still below the slow nullcline. Therefore, (x V , yV ) remains at an inactive fixed point instead of jumping up to the active phase. Conversely, I SCN has an excitatory influence on AMIN through the I ORX term. Therefore, when I SCN is at an elevated level, the AMIN cubic nullcline is shifted upward such that both knees are between the upper and lower values of the slow

123

Mathematical model of sleep

633

(B1) yN

(B2)

NREM

6

6

4

yR

2

−2 −1 0

0.2 0

3

2

0

4

yN 8

12

time (hours)

16

20

1

−2 −1

(C2) 6

4

yR

2

−2 −1

(D1) yN

4

0

xR

1

2

1

2

REM

4 2 0

0

6

2

2

NREM

6

sleep bout

1

x e 0.6

xN

(C1)

1

4

0

0

(A)

REM

0

xN

1

2

−2 −1

(D2)

NREM

6

IREM=0

yR

2

IREM=gREM

0 −2 −1

0

xN

1

2

4

0

xR

REM INREM=0

2 INREM=gNREM

0 −2 −1

0

xR

1

2

Fig. 10 a The activity of eVLPO in the model throughout a sleep bout. b The nullclines for NREM b1 and REM b2 corresponding to the location marked 1 in panel a. c The nullclines for NREM c1 and REM c2 corresponding to the location marked 2 in panel a. d The nullclines for NREM d1 and REM d2 corresponding to the location marked 3 in panel a

nullcline as in Fig. 4b. As a result, (x A , y A ) makes several oscillations. If we interpret the system to be awake whenever AMIN is active, then these oscillations represent short wake bouts separated by very short sleep episodes, even though VLPO is not able to become active. 4.6 The REM-NREM switch To understand how the model demonstrates the REM/NREM rhythm, we first consider the behavior of xe , the eVLPO activity. This is shown in Fig. 10a. Note that xe is relatively high during wake (label 1), falls just before sleep onset (label 2) and then grows again as the sleep bout progresses (label 3). The reason why xe decreases just before sleep onset is because there is an increase in xv , the clVLPO activity, during this time (see Fig. 2a) and we assume that clVLPO inhibits eVLPO. While asleep, xv slowly decreases and this causes xe to increase. While awake, both the NREM and REM systems are silent and lie at stable fixed points as shown in Fig. 10b. During this time, REM receives strong inhibitory input from the active AMIN cells, while NREM receives strong inhibitory input from eVLPO. Just before sleep onset, xe drops, decreasing the inhibitory input to NREM. This causes the NREM cubic nullcline to rise. Once the left knee of the NREM cubic

123

634

M. J. Rempe et al.

rises above the slow nullcline, the stable fixed point disappears via a saddle-node bifurcation and (x N , y N ) jumps up to the active phase. Through this mechanism, the REM-off population activates just before sleep onset and remains active during the initial portion of the sleep episode. Note that when sleep begins, the wake-active cells, AMIN, stop firing and, therefore, stop sending inhibition to the REM-active cells. However, the REM-active cells do not rebound and begin firing at sleep onset because of the inhibition they receive from the already active NREM cells. This is illustrated in Fig. 10c. While asleep, the NREM and REM cell populations take turns firing. A NREM episode ends when (x N , y N ) reaches the right knee of its corresponding cubic and jumps back to the silent phase. This releases the REM cells so that (x R , y R ) jumps up to the active phase. Since the inhibition from REM to NREM is weaker than the inhibition from NREM to REM, when REM activates, it inhibits NREM, shifting the NREM cubic nullcline down, but not enough to create a stable fixed point there. So NREM cycles around its cubic and drives the REM population by post-inhibitory rebound to oscillate as well. Note that the NREM episodes become shorter and REM episodes become longer while asleep. The reason why this happens is because during this time, eVLPO activity slowly rises. This increases inhibition onto NREM which, in turn, lowers the NREM cubic nullcline. When NREM is active, its solution moves up the active branch of this nullcline, but this nullcline itself is moving downward. This results in NREM spending less time in the active phase. The REM/NREM cycling continues until the system wakes up. Once AMIN activates, it inhibits the REM-active group, causing the REM-active population to stop oscillating. Note that the frequency of the REM/NREM cycling is dictated by the parameter σ , which can be tuned to give more or fewer cycles between REM and NREM. Similar to the sleep/wake switch, the transitions from NREM to REM and from REM to NREM are driven by the NREM population. When NREM inactivates it releases the REM population, allowing (x R , y R ) to move toward the stable fixed point along the right branch of the x R -nullcline. Once (x N , y N ) jumps down it takes 1–2 min for (x R , y R ) to reach its active phase. In the transition from REM to NREM, the (x N , y N ) jumps up to its active phase first and then inhibits REM, causing it to jump down to its inactive phase. Once NREM activates, it takes about 1–2 min for REM to jump down to its inactive phase. So each of the transitions between sleep and wake and between REM and NREM sleep take less than 3 min, which is consistent with the rapid transitions implied by the flip–flop conceptual model. 4.7 The role of noise Including a noise term in each of the models of the four main cell groups does not qualitatively alter the sleep–wake behavior and we conclude that the behavior of the model is robust to noise. However, introducing noise is an important element in making the model correctly reproduce the SOREM data; without noise the model does not exhibit SOREM at any circadian phase. Without noise, sleep is always initiated by

123

Mathematical model of sleep

635

(B)

(A)

high circadian input (SOREM less likely)

low circadian input (SOREM more likely)

6

6

4

4

yA

yA 2

2

0

0 −2

−1

0

xA

1

2

−2

−1

0

xA

1

2

Fig. 11 Noise plays an important role in SOREM. a When circadian input is high, the AMIN cubic is well above the slow nullcline even with the jitter caused by noise (the upper and lower bounds of its location are shown in dashed lines). In this case sleep starts by VLPO activating (as opposed to AMIN inactivating). b When circadian input is low, the addition of noise may cause the AMIN cubic curve to fall below the slow nullcline (lower dashed line). In this case, sleep would initiate by AMIN jumping down rather than VLPO jumping up. The mechanism by which sleep begins determines if SOREM occurs or not (see text for explanation)

VLPO becoming active and, as described previously, this allows eVLPO activity to drop just before sleep-onset. This, in turn, releases NREM from inhibition and sleep begins with NREM. Since the noise terms are added to the x ! equation for each of the four cell groups, they have the effect of shifting the cubic nullclines up or down by a small amount each time step. Although without noise sleep is always initiated by the VLPO group activating and not by the AMIN group inactivating, with noise sleep may be initiated by AMIN jumping down from its active phase. This is because when circadian input is low the right knee of the AMIN cubic is close to the slow nullcline. With the addition of noise, the AMIN cubic may fall below the slow nullcline allowing AMIN to deactivate (see Fig. 11). When AMIN jumps down, it releases REM from inhibition before the eVLPO has time to release NREM from inhibition, so REM activates before NREM. When the circadian input is high, the AMIN cubic nullcline is not close enough to the slow nullcline for this mechanism to work. In this manner the model exhibits SOREM only at certain circadian phases and then only a certain percentage of the time. 5 Discussion We have presented a model that accounts for several features of the human sleep/wake cycle. These features include the timing of sleep and wakefulness under normal and sleep-deprived conditions, ultradian rhythms, rapid switching between sleep and wakefulness due to the loss of orexin, and the circadian dependence of total sleep time and sleep-onset REM. The model demonstrates how these features depend on

123

636

M. J. Rempe et al.

interactions between a circadian pacemaker and a sleep homeostat and shows consistency between the two-process model for sleep regulation and the flip–flop models. A major goal of this study is to determine whether the neuronal components in the flip–flop models, with the addition of a sleep-homeostatic process and a circadian pacemaker, are sufficient to account for these features of the human sleep/wake cycle. The model is based on flip–flop models for sleep/wake and REM/NREM due to Saper and colleagues. It is a minimal model in the sense that, besides the sleep homeostat and constant cortical drives, the model includes only those nuclei described in the flip–flop models. Each of the cell groups is modeled by at most two differential equations for the evolution of the total population activity. Using a simple cubic-nullcline in the model for each cell group allows us to more easily choose parameters based on the analysis presented in this paper. The mathematical model used for each of the four main cells groups is a model that was originally developed to model the voltage of an individual cell. Here we are using a similar formulation, but with the intention of representing the activity within populations as opposed to the voltage of an individual cell. For each population i, the variable xi represents a normalized average voltage of cells within that population. While this approach has been employed by others to represent populations of cells before (Diniz Behn et al. 2007, 2008; Rubin et al. 2009), it does constitute a level of phenomenology not present in firing rate models or conductance-based models. We chose to use an oscillator model for the cell populations because we wanted to investigate whether typical sleep patterns could arise from the interactions of simple oscillators and because we wanted to use a model that is tractable with phase plane analysis. The synaptic connections are consistent with those described in the flip–flop models. However, in order to account for certain features of the ultradian rhythms, we found it necessary to include one additional hypothesis about the connections: that VLPO cluster inhibits eVLPO. This assumption is needed to ensure that under normal conditions, sleep typically begins with a NREM episode and that later sleep exhibits less NREM activity. The model addresses human sleep/wake behavior and the underlying neural interactions. While both the two-process and the flip–flop models have been applied to human sleep, it is important to note that the two-process model is based on human EEG recordings while the flip–flop models are based on electrophysiology data from rodents. Rodent sleep has a different structure than human sleep since mice are nocturnal and exhibit polyphasic sleep rather than sleep that is primarily consolidated into a single, daily episode. Nevertheless, data from mice and rats have proved useful in making inferences concerning human physiology; for example, genetically orexin-deficient mice express some symptoms found in narcoleptic humans who are also orexin deficient. The early observations of von Economo (1930) of altered sleep– wake behavior of humans with posthumously-confirmed brain lesions are consistent with the hypothalamic switch hypothesis. Nonetheless, species differences have been identified among mice, rats, and cats (Siegel 2008) and equal or greater as yet unidentified species differences surely exist between each of these species and humans. The model developed here, and its ongoing refinements, may be useful in identifying likely differences of human neural circuitry for which direct, human data are not available.

123

Mathematical model of sleep

637

The model presented here differs in a number of ways from other recent mathematical models of sleep/wake cycles which have either not included REM/NREM dynamics (Chou 2003; Phillips and Robinson 2007) or have included other hypothetical interactions (Yin 2007). There have been previous papers motivated, in part, by the hypothalamic sleep switch (Diniz Behn et al. 2007; Phillips and Robinson 2007; Tamakawa et al. 2006). Diniz Behn et al. (2007) use a similar approach to model sleep in the mouse brain. They use a relaxation oscillator of a very similar form to the one used here. Also, they employ phase-plane analysis to analyze behavior, much as we do. They include three rather than four distinct neural populations representing a sleep-, wake-, and REM-promoting nuclei. Their model has the ability to produce polyphasic sleep rhythms with frequent brief awakenings, features more prominent in the sleep of mice than of humans. Orexin appears in their model as a factor scaling the inhibition from wake-active cells to VLPO; however, they do not consider the removal of the orexinergic drive. In contrast to the model presented here, their model does not include a REM-off population or circadian drive, but does include a phenomenological homeostatic REM drive allowing for increased REM activity through the sleep period. A key difference between their model and ours is that our model produces an ultradian oscillation without the addition of an ultradian drive, although this behavior relies on our assumption that clVLPO inhibits eVLPO. This inhibitory connection is a novel prediction of our model. Our method of modeling the orexin system is simpler than that used by Diniz Behn and colleagues in a recent study (Diniz Behn et al. 2008). Like our approach, their model of orexin is only active during wake, but to model the effect of orexin on polyphasic mouse sleep they used a saturating sigmoidal function of the time awake. This ensured that orexin affected long wake bouts, but not short ones. In contrast to the model presented here, in their model orexin is not modeled as a separate population, but is assumed to affect the homeostat, and the strength of connections between the sleep- and wake-active populations. This approach is necessary for modeling rodent, but not human sleep. The model of Phillips and Robinson (2007), based on the hypothalamic switch hypothesis without REM/NREM ultradian dynamics, describes hysteresis in sleep and wake states. Their model does not contain orexinergic cells, but the authors describe a possible interpretation of how orexin could be incorporated. A later study by Phillips and Robinson (2008) demonstrates consistency between their biologically-based model and the two-process model. They plot the output of the two-process model after the circadian component has been subtracted from each of the three components (the upper threshold, the lower threshold, and the homeostat). For the standard two-process model this results in the total sleep drive (the homeostat minus the circadian component) staying between two constant thresholds, one for waking up and one for falling asleep. However, when their model is put in this form the total sleep drive does not stay between the two constant thresholds, since their model includes hysteresis. Our model does not demonstrate hysteresis, and therefore the total sleep drive remains between the two thresholds, like the original two-process model. Some earlier models have relied on McCarley’s hypothesis that REM/NREM cycles are regulated by reciprocal interaction between cholinergic neurons and brainstem

123

638

M. J. Rempe et al.

monoaminergic populations (Ferrillo et al. 2007; Massaquoi and McCarley 1992; McCarley and Massaquoi 1992; McCarley and Massaquoi 1986; Nobili et al. 2001); these models have generally focused on REM/NREM dynamics with little or no detail concerning the generation of sleep–wake behavior. While McCarley’s models may provide a better match to neuronal data in which firing rates change more gradually than predicted by a flip–flop model (McCarley 2007; Siegel 2006), some experimental studies have failed to support these models since lesions of monoaminergic populations have little effect on REM/NREM cycles (Fuller et al. 2007). In contrast, some recent studies implicate GABAergic and glutamatergic interactions in the dynamics of REM sleep (Boissard et al. 2002; Lu et al. 2006) consistent with the REM flip–flop hypothesis. The nature of the sleep/wake homeostat is poorly understood. In our model, we simply assumed that the homeostat, represented by the function h(t), decays exponentially while the system is asleep and increases exponentially towards some limiting value while awake. This model for the sleep/wake homeostasis is similar to that proposed by the original makers of the two-process model of sleep regulation (Achermann and Borbély 2003; Daan et al. 1984). Hence, our model is consistent with the two-process conceptual model, although in our model the levels of the thresholds are determined from sleep- and wake-active populations. This was discussed in Sect. 4.4. There have been numerous papers that have proposed a physiological basis for the sleep/wake homeostasis. For example, Krueger and Obal (1993) (see also Hill and Tononi 2005) have proposed that the homeostasis is a local and use-dependent process that is based primarily at the level of the cortical column. In order to interpret this within the context of our model, we could set the total cortical drives to VLPO and AMIN to be ICTX A = I0A − I HOM and ICTX V = I0V + I HOM . Then, while awake, the total cortical drive to the wake-promoting cells decreases and the total cortical drive to the sleep-promoting cells increases. The opposite holds when the system is asleep. Others have suggested that the sleep/wake homeostasis is associated with the accumulation and decay of adenosine (Porkka-Heiskanen et al. 1997), so another interpretation is that I HOM corresponds to a homeostatic input associated with the accumulation of adenosine. Multiple factors with differing mechanisms likely contribute to sleep homeostasis. Clearly, further research into the nature of the sleep/wake homeostat is needed. Mathematical modeling, including the analysis presented in this paper, should prove very useful in guiding these experiments. Under certain circumstances this model exhibits a phenomena of simultaneously inactive VLPO and monoaminergic groups. In Sects. 3.2 and 4.5 this state is reported within the context of narcolepsy. In Sect. 4.3 this phenomenon is mentioned in the case of a normal transition from sleep to wake. In the narcoleptic case, the authors are not aware of data that report simultaneously inactive AMIN and VLPO, but this would be a hard phenomenon to test. A common tool used to determine whether or not a region is active is fos immunostaining, but since the peak of the synthesis occurs 30–90 min after stimulation and the half-life is on the order of hours (see Deurveilher and Semba 2006) it does not have the resolution to capture these brief moments when both AMIN and VLPO may be simultaneously inactive. Note

123

Mathematical model of sleep

639

that an analogous state (simultaneous activity in AMIN and VLPO) has been observed in another model (Diniz Behn et al. 2007). In the case of the normal transition this “intermediate” state is short (1–2 min) and occurs because when VLPO inactivates it is not assumed that it instantaneously releases the monoaminergic population from inhibition. On the time scale of the sleep/wake structure this transition is fast, but it is not assumed to be instantaneous. In our model the activity of the aminergic cell groups is inactive during both NREM and REM sleep and is active during wakefulness. In our modeling framework a group of cells can either be active or inactive. Often, it is reported that the firing rate of aminergic cells decreases during transitions from wake to NREM sleep, then becomes nearly quiescent during REM sleep (Lu et al. 2002). It is not clear if aminergic cell groups re-activate during transitions from NREM to REM sleep. Many reports of the activity of aminergic nuclei activity are based on sleep data in rodents, that unlike humans, rarely transition directly from REM to NREM (Diniz Behn et al. 2007; Tamakawa et al. 2006). If the activity of aminergic groups wanes during NREMS that occurs between wakefulness and REM sleep but remains low during REM-NREM cycling then our model would not differ significantly from data in this regard. In order to facilitate the possibility of analysis we have chosen a model that takes on only two distinct states. Future refinements of the model will reflect more nuanced features of the data. Our model leads to several insights and predictions concerning the human sleep/ wake cycle. For example, the model predicts that, during the normal sleep/wake cycle, both waking up and falling asleep are driven by the activity of VLPO. That is, the system falls asleep when activity of sleep-active cells rises above some threshold (instead of the activity of wake-active cells falling below some threshold). Moreover, waking up happens when the activity of sleep-active cells fall below a threshold. The model also helps explain the role of orexin in stabilizing the sleep/wake cycle. Saper and colleagues have suggested that an important role of the excitatory input from orexin neurons to the wake-promoting cells is to stabilize the flip–flop switch. The model and analysis presented here provides a mathematical interpretation for how orexin stabilizes the switch and why the loss of orexin leads to rapid switching between sleep and wakefulness. In our model, without excitatory inputs from ORX to AMIN, there is no longer a stable fixed point corresponding to the wake state. This leads to a short wake episode and, therefore, the sleep homeostat does not have time to build up sufficiently. This, in turn, leads to a short sleep episode. Our model also predicts how the short sleep/wake episodes should depend on circadian phase. Finally, as we have already discussed, in order to get a model based on the REM/ NREM switch hypothesis to generate features of the ultradian rhythms we need to posit an inhibitory connection from clVLPO to eVLPO and a drop in the activity of the eVLPO just before sleep onset. Our hope is that each of these predictions will motivate future experimental studies. Our model is intended to further the development of mathematical models for human sleep–wake cycles by incorporating structures with more direct biological interpretations than previously available. This work explores the explanatory power of the flip–flop model by testing its ability to account for the principal characteristics of

123

640

M. J. Rempe et al.

human sleep/wake regulation including consolidated wakefulness and REM/NREM cycling. With this focus on the neuronal interactions, we have used simple models for the circadian cycle and for homeostatic sleep pressure, adding these processes in a modular manner as suggested by Borbély and Achermann (1992). Moreover, we have neglected important factors likely to influence network dynamics. We have for instance treated the monoaminergic cell groups as a single population with coordinated activity that transitions between high and low states while spending little time at intermediate activity levels. The current model is therefore unable to address experimental data indicating that narcoleptics may experience states during which the tuberomammillary nucleus is active yet the locus coeruleus and dorsal raphe nucleus are not (John et al. 2004; Wu et al. 1999, 2004). Future refinement of the model will investigate the role of individual monoaminergic nuclei and will more flexibly incorporate firing patterns as well as a range of activity levels.

Appendix Circadian input To model the 24-h circadian input, we used the same model as in Daan et al. (1984): C(t) = [0.97 sin(ωt) + 0.22 sin(2ωt) + 0.07 sin(3ωt) + 0.03 sin(4ωt) + 0.01 sin(5ωt)] 2π and τ = 24. ω= τ Noise A separate input is added to each of the four cell groups (wake-active, sleep-active, REM-active and NREM-active) to represent noisy inputs from other brain regions. In each case the noise is modeled as follows: Ini = gnoise n i ,

(10)

where i represents AMIN, VLPO, REM, or NREM. The term n i is modeled with the following differential equation: n i! = wn (−n i + an N ),

(11)

where N is a normally distributed random variable with mean 0 and standard deviation 1. Both an and wn are constants. Each time step a new value is chosen for N from the normal distribution. If an N > n i at that time step n i will increase. Conversely, if an N < n i then n i will decrease. A time step of 3.6 seconds was used for all simulations. The values of gnoise , an , and wn can be found in Table 1.

123

Mathematical model of sleep Table 1 Parameter values used in the models for the wake-active, sleep-active, REM-on, and REM-off populations

641 Sleep/wake model Parameter

Value

Parameter

Value

εA

3

εR

0.3

εV

3

εN

0.3

γA

5.7

γR

6.2

γV

3.77

γN

6

τ1A

1h

τ1R

1h

τ2 A

2h

τ2R

2h

τ1V

1h

τ1N

0.5 h

τ1V

2h

τ2N

1.7 h

δA

0.01

δR

0.5

δV

0.01

δN

0.5

gvlpo

5

σ

11

gamin

2

gamin (A to REM)

2.5

gscn

1

gR E M

0.4

I0A

3.3

gN R E M

5

I0V

0.45

ae

2

ghom

5.5

be

1

αh

18.2 h −1

ce

1.82

βh

4.2 h−1

cv

hmax

1

kX V

−0.3

KU

1.55

kX

0.2

KL

−0.68

gevlpo

6.2

I0R

1.1

gnoise

All parameters without units listed are dimensionless

REM/NREM model

5

an

20

wn

0.01

*t

3.6 s

I0N

1

1.9

References Achermann P (2004) The two-process model of sleep regulation revisited. Aviat Space Environ Med 75(3):A37–A43 Achermann P, Borbély AA (1990) Simulation of human sleep-ultradian dynamics of electroencephalographic slow-wave activity. J Biol Rhythms 5(2):141–157 Achermann P, Borbély AA (1994) Simulation of daytime vigilance by the additive interaction of a homeostatic and a circadian process. Biol Cybern 71(2):115–121 Achermann P, Borbély AA (2003) Mathematical models of sleep regulation. Front Biosci 8:S683–S693 Bes FW, Jobert ML, Muller LC, Schulz H (1996) The diurnal distribution of sleep propensity: experimental data about the interaction of the propensities for slow-wave sleep and rem sleep. J Sleep Res 5:90–98 Boeve BF, Silber MH, Saper CB, Ferman TJ, Dickson DW, Parisi JE, Benarroch EE, Ahlskog JE, Smith GE, Caselli RC, Tippman-Peikert M, Olson EJ, Lin SC, Young T, Wszolek Z, Schenck CH, Mahowald MW, Castillo PR, Del Tredici K, Braak H (2007) Pathophysiology of REM sleep behaviour disorder and relevance to neurodegenerative disease. Brain 130:2770–2788

123

642

M. J. Rempe et al.

Boissard R, Gervasoni D, Schmidt MH, Barbagli B, Fort P, Luppi PH (2002) The rat ponto-medullary network responsible for paradoxical sleep onset and maintenance: a combined microinjection and functional neuroanatomical study. Eur J Neurosci 16(10):1959–1973 Borbély A (1982) A two process model of sleep regulation. Hum Neurobiol 1:195–204 Borbély A, Achermann P (1999) Sleep homeostasis and models of sleep regulation. J Biol Rhythms 14:557– 568 Borbély AA, Achermann P (1992) Concepts and models of sleep regulation, an overview. J Sleep Res 1:63–79 Campbell S, Czeisler C, Dijk D, Kronauer R, Brown E, Duffy J, Allan J, Shanahan T, Rimmer D, Ronda J, Mitchell J, Silva E, Emens J (2000) Is there an intrinsic period of the circadian clock? Science 288:1174–1175 Chemelli RM, Willie JT, Sinton CM, Elmquist JK, Scammell T, Lee C, Richardson JA, Williams SC, Xiong YM, Kisanuki Y, Fitch TE, Nakazato M, Hammer RE, Saper CB, Yanagisawa M (1999) Narcolepsy in orexin knockout mice: molecular genetics of sleep regulation. Cell 98(4):437–451 Chou T (2003) Regulation of wake-sleep timing: circadian rhythms and bistability of sleep–wake sates. Ph.D. thesis, Harvard University Chou TC, Bjorkum A, Gaus S, Lu J, Scammell T, Saper C (2002) Afferents to the ventrolateral preoptic nucleus. J Neurosci 22:977–990 Chou TC, Scammell TE, Gooley JJ, Gaus SE, Saper CB, Lu J (2003) Critical role of dorsomedial hypothalamic nucleus in a wide range of behavioral circadian rhythms. J Neurosci 23(33):10691–10702 Colwell CS, Michel S (2003) Sleep and circadian rhythms: do sleep centers talk back to the clock? Nat Neurosci 6(10):1005–1006 Czeisler C, Duffy J, Shanahan T, Brown E, Mitchell J, Rimmer D, Ronda J, Silva E, Allan J, Emens J, Dijk D, Kronauer R (1999) Stability, precision, and near-24-hour period of the human circadian pacemaker. Science 284:2177–2181 Daan S, Beersma DGM, Borbely AA (1984) Timing of human sleep—recovery process gated by a circadian pacemaker. Am J Phys 246(2):R161–R178 Dauvilliers Y, Amulf I, Mignot E (2007) Narcolepsy with cataplexy. Lancet 369(9560):499–511 Deurveilher S, Semba K (2006) Immediate early genes in sensory processing, cognitive performance and neurological disorders, chap. Mapping Sleep–wake control with the transcription factor c-Fos, pp 113–136. Springer, USA Dijk DJ, Czeisler CA (1994) Paradoxical timing of the circadian-rhythm of sleep propensity serves to consolidate sleep and wakefulness in humans. Neurosci Lett 166(1):63–68 Dijk DJ, Czeisler CA (1995) Contribution of the circadian pacemaker and the sleep homeostat to sleep propensity, sleep structure, electroencephalographic slow waves, and sleep spindle activity in humans. J Neurosci 15(5):3526–3538 Diniz Behn CG, Brown EN, Scammell TE, Kopell NJ (2007) Mathematical model of network dynamics governing mouse sleep–wake behavior. J Neurophysiol 97(6):3828–3840 Diniz Behn CG, Kopell N, Brown EN, Mochizuki T, Scammell TE (2008) Delayed orexin signaling consolidates wakefulness and sleep: physiology and modeling. J Neurophysiol 99(6):3090–3103 Economo Cvon (1930) Sleep as a problem of localization. J Nerv Ment Dis 71:249–259 Ferrillo F, Donadio S, De Carli F, Garbarino S (2007) A model-based approach to homeostatic and ultradian aspects of nocturnal sleep structure in narcolepsy. Sleep 30(2):157–165 FitzHugh R (1961) Impulses and physiological states in theoretical models of nerve membrane. Biophys J 1:445–466 Fuller PM, Saper CB, Lu J (2007) The pontine REM switch: past and present. J Phys Lond 584(3):735–741 Hill S, Tononi G (2005) Modeling sleep and wakefulness in the thelamocortical system. J Neurophysiol 93(3):1671–1698 Hodgkin A, Huxley A (1952) A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 117:500–544 Jewett ME, Kronauer RE (1999) Interactive mathematical models of subjective alertness and cognitive throughput in humans. J Biol Rhythm 14(6):588–597 John J, Wu M, Boehmer L, Siegel J (2004) Cataplexy-active neurons in the hypothalamus: implications for the role of histamine in sleep and waking behavior. Neuron 42(4):619–634 Kleitman N (1987) Sleep and wakefulness. University of Chicago Press, Chicago Krahn LE, Black JL, Silber MH (2001) Narcolepsy: new understanding of irresistible sleep. Mayo Clin Proc 76(2):185–194 Krueger JM, Obal F (1993) A neuronal group-theory of sleep function. J Sleep Res 2(3):186–186

123

Mathematical model of sleep

643

Lavie P (1986) Ultrashort sleep-waking schedule. III. ‘gates’ and ‘forbidden zones’ for sleep. Elect Clin Neuro 63:414–425 Lu J, Bjorkum AA, Xu M, Gaus SE, Shiromani PJ, Saper CB (2002) Selective activation of the extended ventrolateral preoptic nucleus during rapid eye movement sleep. J Neurosci 22(11):4568–4576 Lu J, Sherman D, Devor M, Saper CB (2006) A putative flip–flop switch for control of rem sleep. Nature 441(7093):589–594 Massaquoi S, McCarley R (1992) Extension of the limit cycle reciprocal interaction model of REM cycle control. an integrated sleep control model. J Sleep Res 1:138–143 McCarley R, Massaquoi S (1992) Neurobiological structure of the revised limit cycle reciprocal interaction model of REM cycle control. J Sleep Res 1:132–137 McCarley RW (2007) Neurobiology of REM and NREM sleep. Sleep Med 8(4):302–330 McCarley RW, Hobson JA (1975) Neuronal excitability modulation over sleep cycle—structural and mathematical-model. Science 189(4196):58–60 McCarley RW, Massaquoi SG (1986) A limit-cycle mathematical-model of the rem-sleep oscillator system. Am J Physiol 251(6):R1011–R1029 McGinty D, Szymusiak R (2003) Hypothalamic regulation of sleep and arousal. Front Biosci 8:S1074– S1083 Merica H, Fortune RD (2004) State transitions between wake and sleep, and within the ultradian cycle, with focus on the link to neuronal activity. Sleep Med Rev 8(6):473–485 Mochizuki T, Crocker A, McCormack S, Yanagisawa M, Sakurai T, Scammell TE (2004) Behavioral state instability in orexin knock-out mice. J Neurosci 24(28):6291–6300 Morris C, Lecar H (1981) Voltage oscillations in the barnacle giant muscle fiber. Biophys J 35:193–213 Nagumo J, Arimoto S, Yoshizawa S (1962) An active pulse transmission line simulating nerve axon. Proc IRE 50:2061–2070 Nobili L, Beelke M, Besset A, Billiard M, Ferrillo F (2001) Nocturnal sleep features in narcolepsy: a model-based approach. Rev Neurol 157(11):S82–S86 Phillips A, Robinson P (2008) Sleep deprivation in a quantitative physiologically based model of the ascending arousal system. J Theor Biol 255(4):413–423 Phillips AJK, Robinson PA (2007) A quantitative model of sleep–wake dynamics based on the physiology of the brainstem ascending arousal system. J Biol Rhythms 22(2):167–179 Plazzi G, Serra L, Ferri R (2008) Nocturnal aspects of narcolepsy with cataplexy. Sleep Med Rev 12:109– 128 Porkka-Heiskanen T, Strecker RE, Thakkar M, Bjorkum AA, Greene RW, McCarley RW (1997) Adenosine: a mediator of the sleep-inducing effects of prolonged wakefulness. Science 276(5316):1265– 1268 Rubin JE, Shevtsova NA, Ermentrout GB, Smith JC, Rybak IA (2009) Multiple rhythmic states in a model of the respiratory central pattern generator. J Neurophysiol 101(4):2146–2165 Sakurai T (2007) The neural circuit of orexin (hypocretin): maintaining sleep and wakefulness. Nat Rev Neurosci 8(3):171–181 Saper CB, Chou TC, Scammell TE (2001) The sleep switch: hypothalamic control of sleep and wakefulness. Trends Neurosci 24(12):726–731 Saper CB, Scammell TE, Lu J (2005) Hypothalamic regulation of sleep and circadian rhythms. Nature 437(7063):1257–1263 Sasaki Y, Fukuda K, Takeuchi T, Inugami M, Miyasita A (2000) Sleep onset REM period appearance rate is affected by REM propensity in circadian rhythm in normal nocturnal sleep. Clin Neurophysiol 111(3):428–433 Siegel JM (2006) The stuff dreams are made of: anatomical substrates of rem sleep. Nat Neurosci 9(6):721– 722 Siegel JM (2008) Do all animals sleep? Trends Neurosci 31(4):208–213 Strogatz S (1986) The mathematical structure of the human sleep–wake cycle. Springer, Heidelberg Strogatz SH (1987) Human sleep and circadian-rhythms—a simple-model based on 2 coupled oscillators. J Math Biol 25(3):327–347 Swick TJ (2005) The neurology of sleep. Neurol Clin 23(4):967–989 Szymusiak R, Gvilia I, McGinty D (2007) Hypothalamic control of sleep. Sleep Med 8(4):291–301 Takeuchi T, Ogilvie R, Murphy T, Ferrelli A (2003) EEG activities during elicited sleep onset REM and NREM periods reflect different mechanisms of dream generation. Clin Neurophysiol 114:210–220 Tamakawa Y, Karashima A, Koyama Y, Katayama N, Nakao M (2006) A quartet neural system model orchestrating sleep and wakefulness mechanisms. J Neurophysiol 95(4):2055–2069

123

644

M. J. Rempe et al.

Wright K, Hughes R, Kronauer R, Dijk D, Czeisler C (2001) Intrinsic near-24-h pacemaker period determines limits of circadian entrainment to a weak synchronizer in humans. Proc Natl Acad Sci 98:14027– 14032 Wu M, Gulyani S, Yau E, Mignot E, Phan B, Siegel J (1999) Locus coeruleus neurons: cessation of activity during cataplexy. Neuroscience 91(4):1389–1399 Wu M, John J, Boehmer L, Yau D, Nguyen G, Siegel J (2004) Activity of dorsal raphe cells across the sleep-waking cycle and during cataplexy in narcoleptic dogs. J Physiol 554(1):202–215 Yin W (2007) A mathematical model of the sleep–wake cycle. Master’s thesis, Georgia Institute of Technology

123