Sleep deprivation in a quantitative physiologically based model of the ascending arousal system

ARTICLE IN PRESS Journal of Theoretical Biology 255 (2008) 413–423

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Sleep deprivation in a quantitative physiologically based model of the ascending arousal system A.J.K. Phillips a,b,, P.A. Robinson a,b,c a b c

School of Physics, University of Sydney, New South Wales 2006, Australia Brain Dynamics Center, Westmead Millennium Institute, Westmead Hospital and Western Clinical School, University of Sydney, Westmead, New South Wales 2145, Australia Faculty of Medicine, University of Sydney, New South Wales 2006, Australia

a r t i c l e in fo

abstract

Article history: Received 17 December 2007 Received in revised form 18 August 2008 Accepted 20 August 2008 Available online 29 August 2008

A physiologically based quantitative model of the human ascending arousal system is used to study sleep deprivation after being calibrated on a small set of experimentally based criteria. The model includes the sleep–wake switch of mutual inhibition between nuclei which use monoaminergic neuromodulators, and the ventrolateral preoptic area. The system is driven by the circadian rhythm and sleep homeostasis. We use a small number of experimentally derived criteria to calibrate the model for sleep deprivation, then investigate model predictions for other experiments, demonstrating the scope of application. Calibration gives an improved parameter set, in which the form of the homeostatic drive is better constrained, and its weighting relative to the circadian drive is increased. Within the newly constrained parameter ranges, the model predicts repayment of sleep debt consistent with experiment in both quantity and distribution, asymptoting to a maximum repayment for very long deprivations. Recovery is found to depend on circadian phase, and the model predicts that it is most efficient to recover during normal sleeping phases of the circadian cycle, in terms of the amount of recovery sleep required. The form of the homeostatic drive suggests that periods of wake during recovery from sleep deprivation are phases of relative recovery, in the sense that the homeostatic drive continues to converge toward baseline levels. This undermines the concept of sleep debt, and is in agreement with experimentally restricted recovery protocols. Finally, we compare our model to the two-process model, and demonstrate the power of physiologically based modeling by correctly predicting sleep latency times following deprivation from experimental data. & 2008 Elsevier Ltd. All rights reserved.

Keywords: Mathematical model Brainstem Computational model Sleep latency Sleep debt

1. Introduction Sleep deprivation has been studied widely with the intention of discovering the function of sleep, and the nature of the recovery process. Such studies also provide useful insights into the mechanisms governing the human sleep cycle, and can thus assist in the building and calibration of models. Better understanding recovery from sleep deprivation has a myriad of applications, including medicine, aviation, and general public health. In this instance, we are interested in improving an existing model, based on dynamics known from experiment. We use a small number of experimentally derived criteria to calibrate the model, and the refined model is then used to make more general predictions about other experiments, and the recovery process as a whole.

 Corresponding author at: School of Physics, University of Sydney, New South Wales 2006, Australia. Tel.: +61 9036 7960; fax: +61 93517726. E-mail address: [email protected] (A.J.K. Phillips).

0022-5193/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2008.08.022

The overall arousal state of the brain is controlled by a series of brainstem nuclei, collectively termed the ascending arousal system (AAS), which diffusely projects neuromodulators to the cerebrum (Kandel et al., 1991; Jones, 2000; Pace-Schott and Hobson, 2002). The AAS is regulated by homeostatic and circadian sleep drives, which are integrated in the hypothalamus, projecting to the AAS from there (Pace-Schott and Hobson, 2002). Based on common effects and temporal activities, AAS nuclei can be broadly classified into two groups: (i) monoaminergic (MA) and (ii) acetylcholine-related (ACh). The MA group includes nuclei that use MA neurotransmitters (Aston-Jones et al., 1991; Jones, 2000; Webster, 2002), while the ACh group includes nuclei which express acetylcholine as well as nuclei excited cholinergically to share a common activity pattern (Thakkar et al., 1998; Jones, 2000; Saper et al., 2001; Pace-Schott and Hobson, 2002), as shown in Fig. 1(a). Circadian and homeostatic drives affect the AAS via GABAergic inhibitory projections from the ventrolateral preoptic area (VLPO) of the hypothalamus (Sherin et al., 1996; Pace-Schott and Hobson, 2002; Saper et al., 2005), as shown in Fig. 1(a). The 24-h circadian

ARTICLE IN PRESS 414

A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

(a)

(b) arousal state

LC DR TMN

BRF

LDT/ PPT

MA

ACh

VLPO VLPO/ eVLPO

In this process, we demonstrate the wide scope of application of the model rather than fitting the model to an individual study, starting with repayment of sleep debt in Section 4. We consider circadian phase dependence of recovery, and optimal recovery plans, in Section 5. We then compare the model to the phenomenological two-process model of sleep regulation (Borbe´ly and Achermann, 2000) in Section 6, and demonstrate the advantages of physiologically based modeling by simulating a set of experimental sleep latency data. Finally, the results are discussed in Section 7.

D

Fig. 1. Parts (a) and (b) show schematics of the actual AAS populations, and our sleep model, respectively. Excitatory inputs are represented by solid arrow heads, and inhibitory by open ones. In each case, the top left box is the MA group, and the top right box is the ACh group. In (a), the MA group consists of the locus ceruleus (LC), dorsal raphe (DR) and tuberomamillary nucleus (TMN); and the ACh group consists of cholinergic laterodorsal and pedunculopontine tegmentum (LDT/PPT), and glutamergic brainstem reticular formation (BRF). The ventrolateral preoptic area (VLPO/eVLPO) GABAergically inhibits other AAS nuclei. In (b) the drive D is shown, which consists of circadian (C) and homeostatic (H) components. In our model the thick-lined interactions in (b) are used.

2. Sleep model We begin by briefly reviewing the structure of the Phillips and Robinson (2007) model of the AAS in Section 2.1. The circadian and homeostatic drives to the model are then discussed in greater detail in Section 2.2 due to their high relevance in modeling sleep deprivation. 2.1. Model structure

signal keeps the sleep cycle entrained with the day–night light cycle (Pace-Schott and Hobson, 2002), and originates in the suprachiasmatic nucleus (SCN), being projected to the VLPO via the dorsomedial nucleus of the hypothalamus (Saper et al., 2005; Gooley et al., 2006). The homeostatic drive is a somnogenic drive which increases with time spent awake, one effect of which is presumably to prevent excessive sleep deprivation. It is believed to result from net build up of somnogens during waking, including extracellular adenosine in the basal forebrain (Jones, 2000; McKenna et al., 2003), resulting in disinhibition of the VLPO (Pace-Schott and Hobson, 2002; Saper et al., 2005). During sleep, metabolic rates are low and adenosine is cleared faster than it is produced (Dı´az-Mun´oz et al., 1999; McKenna et al., 2003). Mutual inhibition between the MA group and the VLPO gives rise to the sleep–wake ‘flip-flop’, with each group disinhibiting its own activity. Thus, during wake the MA group is active and the VLPO suppressed, with the converse true during sleep. The model used in this paper represents the above physiology, focusing on the ‘flip-flop’, and is briefly reviewed in Section 2. In previous work, our model was introduced in the context of the normal sleep–wake cycle, allowing constraint of most model parameters (Phillips and Robinson, 2007). However, the previous version of the model was not calibrated to accurately capture recovery from sleep deprivation. Now by applying a small set of constraints on how the model responds to sleep deprivation, we will constrain the circadian and homeostatic drive parameters accurately, improving the utility of the model, and allowing predictions to be made about protocols which lie beyond those used for constraint. Biomathematical models have studied sleep deprivation in some depth, and have used it as a calibration method (Daan et al., 1984; Achermann et al., 1993; Mallis et al., 2004). However, because these models are phenomenological in derivation, it is difficult to make direct inferences regarding the underlying physiology. A few recent models which have been based more closely on physiology (Tamakawa et al., 2006; Yin, 2007; Diniz Behn et al., 2007) provide potential for such investigations; but they contain significantly more parameters than our model, making the process of constraining parameters more difficult. In Section 2 we briefly review the sleep model used, and in Section 3 we use a handful of experimentally derived criteria to calibrate the model. We then simulate several other protocols to investigate experiments that were not part of the calibration set.

We use a continuum neuronal population model, in which we ascribe average properties to large populations of neurons and their interactions (Nunez, 1995; Robinson et al., 1997). To model the AAS in tractable form, we divided the nuclei into MA and ACh groups, based on common activity patterns and neuromodulatory effects (Phillips and Robinson, 2007). The third population modeled was the VLPO, where circadian (C) and homeostatic (H) drives are integrated and projected to brainstem nuclei. The model incorporates the sleep–wake flip-flop between the VLPO and MA populations. The ACh group is held at a constant intermediate activity level, smoothing out the 90 min ultradian oscillation between REM and NREM, while retaining the main features of the sleep–wake dynamics. Non-ultradian dynamics on a timescale shorter than 90 min can still be considered, such as sleep–wake transitions, since the MA group is chiefly responsible for such transient events (Aston-Jones et al., 1991; Jones, 2000). A schematic of the model and its connectivities is seen in Fig. 1(b), including only connections for which there is strong evidence. Each population j ¼ v; m, where v is VLPO and m is MA, has a mean cell body potential V j ðtÞ relative to resting and a mean firing rate Q j ðtÞ. Due to the small sizes of the nuclei, we assume spatial homogeneity of each population and neglect propagation delays. The relationship of Q j to V j is approximated by a sigmoid (Freeman, 1975), with Q j ¼ SðV j Þ ¼

Q max , 1 þ exp½ðV j  yÞ=s0 

(1)

y is the mean firing where Q max is the maximum possibleprate; ffiffiffi threshold relative to resting, and s0 p= 3 is its standard deviation (Robinson et al., 1997, 2004), determining the sigmoid width. Neuronal dynamics are represented by

tv dV v =dt þ V v ¼ nvm Q m þ D,

(2)

tm dV m =dt þ V m ¼ A þ nmv Q v ,

(3)

where A represents constant input from the acetylcholine and orexin groups, the njk weight the input from population k to j, tj is the decay time for the neuromodulator expressed by group j, and the total drive is D ¼ nvc C þ nvh H.

(4)

The nvc and nvh weight the strengths of the circadian and homeostatic drives to the VLPO, respectively. At present, circadian and homeostatic drives to the MA group are neglected, since their

ARTICLE IN PRESS A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

exact nature is still being verified experimentally (Mistlberger, 2005). During normal functioning the model evolves on two timescales, one being the timescale of D (about 24 h), the other being the rapid transitions between wake and sleep states (about 10 min). Considering D to be slowly varying, the model has up to two stable steady state solutions, corresponding to wake and sleep. During wake, V m is high, V v is low, and H is growing; whereas during sleep, V m is low, V v is high, and H is falling (Phillips and Robinson, 2007). 2.2. Model drives As a first approximation the circadian drive is considered well entrained with the daily fluctuation in light intensity, so is modeled by a sinusoidal function of time, CðtÞ ¼ c0 þ cos½o ðt  aÞ,

(5)

1

where o ¼ ð2p=24Þ h is the angular frequency corresponding to the 24 h period of Earth’s rotation, and a is constant, defining the initial phase of the circadian oscillation. The phase of the sleep–wake cycle relative to the circadian drive depends on the model parameters, particularly the drive parameters nvc , nvh , and c0 . With the new parameter set derived in Section 3, t ¼ a corresponds to approximately 4 pm. Note that because C is scaled by nvc in its interaction with the VLPO, we can use an oscillation amplitude of 1, without loss of generality. In previous work we assumed c0 ’ 1 on the basis of low minimum SCN firing rates corresponding to C ’ 0. Here we allow c0 a1, noting that c0 also encapsulates any other net input to the VLPO from nuclei that are not explicitly modeled. In general, a dynamic oscillator model of the SCN could be used, but in this paper we study sleep deprivation in isolation from other phenomena. We model homeostatic sleep drive by the somnogen level H, which we consider to be the concentration of extracellular adenosine in the basal forebrain. This is discussed in more detail in Section 7, and further details can be found in Phillips and Robinson (2007). The exact physiological pathway for the homeostatic drive is not yet clear, but its dynamical characteristics for humans have been sufficiently well explored for us to model it, with the form motivated by the two-process model (Borbe´ly and Achermann, 2000). In our model, adenosine levels increase in wake due to higher metabolic activity, and decrease in sleep when clearance exceeds production. Presently, we use the simplest realistic forms for production and clearance rates. Clearance rate is assumed to be proportional to concentration, and production is approximated as a linear function of Q m alone, since MA activity is well correlated with arousal (Aston-Jones et al., 1991). As the vast majority of the sleep–wake cycle is spent in either sleep or wake, for each of which Q m is well constrained, only two points are spanned, justifying a linear function of Q m as an adequate fit. In cases where significant time is spent transitioning between wake and sleep, such as sleep apnea, this may not be a suitable assumption, and the form of the production rate will need to be investigated further. Assuming that Q m ¼ 0 results in negligible production yields

wdH=dt þ H ¼ mQ m ,

(6)

where w is the characteristic time for somnogen clearance and m is constant. The value of H grows during wake to a maximum value, Hhi , at the time of sleep onset, and decreases during sleep to a minimum value, Hlo , at the time of waking. To model sleep deprivation in this paper, we hold dV m =dt ¼ dV v =dt ¼ 0 with V m and V v at their waking values, with H ¼ Hlo initially. Deprivation can begin at any time of day we choose to model. During the deprivation, H grows. To model recovery, these constraints on V v and V m are then removed, being replaced by Eqs. (2) and (3), and

415

the excess H causes recovery sleep. Solving piecewise in this fashion allows us to graph the entire simulation of the experimental deprivation protocol.

3. Model calibration We now use a small number of experimentally derived criteria to calibrate the model, so as to yield a region in parameter space within which the model behaves realistically in response to sleep deprivation and recovery. These criteria are defined in Section 3.1. Homeostatic drive parameters are approximated based on adenosine concentration data in Section 3.2, and the form of the total drive is constrained in Section 3.3. Using a nominal drive amplitude and mean, the calibration criteria are applied in Section 3.4, yielding a better set of constraints on the drive parameters. Ultradian oscillations between REM and NREM are not yet considered, so for model comparison we only use total sleep deprivation experiments. 3.1. Calibration criteria With its original parameter set the model replicated the normal human sleep–wake cycle, with H varying over ranges consistent with experimentally measured concentrations of adenosine. Parameters governing the neuronal dynamics and connections were well constrained based on both physiology and dynamics, but the drive parameters were not so well constrained, since the model’s normal sleep–wake dynamics are not very sensitive to the relative sizes of the circadian and homeostatic drives. The form of the total drive D was estimated based on normal sleeping time and time spent in a region of bistability of wake and sleep states (Phillips and Robinson, 2007), and it is better constrained in Section 3.3. The homeostatic drive’s form and size relative to the circadian drive were not well constrained, since the normal sleep cycle is relatively insensitive to these factors. The parameters controlling somnogen production and clearance—m and w, respectively—were estimated by assuming that somnogen concentration ranges from 8 to 16 nM, using rat data, since similar cellular mechanisms and hence similar concentration ratios, are expected across mammals. The concentration of H may need to be rescaled to fit human data when it becomes available, but only the ratio of maximum to minimum concentrations is important; a point which is addressed in Section 3.2. Meanwhile, the relative weightings of the circadian and homeostatic drives, controlled by nvc and nvh , were estimated by assuming homeostatic pressure is sufficient to induce sleep after 60 h of waking, without additional stimulus. Bounds on these parameters were provided purely by the constraints on the total drive D. To derive improved estimates of the drive parameters we consider the model’s response to sleep deprivation, since it depends sensitively on these parameters. The model response should be consistent with dynamics of the human sleep cycle known from experiment. Specifically, we require that: (i) The concentration of somnogens H remains physiologically plausible. (ii) The amount of recovery sleep is consistent with experiment. (iii) The model response depends on the length of deprivation up to about 100 h, since significant changes in EEG continue to be observed beyond this point, including loss of normal alpha activity on eye closure after about 115 h without sleep (Bonnet, 2000). (iv) Without any increase in stimulus, the homeostatic drive is sufficient to induce sleep at a circadian peak (i.e., when the SCN is most active), after 40–80 h of continuous wake, since the ability to concentrate decreases markedly after about 60 h of waking (Bonnet, 2000; Williams et al., 1959).

ARTICLE IN PRESS A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

We first restrict H to physiologically plausible values, to satisfy constraint (i). In earlier work we found Q m to vary little within wake and sleep, with rapid transitions between the states. We can thus assume approximately constant Q m in wake and sleep, labeling these values Q w and Q s , respectively. This results in exponential dynamics during each state, HðtÞ ¼ mQ m þ ½Hð0Þ  mQ m et=w ,

(7)

where Hð0Þ is the initial value of H. During wake, mQ m 4H so HðtÞ exponentially saturates, whereas during sleep mQ m oH so HðtÞ exponentially decays. If we label the maximum and minimum values of H as Hhi and Hlo , respectively, (7) gives Hhi ¼ mQ w þ ðHlo  mQ w Þetw =w ,

(8)

Hlo ¼ mQ s þ ðHhi  mQ s Þets =w ,

(9)

where t w and t s are the durations of wake and sleep in a normal day, respectively, so that T day ¼ t w þ t s ¼ 24 h. Using the above approximation of constant Q m in each phase, ¯ w and H ¯ s , the mean levels of H in wake and we can calculate H sleep, respectively: ¯ w ¼ mQ w þ wðHlo  mQ w Þð1  etw =w Þ=t w , H

(10)

¯ s ¼ mQ s þ wðHhi  mQ s Þð1  ets =w Þ=t s . H

(11)

¯ ¼ ðt w H ¯ w þ ts H ¯ s Þ=T day yields Then, H ¯ ¼ ½t s mQ s þ wðHhi  mQ s Þð1  ets =w Þ þ t w mQ w H þ wðHlo  mQ w Þð1  etw =w Þ=T day .

(12)

Parameter

Physiological range

Constrained range

¯ D DA nvc

– 40 40 40

1.8–1.7 1.4–4.9 2.7–3.2 1.0–2.9

10–72 40 –

34–72 2.1–7.2 4.1–12.7

o150 10–20 2–6 40 40 40 0.1–600 0.1–600

o150 10–20 2–6 0–1.5 1.6–2.3 1.6–2.3 1.0–20 1.0–20

nvh w m c0 Q max

y

s0 A nvm nmv

tm tv

(13)

Based on model behavior we allow the ranges 0:1 s1 oQ s o0:5 s1 and 3:5 s1 oQ w o5 s1 , and for a healthy sleep duration we use 6 hot s o10 h. Based on adenosine concentrations in the basal ¯ forebrain, we estimate the range 8 nMoHo15 nM (McKenna et al., 2003). These yield 2:1 nM somo7:2 nM s, which overlaps the previous range 3:4 nM somo4:0 nM s. Based on model firing rates with nominal parameters we take Q s ¼ 0:3 s1 and Q w ¼ 4:2 s1 , and as nominal values from near the middle of the ¯ ¼ 12 nM, yielding m ¼ 4:4 nM s. We range we take t s ¼ 8 h and H note that the relative shape of H (i.e., the ratio of maximum to minimum concentrations across the daily cycle) is more important than the precise values of the concentration—which are hard to obtain for humans—since renormalizing H and m by any multiplicative factor a, and nvh by 1=a, yields identical dynamics. 3.3. Weighting of drives The relative weighting of the homeostatic and circadian drives will determine how sensitive the model is to sleep deprivation. Using the old model parameters in Table 1 (Phillips and Robinson, 2007) Fig. 2 shows a recovery from a 200 h sleep deprivation. It results in only 1.5 h of recovery sleep above baseline sleep duration, with almost all of this on the first night. By the second night of recovery, sleep duration is increased by less than 1%

New nominal

Unit

0.77 0.42 3.2 0.19

0.6 3.4 2.9 1.0

mV mV mV

11 3.6 1.0

45 4.4 4.5

100 10 3 1.0 1.9 1.9 10 10

100 10 3 1.3 2.1 1.8 10 10

mV nM1 h nM s 1 s1 mV mV mV mV s mV s s s

5 0 -5 -10 0 -5 -10 15

Substituting (8) and (9) into (12), we obtain ¯ day =ðQ w t w þ Q s t s Þ. m ¼ HT

Old nominal

Use of dynamical and physiological constraints in tandem yields the constrained range. Nominal values are chosen within the constraints. Parameters above the line are those that have been better constrained in this paper. Others are discussed in Phillips and Robinson (2007). The methods of constraint are discussed in Section 3.1, with parameters subject to improved constraints above the line. Note that the ¯ and DA , with constraints on nvc , nvh , w, and c0 are subject to the values chosen for D those above being for the new nominal values. In previous work (Phillips and Robinson, 2007), we used C ¼ ½1 þ cosðotÞ=2, so nvc and c0 were defined slightly differently; the old values have been re-expressed in the current form in this table.

Vm (mV)

3.2. Somnogen concentrations

Table 1 Parameter constraints and nominal values

Vv (mV)

Constraint (i) is applied in Section 3.2, and the form of the overall drive is constrained in Section 3.3. The response of the model to sleep deprivation is then studied in Section 3.4, and constraints (ii), (iii), and (iv) define the physiologically realistic regions of parameter space more accurately than in previous work, giving an improved set of nominal parameters.

H (nM)

416

10 5

0

24

48

72

t (h) Fig. 2. Time series for (a) V m , (b) V v , and (c) H, using old nominal parameters. Recovery from 200 h of wake is plotted (solid), along with the normal cycle (dashed). V m is high during wake and low during sleep, with the converse true for V v . Periods of sleep are indicated by shading for the recovery. The approximate clock time at t ¼ 0 is 11 pm.

above normal levels, suggesting that somnogen clearance is too fast with the old value of w ¼ 10 h. The distribution of recovery sleep is almost purely dictated by circadian phase, with sleep occurring only during the night. This is unrealistic for modeling long deprivation experiments, where sufficient tiredness causes recovery sleep to occur at any circadian phase. These results suggest that the model dynamics are circadian dominated for the old parameters, with the homeostatic drive making too small a contribution. We thus investigate varying the relative sizes of the drive couplings nvc and nvh , under the constraint that the normal sleep cycle is unaffected. This leads us to a reduction of the parameter space. ¯ of the total drive D, To estimate the amplitude DA and mean D we approximate H as a sinusoid, with H leading C by a quarter

ARTICLE IN PRESS A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

period, as in previous work (Phillips and Robinson, 2007), because C peaks near the middle of wake, while H peaks at the onset of sleep. We stress that this approximation is only used here for analytic purposes; all our numerical modeling uses (5) and (6) as ¯ and the drives. The mean and amplitude of H are H,

DH ¼ ðHhi  Hlo Þ=2 ¼ mðQ w  Q s Þð1  etw =w  ets =w þ eT day =w Þ=½2ð1  eT day =w Þ, (14) respectively. This yields DA ¼ ðn2vc þ n2vh DH2 Þ1=2 ,

(15)

¯ ¼ nvc c0 þ nvh H ¯ ¼ nvc c0 þ mðQ w t w þ Q s t s Þ=T day . D

(16)

Rearranging (15) gives nvh as a function of nvc and w,

nvh ¼ ðD2A  n2vc Þ1=2 =DH,

(17)

which yields the constraint jnvc joDA . Substituting (17) in (16) gives c0 as a function of nvc and w: ¯ 2A  n2vc Þ1=2 ¯ DH  HðD D . nvc DH

c0 ¼

(18)

¯ can be constrained by The overall drive parameters DA and D requiring sleep duration to lie in the range 6–10 h and hysteresis loop crossing time to lie in the range 1–3 h, based on the adaptability of the human sleep–wake cycle (Phillips and Robinson, 2007). This ensures that changes we make to the model parameters in improving the sleep deprivation response do not adversely affect normal model behavior. These constrained regions in parameter space are shown in Fig. 3, with the overlapping region satisfying both constraints. As nominal values from near the middle of the range we choose DA ¼ 3:4 mV and ¯ ¼ 0:6 mV. We note that the old nominal value of DA ¼ 0:42 mV D lies outside the more tightly constrained range, because slightly different values were used for nvm, nmv , and A to adjust the sleep duration to 8.5 h, although the constrained ranges of these variables remain the same. Using Eqs. (17) and (18) to define nvh and c0 , respectively, in terms of w and nvc , reduces the dimension of the parameter space. The response of the system to sleep

5

D A (mV)

4

3

2

1

0 −2

−1.5

−1

−0.5

0

¯ (mV) D

0.5

1

1.5

¯ with darker Fig. 3. Constraints on the parameter space spanned by DA and D, regions satisfying more constraints. The region between the solid lines corresponds to a sleep duration in the range 6–10 h. The region between the dotted lines corresponds to a hysteresis loop crossing time of 1–3 h. Nominal values (heavy dot) are chosen to satisfy all constraints.

417

deprivation can then be examined in the two dimensional space ¯ will yield slightly spanned by w and nvc . Each selection of DA and D different constraints on w and nvc , so for tractability we settle on ¯ We then find constraints on the the nominal values for DA and D. model drive parameters w, nvc , nvh , and c0 , and explore the sensitivity of the model to these parameters afterwards. 3.4. Calibration of model response to sleep deprivation We now examine how the model responds to sleep deprivation ¯ The other as we vary w and nvc , given nominal values for DA and D. drive parameters, nvh and c0 , are varied as a function of these parameters, using (17) and (18). By applying constraints (ii)–(iv) of Section 3.1, we find the regions of parameter space that correspond to realistic behavior. Since all deprivations—including over 200 h—typically require only 1–3 nights of recovery (Bonnet, 2000), we have an approximate physiological range of 10 howo72 h, which we next improve based on dynamics in a restricted set of experiments. For constraint (ii), we make a reasonable estimate to account for the large variability between studies, that in response to 40, 80, and 120 h of continuous waking, the amount of recovery sleep above baseline is 10–30% of the sleep missed. For constraint (iii), we require that the model’s response to 120 h of deprivation be clearly distinguishable from that after 80 h of deprivation. Quantitatively, we require that the amount of recovery sleep above baseline levels be at least 10% different. Finally, constraint (iv) requires that somewhere between 40 and 80 h of waking, the homeostatic drive becomes strong enough to induce sleep even at a circadian peak. The ability to concentrate decreases markedly on this timescale, as signified by an inability to maintain alpha rhythms in EEG (Bonnet, 2000), and a significant increase in maximum reaction times (Williams et al., 1959), and we take this decrease in cognitive performance to be a marker of a greater need for the maintenance of wake. In each case, recovery begins at a ¼ 16 h, corresponding approximately to the normal bedtime, although this will itself depend on the parameters being constrained. These constraints are represented in Fig. 4 by overlapping regions. The region of consistency between all constraints is termed the constrained range, and near its center we choose nominal values of nvc ¼ 2:9 mV and w ¼ 45 h. From there we choose a self-consistent parameter set, c0 ¼ 4:5 and nvh ¼ 1:0 mV nM1 , where we note that the bounds on these two parameters are only approximate, given that relations (17) and (18) are themselves approximate due to the non-sinusoidal form of H. Defining the waking time as 7:00 am, which corresponds to a ’ 8 h, the model bedtime is 10:27 pm, corresponding to a ’ 17 h, with nominal parameters. Using the nominal parameters from Table 1 the behavior of the model in response to 80 h of continuous wake is shown in Fig. 5. During deprivation H continually increases, although there are times when the rate of increase in H is lower than for baseline, resulting in decrease in DH. These correspond to normal wake times, and this point is addressed in more detail in Section 4.2. Initially, despite the circadian drive being maximal, the model stays in the sleep state, with V m low and V v high. The level of H is initially high, and it takes a couple of nights to return to normal. As seen in Fig. 5, sleep duration is greater during recovery, and the sleep is deeper in the sense that there is higher V v and lower V m. It seems reasonable to assume that greater suppression of the MA group will result in greater EEG slow wave activity, consistent with experimental findings (Bonnet, 2000), although this needs further exploration via experimental measures and modeling of the effects of the AAS on the thalamocortical system. The drive parameter values in the previous paper (Phillips and Robinson, 2007) gave too rapid a saturation of H, as well as too much

ARTICLE IN PRESS 418

A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

100

that are independent of the criteria used for constraint. Specifically, it explains saturation of sleep debt for long deprivations, and how recovery is distributed. In Section 4.1, we consider saturation effects, explore how recovery sleep is distributed across the following nights, and examine sensitivity of the dynamics to the parameter w. In Section 4.2, the role of wake in recovery of somnogens to normal levels is investigated.

(ii)

χ (h)

80 (iv)

60

4.1. Repayment of sleep debt

40 20

χ

−3.2

−3

−2.8

ν

vc

(h)

−2.6

−2.4

Δ H (nM) H (nM) Vv (mV) Vm (mV)

Fig. 4. Constraints on the parameter space spanned by nvc and w, with darker regions satisfying more constraints. By definition of the drive amplitude, jnvc joDA ¼ 3:4 mV, while physiological constraints force 10 howo72 h, corresponding to the region between the dashed lines. Satisfying constraint (ii) of Section 3.1 for 40, 80, and 120 h of continuous wake, corresponds to lying between the three lower and three upper solid lines. Constraint (iii) corresponds to the region above the dotted line, and constraint (iv) to the region between the thick solid lines. Nominal values (heavy dot) are chosen to satisfy all constraints.

0 -15 -30 5 -5 -15 20 15 5 0

0

24

48

72

96 t (h)

120

144

168

192

Fig. 5. Time series of (a) V m , (b) V v , (c) H, and (d) the amount by which H exceeds baseline levels in the deprivation protocol. We simulate 80 h of continuous waking by holding V m and V v at their waking values, between a normal waking time at t ¼ 40 and 120 h, followed by recovery (solid), with the normal cycle overlaid (dashed). Recovery begins at circadian phase a ¼ 0 h, corresponding to maximum circadian drive (approximately 2 pm). The new nominal parameters in Table 1 are used. Following recovery, the model returns to a normal sleep–wake cycle. Periods of sleep for the deprivation and recovery case are illustrated by shading.

weighting toward the circadian drive. With these improved parameter constraints and nominal values, the model still behaves consistently with previous work in the normal sleep cycle— as seen following recovery in Fig. 5, but now also responds realistically to sleep deprivation. Having completed calibration of the model on a restricted set of experiments, we now turn to predictions of other experiments, and comparison with published data.

As deprivation increases the amount of sleep required for recovery of somnogens to normal levels also increases (PorkkaHeiskanen et al., 2000). However, experiments show that this need for sleep saturates, since if left undisturbed, young adults usually initially sleep for only 12–15 h even after the most extreme deprivations (Johnson et al., 1965). This is consistent with our model of the homeostatic drive, which saturates exponentially to a maximal value of H ¼ mQ w with time awake, from (7), similar to the dynamics of process H in two-process models (Daan et al., 1984). Our model actually gives a maximum of 13 h sleep on the first night of recovery, as shown in Fig. 6. This is in good agreement with experiment (Johnson et al., 1965), as is the fact that duration of sleep on later nights depends less on the length of deprivation, as sleep need recovers to baseline levels with time. Another point of interest is how recovery sleep is distributed following deprivation. The model response to 60, 120, and 180 h of deprivation is shown in Fig. 7, beginning recovery at a ¼ 17 h (normal bedtime). In each case, the model predicts that recovery is effectively complete after the third night, in agreement with experiment (Bonnet, 2000); and around 60% of recovery sleep above baseline occurs on the first night. In general, the longer the deprivation, the smaller the fractional repayment of accumulated sleep debt, where we define sleep debt as the amount of sleep missed e.g., a 24 h deprivation gives a sleep debt of one night’s sleep duration. As discussed above, this is expected given the saturating growth of H, and is consistent with experiments showing amount of recovery sleep saturates for very long deprivations. To examine the model’s sensitivity to w, we model recoveries from 60, 120, and 180 h of deprivation, while varying w across its

13

12

Sleep (h)

(iii)

11

1st night 10

2nd night 3rd night

9

8 40

60

80

100

120

140

160

180

200

Deprivation (h) 4. Sleep debt The model predicts certain dynamics for the recovery of sleep debt (the amount of sleep missed during deprivation),

Fig. 6. Amount of sleep on the first three nights of recovery as a function of deprivation length, shown solid (1st night), dashed (2nd night), and dotted (3rd night). All three curves asymptote to maximal values. All parameters take on their new nominal values from Table 1, and the baseline sleep duration without deprivation is 8.6 h.

ARTICLE IN PRESS

60 h 0.3

120 h

120 h 180 h

0.2 0.15

4

0.1

3

0.05

2 2 Night

3

Fig. 7. Cumulative recovery of sleep across the first three nights as a fraction of the amount of sleep missed. Points are plotted for deprivations of 60, 120, and 180 h as circles, triangles, and squares, respectively. All parameters take on their new nominal values from Table 1. By the third night, the amount of recovery sleep has saturated to within 1% of its ultimate value.

constrained range, with all other parameters fixed at their nominal values. In each case, recovery begins at t ¼ 0 with a ¼ 17 h, corresponding to the normal bed time. Fig. 8 shows that for a given length of deprivation, recovery sleep increases with increasing w. Larger values of w also give greater discrimination between long duration deprivations, since H saturates more slowly. This is most clearly seen by looking at the difference in recovery times at either extreme of w. For w ¼ 34 h, the recovery time is little different between the three cases, in particular for the 120 and 180 h deprivations. For w ¼ 72 h, the three recoveries are clearly differentiated. 4.2. Recovery during wake While sleep is traditionally viewed as the only phase of recovery from sleep deprivation, this model suggests that wake also plays a role in recovery of H to normal levels, but not necessarily subjective fatigue. The same is true of two-process models with exponential dynamics for Process S (Daan et al., 1984), and is a consequence of the reasonable assumption that H is cleared at an increased rate when its concentration is higher. During the normal sleep–wake cycle, production of H is greater during wake than during sleep, and clearance is proportional to H. During wake, production exceeds clearance, so H grows, with the converse true during sleep. During recovery from sleep deprivation, the rates of production remain more or less the same, but the rate of clearance is significantly increased. Hence, during periods of wake following sleep deprivation, although H is still net increasing, the rate of increase is lower than in normal wake. In this sense, wake following deprivation yields a relative recovery. This can be seen by considering the return of H to baseline levels, as in Fig. 9. The most effective phases for returning H to baseline levels, are those when recovery sleep occurs at a time when wake would normally be expected, due to lower Q m. However, sleep during normal sleeping hours is no more effective at returning H to baseline levels than wake during normal waking hours. In summary, both wake and sleep are more effective at clearing H after deprivation than before. This result is discussed in Section 7.

1

40

50 χ (h)

60

70

Fig. 8. Amount of sleep in excess of baseline levels on the first two nights after (i) 60 h deprivation (solid lines), (ii) 120 h (dashed lines), (iii) 180 h (dotted lines), as a function of homeostatic clearance time constant w, with other parameters taking the nominal values in Table 1. In each case, the top line is the first night, and the bottom line is the second night.

Clock time 23:00

11:00

23:00

11:00

0 −10 −20

Normal 40 h

Δ H (nM)

1

m

0

60 h

5 Sleep (h)

0.25

419

180 h

6

V (mV)

Cumulative fraction of lost sleep recovered

A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

80 h 4

120 h

2 0 0

12

24

36

48

t (h) Fig. 9. (a) Time series for V m for normal (solid), and recoveries from 40 (dashed), 80 (dotted), and 120 h (dash-dotted) of deprivation. (b) Recovery of H to baseline levels after 40 (solid), 80 (dashed), and 120 h (dotted) of deprivation. All time series begin at the phase a ¼ 17 h, with approximate clock time on the top axis. H recovers monotonically, but most rapidly during recovery sleep that occurs during normal waking hours. All parameters have the new nominal values from Table 1.

5. Circadian phase dependence of recovery We now consider how the recovery process depends on the circadian phase at which deprivation is ceased, which is parametrized by a in (5). The amount of recovery sleep in excess of baseline sleep duration that is required to return to a normal cycle, is shown in Fig. 10, as a function of a. We define a recovery to be more ‘efficient’ if there is less total recovery sleep, giving us a measure by which to optimize recovery. Results differ based on the length of deprivation, but the model makes some general predictions. First, the amount of recovery sleep required to return to a normal sleep cycle is shortest when begun at the normal

ARTICLE IN PRESS 420

A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

14:00

Clock time 02:00 08:00

5

20:00

14 12

−5 −10 −15

m

V (mV)

Recovery sleep (h)

0

120 h 80 h 40 h

10

−20 −25

8 −30 −35 0

6

5

10

15

20

25

30

35

t (h)

4 0

6

12 α (h)

18

24

Fig. 10. Amount of recovery sleep as a function of the circadian phase a at which it begins, using the new nominal parameters in Table 1. Recoveries are from 120 h (solid), 80 h (dashed), and 40 h (dotted) of continuous wake. Approximate clock time is given on the top axis, which runs in reverse chronological order because an increased value of a corresponds to an earlier clock time.

bedtime, a ¼ 17 h, since recovery occurs in phase with the circadian cycle. The model predicts that this is the optimal phase at which to begin recovery, for deprivations of all durations, in terms of efficiency defined above. Second, there are two peaks in the amount of recovery sleep as a function of circadian phase, representing the least efficient times at which to begin recovery. One of these corresponds to beginning recovery at the normal waking time of a ¼ 8 h, since in this case recovery is out of phase with the circadian cycle. The other peak, which is a sharp transition, occurs at different values of a for different lengths of deprivation. It corresponds to a very different phenomenon, which occurs when recovery starts during the daytime. Depending on the amount of deprivation, the model predicts that it is possible either to: (a) briefly wake up due to circadian influences after a short restorative period of sleep, then sleep again at night, or (b) to simply sleep right through until the next morning. This peak represents the boundary between these two regimes, with sleeping through to the left of the peak, and brief awakening to the right of the peak. This behavior is demonstrated in Fig. 11, where a very small increase in a about this point, results in a significant change in behavior; noting that this is in absence of noise or extraneous stimuli which might smear the transition. This is because increasing a begins the recovery at an earlier circadian phase, so a brief awakening is possible later in the day due to the lower H value. Fig. 10 shows that for longer deprivations, this peak occurs at larger a, because a longer sleep is required before circadian influences are sufficient to induce wake. Sleeping through is predicted by the model to be a less efficient method of recovery, because wake shortly after deprivation can also act as a restorative phase in terms of H returning to baseline levels, as shown in Section 4.2. For the 40 h deprivation, the least efficient phase to begin recovery is the waking time a ’ 8 h, but for the longer deprivations, sleeping through becomes the least optimal. These predictions should be testable experimentally.

Fig. 11. Time series of V m during recovery from sleep deprivation, using the new nominal parameter values from Table 1. A small change in the initial circadian phase causes a sharp transition between briefly waking during the day at a ¼ 5:148 h (solid), and ‘sleeping through’ at a ¼ 5:147 h (dotted). Note that this is in the absence of noise which might smear the transition.

6. Comparison with phenomenological models Certain aspects of sleep deprivation and recovery have been modeled phenomenologically on numerous occasions, with most such models being based on the two-process model of sleep regulation (Achermann and Borbe´ly, 2003). We here compare the dynamics of our model with that of the two-process model in Section 6.1. We then demonstrate that our model can predict sleep latencies in Section 6.2, as an example of the additional capabilities of physiologically based modeling. 6.1. The two-process model The drives entering the VLPO in our model are very much in the two-process mold; however, there are important differences in their form. In our model, increasing D takes us from a zone where there is only a wake state (DoD1 ), to a bistable zone where both wake and sleep states exist (D1 oDoD2 ), and finally to a zone where there is only a sleep state (D2 oD). The key feature is that transitions from wake to sleep occur at D ¼ D2 , whereas transitions from sleep to wake occur at D ¼ D1 . These different thresholds result in hysteresis behavior, as shown in Fig. 12. An important difference between our model and the two-process model is that these thresholds for transition are a natural consequence of the physiology, and not artificially imposed. The two-process model is similar in that transitions occur when the homeostatic drive—or ‘‘Process S’’ in Achermann and Borbe´ly’s terminology (Borbe´ly and Achermann, 2000)—reaches one of two time varying oscillatory thresholds, as illustrated in Fig. 13(a). These thresholds represent circadian variations in wakefulness, and we call them T 1 and T 2 here. We can put the two-process model in a form analogous to our model by subtracting the circadian oscillation from the three functions T 1 , T 2 , and Process S. This is shown in Fig. 13(c) to yield constant thresholds, D1 and D2 , and a function which we call the total sleep drive, which moves between the thresholds. As in our model, the total sleep drive is effectively the weighted sum of homeostatic (Process S) and circadian drives, where the weightings nvc and nvh are of opposite signs. Alternatively, it is possible to put our model in the two-process form by defining oscillatory thresholds T 1 and T 2 by nvc C plus

ARTICLE IN PRESS A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

6.2. Sleep latencies

Wake

Vm (mV)

0

−5

−10 Sleep −15

D

D

1

0

2

1

2

3

D (mV) Fig. 12. Plot of V m against the total drive, D. As D oscillates across the day, V m makes transitions between the wake and sleep states, evolving around a hysteresis loop.

Two−process

Phillips and Robinson

(a) Process S

Process S

(b) T2 T

T

2

T1

1

(d) D2 D

1

12 24 36 48 t (h)

Total drive

Total drive

(c)

0

421

D

2

D

1

0

The most important distinction between our model and the two-process model is the inclusion of brainstem physiology in the former. This allows us to make direct inferences between model behavior and physiology—such as the thresholds and hysteresis behavior discussed above—and to directly relate model parameters to dynamics. An example is sleep latencies, which are beyond the scope of the two-process model, since it switches instantaneously between wake and sleep states. Previously, sleep latencies were estimated by measuring how long the model took to move between different values of V m (Phillips and Robinson, 2007). This was found to depend strongly on the values of tv and tm . A more precise definition of the sleep latency time, is the time required to move from the stable wake state to the stable sleep state. Here we take the model to be in a stable state when 2 2 vo0:1 mV min1 , where v ¼ ðV_ v þ V_ m Þ1=2 is the velocity in V v 2V m 1 space. Times when v40:1 mV min then correspond to transitions between wake and sleep, as shown in Fig. 14 by the shaded regions. The points at which v intersects 0:1 mV min1 , correspond to the beginning or end of a transition; i.e., leaving or arriving at a fixed point, respectively, and are also illustrated in Fig. 14. For a sleep latency of 15 min, we find tv ¼ tm ’ 20 s, doubling the previous nominal value (Phillips and Robinson, 2007). This change has a negligible effect on the overall recovery process, as the vast majority of time is still spent in either wake or sleep. However, it does allow us to make direct comparison with experimental sleep latency data. Using the above value of tv ¼ tm ¼ 20 s we simulate the sleep deprivation and recovery experiment of Lamond et al. (2007). Deprivations of 39 and 63 h are modeled, followed by a recovery process with a restricted daily duration for time in bed (TIB). This is achieved by solving the problem piecewise, holding the model in a wake state during the other hours of the day, regardless of size of the drive to sleep. We use a fixed waking value of V m ¼ 1:15 mV to produce 8.25 h normal sleep duration when allowed 9 h TIB, so as to match the average reported sleep duration for the subjects in the study. Sleep latency times are measured upon falling asleep at the end of each day, whereas in Lamond et al. they were measured during the morning and evening and averaged. For 9 h TIB the recovery behavior is similar for both deprivations, shown in Fig. 15. Sleep latency is reduced to about 5 min on the first night of recovery due to the loss of stability of the waking

12 24 36 48 t (h)

Fig. 13. An illustration of our model dynamics in a form similar to the two-process model in (a), alongside the two-process model in (b). Oscillatory thresholds for transitions are labeled as T 1 and T 2 . The homeostatic drive (Process S) increases during wake until T 2 is reached, when there is a transition to sleep. It then declines during sleep until a transition to wake occurs at T 1 . Equivalently, we show the total sleep drive for our model in (c), alongside the two-process model in our form in (d). The thresholds for state transitions are now constant, and labeled D1 and D2 .

appropriate constant offsets, and this is shown in Fig. 13(b). The basic dynamics of the two models are therefore qualitatively similar, and they share hysteresis behavior. The intuitive form of the two-process model is therefore justified by our physiological modeling approach. A key difference is that in the two-process model, the total drive remains between the two thresholds at all times, which is dynamically equivalent in our model to remaining in the bistable zone of the hysteresis loop, D1 oDoD2 , at all times, as shown in Fig. 12. In the presence of stimuli this can give unrealistic behavior, such as maintenance of wake without the need for prolonged external stimuli following extreme deprivations, and this is discussed in Section 7.

2 2 Fig. 14. Time series of V m , with circles indicating v ¼ ðV_ v þ V_ m Þ1=2 ¼ 0:1 mV min1 . These points correspond to leaving or arriving at a fixed point, with the time taken to move between the points labeled A and B being the sleep latency time. Times where v40:1 mV min1 are shaded, corresponding to transitions between wake and sleep states.

ARTICLE IN PRESS 422

A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

9 h TIB 6 h TIB

9 h TIB* 7 h TIB

*long deprivation

20

20 Experimental sleep latency (min)

Model sleep latency (min)

(a) 15

10

5

0

B

N1 Day

N2

(b) 15

10

5

0

B

SD Day

R1

Fig. 15. Sleep latencies following deprivation for (a) model simulations, and (b) experimental data (mean  SEM). The triangles are for 9 h TIB following 63 h sleep deprivation. All other curves are following 39 h sleep deprivation, with 9 h TIB (squares), 6 h TIB (circles), plus 7 h TIB (dotted line) for which there is no comparable experimental data. In (a), sleep latencies are measured at the end of each day, for baseline (B), and each night of recovery (N1, N2, etc.). In (b), sleep latencies are measured across the day, for baseline (B), during deprivation (SD), and during each recovery day (R1, R2, etc.). Data in (b) are adapted from Lamond et al. (2007).

state under higher sleep pressure, and most of the recovery to baseline occurs within 3–4 days as H is cleared. For 6 h TIB, H is never fully cleared, so the experiment becomes a chronic deprivation protocol where sleep is deprived every night. In this case, sleep latency does not recover to baseline, and instead asymptotes to a lower value of about 6 min. Intermediate cases are shown for 7 and 8 h TIB. All of these results match the data semiquantitatively (Lamond et al., 2007), demonstrating the ability of our model to reproduce realistic behavior and relate dynamics to physiological parameters; in this case, tm and tv . Sleep latencies are very similar between each experimental protocol and its model simulation, and the time course for recovery is correct to within one day. The failure to fully recover under the 6 h TIB protocol is verified experimentally by the Lamond et al. data, although our 7 h TIB protocol, which is shown for comparison, more closely fits the experimental case. It falls within error bars for all points except the first night, which is consistently higher for the model deprivations than for the experimental data.

7. Discussion We have applied a quantitative model of the ascending arousal system to sleep deprivation and recovery protocols. By calibrating the model response to sleep deprivation using a small number of criteria derived from human experimentation in the literature, the drive parameters have been better constrained. Importantly, these refinements do not significantly affect the normal sleep–wake cycle, so this paper remains consistent with our findings on normal sleep dynamics (Phillips and Robinson, 2007). Once calibrated to give a sensible response, the model was able to make new quantitative predictions about the recovery process, including optimal recovery plans.

The constraints applied to the model in Section 3 give more realistic ranges on the drive parameters m, w, nvh , nvc , and c0 . While only small adjustments were made to the nominal values of m and nvc , the new values of w, nvh , and c0 differ markedly from the old values, as presented in Table 1. The primary change to the model as a result of this calibration process has been an increased dependence on the homeostatic drive. In a previous work the model appeared to be circadian dominated, whereas the improved estimates here show that nvh is significantly higher. This change also implies increased c0 to maintain a realistic wake-to-sleep ratio in the absence of deprivation. Meanwhile, the increase in w causes sleep recovery to occur on the realistic scale of 1–3 days, and prevents H from saturating too rapidly during deprivation, so that deprivations of different lengths are suitably differentiated by model response. With the new parameter set the model response to sleep deprivation was realistic, allowing new model predictions to be made. In Section 4.1, we found that the amount of recovery sleep was consistent with experimental data, including the saturation effect on first night sleep duration after extreme deprivations. The distribution of recovery was also correct, with a return to normal sleep–wake patterns after no more than 3 days. Section 5 demonstrated that the efficiency of the recovery process, in terms of total length of recovery sleep above normal levels, is strongly dependent on the circadian phase at which recovery begins, and the model makes testable predictions of optimal recovery plans. For recovery from 120 h of deprivation, the difference between the most and least efficient phases is approximately 7 h of recovery sleep. Another striking prediction of the model is the sharp transition in total amount of recovery sleep between sleeping through and briefly awakening, where the latter can reduce the total amount of recovery sleep by over 1 h. These results are consistent with experimental findings that recovery sleep at an adverse circadian phase is either extended or shortened, depending upon whether the circadian peak is sufficient to interrupt the sleep period or not (Strogatz et al., 1986; Klerman et al., 1999). These model findings have important practical applications, in particular to optimizing hours for jobs involving long shifts or night shift work. Two results in this paper challenge the notion of sleep debt as a robust concept. In Section 4.1 it was found that a smaller fraction of missed sleep is repayed for longer deprivations. It was then shown in Section 4.2 that H in fact decreases relative to baseline levels during both the wake and sleep phases of recovery from deprivation. This finding is consistent with experiment, where if recovery sleep durations are held to 8 h, recovery is still achieved within a couple of nights, even after extreme deprivation (\200 h). Instead there is a redistribution of sleep stages, consistent with the deeper sleep predicted by the model, and less net production during wake due to the greater rate of clearance. It therefore seems clear that the concept of sleep debt as a fixed amount of sleep that must be repaid in the future is not correct. Indeed, as a result of this greater rate of clearance, recovery is possible even if TIB is restricted to normal durations (i.e., no extra sleep at all). However, it may be that some sleep deprivation-induced changes have a longer timescale of accumulation and recovery than the homeostatic drive, which could explain why some performance measures typically take longer to recover following chronic partial sleep deprivation than the usual 1–3 days for total sleep deprivation recovery (Van Dongen et al., 2003). This is an area for future study, with cortical synaptic density being one possible candidate for longer timescale cognitive changes (Tononi and Cirelli, 2006). The comparison between our model and the two-process model in Section 6.1 showed that the basic qualitative dynamics of the models have many similarities. The two-process model was

ARTICLE IN PRESS A.J.K. Phillips, P.A. Robinson / Journal of Theoretical Biology 255 (2008) 413–423

put into our form, demonstrating that both models have hysteresis behavior; but in the two-process model, this is the result of arbitrary conditions that are inspired by the dynamics. Our model gives a physiological explanation for this behavior, and is able to relate it directly to the pathology narcolepsy (Phillips and Robinson, 2007). Interestingly, the behavior of the two-process model was found to be dynamically equivalent to always remaining in the bistable zone of our model. This would mean that a short stimulus to the MA group could move the model from sleep to wake at any point during sleep, without the need for prolonged stimulus to remain awake. This does not realistically model the difficulties associated with altering one’s sleep cycle from day to day, and it is therefore required that most time be spent outside the hysteresis loop, as is the case with our model parameters. Additionally, our model is able to make predictions that lie beyond the scope of the two-process model, such as sleep latency times. The semiquantitative reproduction of Lamond et al.’s (2007) data in Section 6.2 provides further confirmation of the form of our model. In the future, fitting model parameters to data for each subject will enable even better predictions that cater for inter-individual differences, and a full sensitivity analysis may give insights into the underlying causes for some pathological behaviors. In future development of the model it may be necessary to examine the form of the H dynamics in greater detail to consider other applications. For example, the production rate of H is currently considered proportional to Q m ; this linear approximation being justified by the fact that the majority of time is spent in one of two states (wake and sleep), within both of which Q m is slowly varying, so there are essentially two points to fit. When more complex phenomena are considered, such as apnea, and partial awakenings in response to stimuli, it will become important to understand how the production of H varies as a function of Q m in more detail. In the longer term, the homeostatic drive may also need to be modeled at the pharmacokinetic level of adenosine A1 and A2a receptors; with the VLPO receiving input from nvb Q b ðHÞ, representing the firing rate of the basal forebrain neuronal population that adenosine targets, rather than simply nvh H. Currently, however, such a form is unnecessarily complicated for the dynamics being studied, and it is sufficient to assume the general form of production and clearance processes, without such details of the mechanisms involved.

Acknowledgment The Australian Research Council supported this work. References Achermann, P., Borbe´ly, A., 2003. Mathematical models of sleep regulation. Front. Biosci. 8, 683–693. Achermann, P., Dijk, D., Brunner, D., Borbe´ly, A., 1993. A model of human sleep homeostasis based on EEG slow-wave activity: quantitative comparison of data and simulations. Brain Res. Bull. 31, 97–113. Aston-Jones, G., Chiang, C., Alexinsky, T., 1991. Discharge of noradrenergic locus coeruleus neurons in behaving rats and monkeys suggests a role in vigilance. Prog. Brain. Res. 88, 501–520. Bonnet, M., 2000. Sleep deprivation. In: Kyger, M., Roth, T., Dement, W. (Eds.), Principles and Practice of Sleep Medicine, third ed. Saunders, Phil., London, pp. 53–71 (chapter 5).

423

Borbe´ly, A., Achermann, P., 2000. Sleep homeostasis and models of sleep regulation. In: Kyger, M., Roth, T., Dement, W., (Eds.), Principles and Practice of Sleep Medicine, third ed. Saunders, Phil., London, pp. 377–390 (chapter 29). Daan, S., Beersma, D., Borbe´ly, A., 1984. Timing of human sleep: recovery process gated by a circadian pacemaker. Am. J. Physiol. 246, R161–R178. Dı´az-Mun´oz, M., Herna´ndez-Mun´oz, R., Sua´rez, J., Vidrio, S., Ya´n´ez, L., AguilarRoblero, R., Oksenberg, A., Rosenthal, L., Villalobos, L., Ferna´ndez-Cancino, F., Drucker-Colı´n, R., de Sa´nchez, V., 1999. Correlation between blood adenosine metabolism and sleep in humans. Sleep Res. Online 2, 33–41. Diniz Behn, C., Brown, E., Scammell, T., Kopell, N., 2007. Mathematical model of network dynamics governing mouse sleep–wake behavior. J. Neurophysiol. 97, 3828–3840. Freeman, W., 1975. Mass Action in the Nervous System. Academic Press, NY. Gooley, J., Schomer, A., Saper, C., 2006. The dorsomedial hypothalamic nucleus is critical for the expression of food-entrainable circadian rhythms. Nat. Neurosci. 9, 398–407. Johnson, L., Slye, E., Dement, W., 1965. Electroencephalographic and autonomic activity during and after prolonged sleep deprivation. Psychosom. Med. 27, 415–423. Jones, B., 2000. Basic mechanisms of sleep–wake states. In: Kyger, M., Roth, T., Dement, W. (Eds.), Principles and Practice of Sleep Medicine, third ed. Saunders, Phil., London, pp. 134–154 (chapter 10). Kandel, E., Schwartz, J., Jessell, T., 1991. Principles of Neural Science, fourth ed. McGraw-Hill, NY. Klerman, E., Boulos, Z., Edgar, D., Mistlberger, R., Moore-Ede, M., 1999. Circadian and homeostatic influences on sleep in the squirrel monkey: sleep after sleep deprivation. Sleep 22, 44–59. Lamond, N., Jay, S., Dorrian, J., Ferguson, S., Jones, C., Dawson, D., 2007. The dynamics of neurobehavioral recovery following sleep loss. J. Sleep Res. 16, 33–41. Mallis, M., Mejdal, S., Nguyen, T., Dinges, D., 2004. Summary of the key features of seven biomathematical models of human fatigue and performance. Aviat. Space Environ. Med. 75, A4–A14. McKenna, J., Dauphin, L., Mulkern, K., Stronge, A., McCarley, R., Strecker, R., 2003. Nocturnal elevation of extracellular adenosine in the rat basal forebrain. Sleep Res. Online 5, 155–160. Mistlberger, R., 2005. Circadian regulation of sleep in mammals: role of the suprachiasmatic nucleus. Brain Res. Rev. 49, 429–454. Nunez, P., 1995. Neocortical Dynamics and Human EEG Rhythms. Oxford University Press, NY. Pace-Schott, E., Hobson, J., 2002. The neurobiology of sleep: genetics, cellular physiology and subcortical networks. Nat. Rev. Neurosci. 3, 591–605. Phillips, A., Robinson, P., 2007. A quantitative model of sleep–wake dynamics based on the physiology of the brainstem ascending arousal system. J. Biol. Rhythms 22, 167–179. Porkka-Heiskanen, T., Strecker, R., McCarley, R., 2000. Brain site-specificity of extracellular adenosine concentration changes during sleep deprivation and spontaneous sleep: an in vivo microdialysis study. Neuroscience 99, 507–517. Robinson, P., Rennie, C., Wright, J., 1997. Propagation and stability of waves of electrical activity in the cerebral cortex. Phys. Rev. E 56, 826–840. Robinson, P., Rennie, C., Rowe, D., O’Connor, S., 2004. Estimation of multiscale neurophysiologic parameters by electroencephalographic means. Hum. Brain Mapp. 23, 53–72. Saper, C., Chou, T., Scammell, T., 2001. The sleep switch: hypothalamic control of sleep and wakefulness. Trends Neurosci. 24, 726–731. Saper, C., Scammell, T., Lu, J., 2005. Hypothalamic regulation of sleep and circadian rhythms. Nature 437, 1257–1263. Sherin, J., Shiromani, P., McCarley, R., Saper, C., 1996. Activation of ventrolateral preoptic neurons during sleep. Science 271, 216–219. Strogatz, S., Kronauer, R., Czeisler, C., 1986. Circadian regulation dominates homeostatic control of sleep length and prior wake length in humans. Sleep 9, 353–364. Tamakawa, Y., Karashima, A., Koyama, Y., Katayama, N., Nakao, M., 2006. A quartet neural system model orchestrating sleep and wakefulness mechanisms. J. Neurophysiol. 95, 2055–2069. Thakkar, M., Strecker, R., McCarley, R., 1998. Behavioral state control through differential serotonergic inhibition in the mesopontine cholinergic nuclei: a simultaneous unit recording and microdialysis study. J. Neurosci. 18, 5490–5497. Tononi, G., Cirelli, C., 2006. Sleep function and synaptic homeostasis. Sleep Med. Rev. 10, 49–62. Van Dongen, H., Rogers, N., Dinges, D., 2003. Sleep debt: theoretical and empirical issues. Sleep Biol. Rhythms 1, 5–13. Webster, R., 2002. Neurotransmitters, Drugs and Brain Function. Wiley, Chichester. Williams, H., Lubin, A., Goodnow, J., 1959. Impaired performance with acute sleep loss. Psychol. Monogr. 73, 1–26. Yin, W., 2007. A mathematical model of the sleep–wake cycle. Masters Thesis, Georgia Institute of Technology.