Alertness Optimization for Shift Workers Using a Physiology-based Mathematical Model

Alertness Optimization for Shift W ork ers Using a Ph ysiology-based Mathematical Mo del Zidi T ao 1,2* , A. Agung Julius 1,2 , John T. W en 1,2 , 1 Dept. Electrical, Computer, and Systems Engineering, Rensselaer P olytec hnic Institute, T ro y , NY, USA 2 Ligh ting Enabled Systems and Applications (LESA) Engineering Research Cen ter, Rensselaer P olytechnic Institute, T roy , NY, USA * Corresp onding Author: taoz2@rpi.edu Abstract Sleep is vital for main taining cognitive function, facilitating metabolic waste remo v al, and supp orting memory consolidation. How ev er, mo dern so cietal demands, particularly shift w ork, often disrupt natural sleep patterns. This can induce excessiv e sleepiness among shift w orkers in critical sectors suc h as healthcare and transp ortation and increase the risk of accidents. The primary contributors to this issue are misalignmen ts of circadian rh ythms and enforced sleep-wak e sc hedules. Regulating circadian rhythms that are tied to alertness can b e regarded as a con trol problem with con trol inputs in the form of light and sleep sc hedules. In this pap er, w e address the problem of optimizing alertness b y optimizing light and sleep sc hedules to impro ve the cognitiv e p erformance of shift w ork ers. A key tool in our approac h is a mathematical mo del that relates the con trol input v ariables (sleep and ligh ting sc hedules) to the dynamics of the circadian clo c k and sleep. In the sleep and circadian mo deling literature, the new er ph ysiology-based mo del sho ws b etter accuracy in predicting the alertness of shift w orkers than the phenomenology-based mo del, but the dynamics of ph ysiological-based mo del ha v e 1/35 differen tial equations with different time scales, whic h p ose c hallenges in optimization. T o o vercome the c hallenge, w e prop ose a h ybrid version of the PR model by applying singular p erturbation tec hniques to reduce the system to a non-stiff, differen tiable h ybrid system. This reformulation facilitates the application of the calculus of v ariation and the gradien t descent method to find the optimal light and sleep sc hedules that maximize the sub jective alertness of shift w ork er. Our approac h is v alidated through n umerical simulations, and the sim ulation results demonstrate impro ved alertness compared to other existing sc hedules. 1 In tro duction Sleep is an imp ortan t routine of our daily life. During sleep, our brain remov es metab olic w aste [1, 2] and supp orts the formation of long-term memories [3]. In preindustrial so ciet y , human routines w ere mainly determined by the natural ligh t-dark cycle. In mo dern so ciet y , our habits are c hanged by the a v ailabilit y of artificial lighting, and shift work has become an essential part of our so ciet y that requires 24-hour service. Ho wev er, sleep depriv ation leads to increased sleepiness among shift w ork ers, which in turn results in acciden ts and work-related injuries [4, 5]. This can ha ve detrimen tal consequences not only for the shift w orkers themselv es but also for those around them, esp ecially in fields suc h as healthcare, transp ortation, and public safet y , where shift w ork is common. Therefore, there is a critical need to address sleepiness during shift w ork and increase alertness of shift work ers. The misalignments of circadian rh ythms and enforced sleep-w ake sc hedules are the main sources of increased sleepiness during shift work [6, 7]. Adjusting circadian rhythms can b e regarded as a control problem of a system with nonlinear dynamics [8]. The input of the system can be in the form of ligh t, c hemicals, and sleep schedules [9 – 11]. Optimal control theory pro vides wa ys to find the con trol inputs for a dynamic system such that the alertness of shift w ork ers is maximized [12, 13]. T o find the optimal control input, w e need mathematical mo dels of circadian rh ythms and sleep-wak e cycles. Man y sleep-wak e cycle mo dels are deriv ed from the Borb ely’s phenomenology-based t wo-process mo del [14]. The transition b et ween sleep and w ak efulness is regulated by t wo processes: homeostatic sleep pressure, whic h increases during wak efulness and 2/35 decreases during sleep, and circadian pacemak er, which is mainly affected b y the ligh t receiv ed by the retina. Achermann form ulated a three-pro cess mo del that com bines the t wo-process mo del and sleep inertia [15, 16]. The three-pro cess mo del w as used to predict cognitive function and alertness [17, 18]. Yin et al. used the three-pro cess mo del to calculate the ligh t and sleep schedules as con trol input to optimize the alertness of w orkers [13]. Ho wev er, as the ph ysiology underlying sleep-wak e dynamics is better understo od, newer models that incorp orate neuronal p oten tial-lev el interactions are dev elop ed [6, 19, 20]. Although the three-pro cess phenomenology-based mo del exhibits qualitativ e trends that match the data from the sleep depriv ation study , this mo del lac ks interpretation of ph ysiological asp ects. The Philips-Robinson (PR) mo del is a widely accepted ph ysiology-based mo del [21]. It asso ciates the phenomenological concepts of sleep depriv ations with the physiological mec hanisms of sleep using a mo del that in volv es the p oten tials of tw o brain stems: the VLPO group, which is sleep-activ e, and the MA group, whic h is wak e-activ e. The inclusion of ph ysiological mechanisms mak es the mo del sufficien tly flexible to describ e a wide range of phenomena, including alertness. W e p erform n umerical simulations of the three-process mo del and the PR mo del and compare the prediction of alertness with exp erimen tal data. W e observe that the PR mo del can impro ve the prediction of sub jectiv e alertness ov er the three-pro cess mo del. Ho wev er, the original PR mo del con tains stiff differen tial equations, which means that the pro cesses change on v ery different time scales. As recognized by Philips et al. [6], c hanges in neuronal p oten tials o ccur m uch faster than sleep homeostasis or circadian rh ythms. The stiff equations of system dynamics p ose c hallenges in optimization problems. The first c hallenge is the computational cost. Numerical solv ers for stiff systems often require small time steps to main tain stability and accuracy . Ho wev er, on a large time scale of m ultiple days, the computational costs increase significan tly . The second c hallenge is the calculation of the gradient. F or optimization tec hniques relying on gradien ts (e.g. gradien t descen t), the stiff equations complicate the computation of accurate and smo oth deriv atives. A recen t work b y P apatsimpa and Linnartz seeks to find the optimal light sc hedule that minimizes the difference b et w een natural wak e time and target w ake time with the PR model [22], but the optimality of the greedy algorithm in the pap er has not been prov en. Hong et al. applied the PR model to 3/35 analyze the sleep states of 21 shift work ers and create a personalized plan designed to reduce the sleepiness of shift w orkers [23]. Song et al. incorp orate the sleep sc hedule in [23] and prop ose a new schedule that aims to maximize the alertness of shift work ers. Ho wev er, these studies are based on empirical data rather than mathematical optimization. In order to use the PR mo del to find optimal con trol inputs that maximize the alertness of shift w orkers, w e prop ose a h ybrid v ersion of the PR mo del that has a h ybrid system. W e apply singular p erturbation techniques (see Chapter 11 in [24]). The system is reduced to a hybrid system that exhibits both contin uous and discrete dynamic b eha viors. The discrete sleep state determines the mo de of the system and the reduced mo del has con tin uous dynamics. The proposed hybrid PR model is non-stiff and differentiable so that w e can use calculus of v ariation tech niques to obtain the optimal solution. In this pap er, w e use the PR mo del to study the problem of optimizing the alertness of a shift w orker using ligh t and sleep sc hedules. W e also compare our optimal sc hedules with existing sleep sc hedule in [23] and demonstrate that our metho d can impro v e alertness in v arious scenarios. The main con tributions of this pap er are: • W e prop ose a hybrid PR model that is non-stiff and differentiable while k eeping the accuracy of alertness prediction of shift w orkers. • W e extend the optimization metho d in [13] with the hybrid PR mo del to optimize the alertness of shift w orkers using ligh t and sleep sc hedules as input. • W e compare our optimal schedules with other existing sc hedules that are based on empirical data and exp erimen tal studies. Through numerical sim ulations, w e demonstrate that our schedules can further improv e alertness during shift work or during the full optimization horizon. 4/35 2 Mo del Deriv ation and V alidation 2.1 Existing Mo dels 2.1.1 Three-Pro cess Mo del Ac hermann developed a widely recognized three-pro cess sleep regulation mo del, comprising three k ey comp onen ts: the circadian pacemak er (Pro cess C), the sleep homeostasis (Pro cess S), and sleep inertia (Process W) [15]. The new er mo del updates comp onen ts of the three-pro cess model to incorp orate new er discov eries [13]. Circadian pacemaker comp onent: The circadian pacemak er is represented b y the Jew ett-F orger-Kronauer (JFK) mo del [25]. The states x ( t ) and x c ( t ) are states of the core b ody tem perature oscillator. The state n ( t ) is the state of the light receptor in the retina. They follow the dynamics: dx dt = π 12 [ x c + µ ( 1 3 x + 4 3 x 3 − 256 105 x 7 ) + (1 − 0 . 4 x )(1 − k c x c ) u ] , (1) dx c dt = π 12 [ − ( 24 0 . 99729 τ ) 2 x + ( q x c − k x )(1 − 0 . 4 x )(1 − k c x c ) u ] , (2) u = Gα 0 ( I I 0 ) p (1 − β )(1 − n ) , (3) dn dt = 60[ α 0 ( I I 0 ) p (1 − β )(1 − n ) − γ n ] , (4) where all the parameters are listed in T able 1. Sleep Homeostasis comp onen t: The sleep homeostasis represents the accum ulation of sleep pressure that forms the sleep drive [26]. The Process S has state H ( t ) and it follo ws the dynamics: dH dt =        − H /τ d , β ( t ) = 1 , (1 − H ) /τ r , β ( t ) = 0 , (5) where β ( t ) = 1 represents the sub ject is asleep and β ( t ) = 0 represents the sub ject is a wak e and τ d = 4 . 2 h, τ r = 18 . 2 h are time constan ts. Sleep Inertia comp onent: The sleep inertia represen ts a decrease in the level of 5/35 men tal alertness immediately after aw ak ening [15, 16]. The process W has dynamics: dW dt = − 1 − β τ W W , (6) where the time constan t τ W = 0 . 662 h . The sleep w ake transition is go v erned by sleep propensity defined as: Φ( t ) = H ( t ) − A c x ( t ) , (7) where A c = 0 . 1333 . When Φ reac hes the low er threshold L m , the sub ject wak es up sp on taneously . Likewise, when Φ reac hes the upp er threshold H m , the sub ject falls asleep sp on taneously . Sleep state β ( t ) follows the discrete dynamics: β ( t ) = F β ( t, β ( t − )) =                1 , Φ( t ) = H m , 0 , Φ( t ) = L m , β ( t − ) , otherwise , (8) where β ( t − ) = lim τ − → t − β ( τ ), t − represen ts the time just b efore t , and constan ts H m = 0 . 67 , L m = 0 . 17 . The alertness is defined as A T P ( t ) = [1 − β ( t )][1 + A c x ( t ) − H ( t ) − W ( t )] , (9) and sleepiness lev el is defined as B ( t ) = 1 − A ( t ) . (10) T able 1. Three-process mo dels P arameter V alue P arameter V alue µ 0 . 13 h − 1 k c 0 . 4 h q 1 / 3 τ x 24 . 2 h k 0 . 55 h − 1 G 33 . 75 α 0 0 . 05 h − 1 p 0 . 5 I 0 9500 lux γ 0 . 0075 h − 1 τ d 4 . 2 h τ r 18 . 2 h τ W 0 . 662 h A c 0 . 1333 H m 0 . 67 L m 0 . 17 6/35 The states of the three-pro cess model are denoted as ξ ( t ) = [ x ( t ) , x c ( t ) , n ( t ) , H ( t ) , W ( t ) , β ( t )] T ∈ R 6 . (11) W e consider a 24-hour p eriodic light sc hedule as: I ref ( t ) =        150 lux , mo d( t, 24) ∈ [0 , 16) , 0 , mo d( t, 24) ∈ [16 , 24) , (12) whic h simulate the indo or en vironmen t like hospital. The initial time t = 0 corresp onds to 6 AM when the light is turned on and t = 16 corresp onds to 10 PM when the light is turned off. The sub ject go es to sleep and wak es up sp ontaneously . A t steady state, ξ ( t ) forms a stable p eriodic solution, denoted as: ξ ref ( t ) = [ x ref ( t ) , x c, ref ( t ) , n ref ( t ) , H ref ( t ) , W ref ( t ) , β ref ( t )] T , (13) in Fig. 1. Fig 1. P erio dic solution of three-process mo del accommo dates to the perio dic light sc hedule of 16-hour ligh t 8-hour dark. 48 hours of the state tra jectories are plotted here. T o sim ulate forced wak efulness, e.g., to study the effects of nigh t shift work and sleep depriv ation, w e ignore the sleep-switc hing dynamics in Eq. 8 and force the sleep state β ( t ) = 0 during forced wak efulness. 7/35 2.1.2 Philips-Robinson (PR) Mo del A widely used and more recen tly published sleep mo del is the Philips-Robinson (PR) mo del [21]. This mo del consists of three parts: a circadian pacemak er comp onen t, a ph ysiological comp onen t that mo dels t wo m utually inhibitory neuron cell potentials, and a sleep homeostasis comp onen t. Circadian pacemaker comp onent: The circadian pacemak er of the published PR mo del uses a circadian pacemaker model published in [27]. In this pap er, w e up date the PR mo del with the circadian pacemak er in [25]. The circadian pro cess has the same dynamics and parameters as Eq. 1-4 in three-pro cess mo del. Ph ysiological comp onen t: The physiological components of the PR model are tw o m utually inhibitory neuron cell p oten tials. The monoaminergic (MA) potential that promotes w akefulness is represen ted b y V m and the v entrolateral preoptic (VLPO) p oten tial that promotes sleep is represen ted by V v . The dynamics are defined by: dV m dt = 1 τ m [ − V m − v mv Q v + A m ] , (14) dV v dt = 1 τ v [ − V v − v v m Q m + D v ] , (15) D v = − v v c C + v v h H + A v , (16) where the firing rates Q m and Q v are sigmoidal functions of V m and V v : Q m = Q max 1 + exp( − V m − θ σ ) , Q v = Q max 1 + exp( − V v − θ σ ) . The term D v is the sleep driv e and C = 1 2 (1 + c x x + c x c x c ) represen ts the circadian comp onents. V m ( t ) is at a low lev el and V v ( t ) is at a high level during the sleep state, and vice v ersa during the wak e state. Sleep homeostasis comp onent: The sleep homeostasis in PR model has different 8/35 dynamics than the three-pro cess model: dH dt = 1 χ ( − H + µ H Q m ) . (17) Sleep inertia is omitted in the PR mo del for t w o reasons: the effect is short-liv ed and the mo del is directly based on ph ysiology , so the introduction of sleep inertia may lead to arbitrary fitting that degenerates the mo del. All parameters of the PR mo del are listed in T able 2. In [28], Postno v a et al. used a linear fitting pro cedure to align C and H with v arious cognitive performance metrics, using a distinct set of parameters for each metric. T o av oid ov erfitting to a single metric, w e use the definition of alertness in [29]: A ( t ) = [1 − β ( t )][ H + ( t ) − H ] . (18) All the states of PR mo del can be denoted as ζ ( t ) = [ x ( t ) , x c ( t ) , n ( t ) , V m ( t ) , V v ( t ) , H ( t )] T ∈ R 6 (19) Similar to three-pro cess model, we apply the same perio dic light sc hedule in Eq. 12. The resulting p eriodic solution after entrainmen t are denoted as ζ ref ( t ) = [ x ref ( t ) , x c, ref ( t ) , n ref ( t ) , V m, ref ( t ) , V v , ref ( t ) , H ref ( t )] T ∈ R 6 (20) In the PR mo del, whether the sub ject is asleep or aw ak e is determined by V m . If V m is higher than a threshold V m, th = − 3 . 785 mV that corresp onds to Q m > 1 s − 1 , the sub ject is aw ake. Since the time constants of the MA and VLPO dynamics in Eq. 14-15 are m uch faster than the circadian clock and sleep homeostasis, w e can use a quasi-static analysis and assume D v con verges to the steady state equilibrium. The equilibrium p oin t of V m as a function of D v can b e obtained b y setting the left-hand 9/35 (a) Circadian part (b) MA,VLPO and sleep homeostasis Fig 2. P erio dic solution of PR model accommo dates to the p eriodic light sc hedule in Eq.12. Fig. 2(a) includes the circadian states x, x c , n and Fig. 2(b) includes the neuron p oten tials V m , V v and sleep homeostasis H . side of Eq. 14 and Eq. 15 to 0. The analytical form of V m as a function of D v is: V m = − V mv ·       Q max 1 + exp V vm Q max 1+exp( V m − θ σ ) − D v − θ σ !       + A m (21) By analyzing the equilibrium p oin ts of V m for differen t D v in Eq. 14-15, we ha v e the bifurcation curv es in Fig. 3. Detailed analysis can b e found in [30]. If D v > 2 . 46, there is a unique equilibrium with lo w V m and therefore the sub ject is asleep. If D v < 1 . 45, there is a unique equilibrium with high V m and the sub ject is aw ak e. If 1 . 45 ≤ D v ≤ 2 . 46, there are t wo equilibria. The sub ject can b e a wak e or asleep. This region is kno wn as the bistable region. F rom the definition of D v in Eq. 16, the sleep and w ake threshold of H can b e deriv ed: H + = 2 . 46 − A v + v v c C v v h , (22) H − = 1 . 45 − A v + v v c C v v h . (23) T o sim ulate forced w ak efulness in the PR mo del, a driv e to V m is needed to main tain the w ake state, known as wake effort in [30 – 32]. W ake effort is an external stimulus that can take the form of a pharmacological agent or sensory input. With an additional w ake effort to V m , the sub ject will stay in the wake ghost due to the characteristics of the system dynamics, k eeping the sub ject a wak e. V m will b e set to the maximum of 10/35 V m, ref during p eriodic solution during forcedwak efulness. If H < H + , no w ake effort is needed to mak e the sub ject sta y aw ak e. If H ≥ H + and the sub ject is required to main tain wak e state, we keep V m at a fixed v alue, as the y elllow curve in Fig. 3. F orced sleep is also rep orted in [23, 29]. Ho w ever, it is generally difficult to induce immediate sleep in an individual. Therefore, we do not consider the case of forced sleep. Fig 3. Bifurcation curves of PR mo del. The x-axis shows the sleep drive D v and y-axis sho ws the corresp onding V m at equilibrium. In the wak e and sleep region, the wak e and sleep branc h ov erlap b ecause there is only one stable equilibrium. In the bistable branc h, there are tw o stable equilibrium p oints for each D v . When simulating sleep depriv ation, we keep V m = 1 . 04 mV . The main differences b et w een the three-pro cess mo del and the PR mo del are summarized in T able 4 in the next subsection. T able 2. PR mo dels P arameter V alue Parameter V alue µ 0 . 13 h − 1 k c 0 . 4 h q 1 / 3 τ x 24 . 2 h k 0 . 55 h − 1 G 33 . 75 α 0 0 . 05 h − 1 p 0 . 5 I 0 9500 lux γ 0 . 0075 h − 1 τ m 1 / 360 h v mv 1 . 8 mV s A m 1 . 3 mV τ v 1 / 360 h v v m 2 . 1 mV s v v c 3 . 37 mV s v v h 1 . 01 mV nM − 1 A v − 10 . 2 mV Q max 100 s − 1 θ 10 mV V m, th − 3 . 785 mV σ 3 mV c x 0 . 8 c x c − 0 . 16 µ H 4 . 2 nM s χ 45 h 11/35 2.2 Prop osed Mo del T o address the issues with the PR mo del, we prop ose a new v ersion of the PR mo del that has h ybrid b eha vior. T o distinguish the difference b et w een these tw o mo dels, we refer to the original PR mo del as ful l PR mo del and the hybrid version as hybrid PR mo del . W e use singular p erturbation to reduce the mo del. The slow timescale comp onen ts, i.e. circadian rh ythms and sleep homeostasis, are still mo deled b y the con tinuous dynamics in the JFK mo del. The fast timescale comp onents, suc h as V m and V v , are assumed to reac h steady states. The v alue of V m is deriv ed from Eq. 21 as a function of sleep driv e D v : V m = F ( D v ) and V v is neglected. Similar to the three-pro cess mo del, the sleep state β ( t ) is discrete and change state sp on taneously with the follo wing dynamics: β ( t ) = F β ( t, β ( t − )) =                1 , H ( t ) = H + , 0 , H ( t ) = H − , β ( t − ) , otherwise , (24) Dep ending on which of the three regions D v falls in, V m can b e determined using the bifurcation analysis in the previous subsection. F or the sleep or w ake region, V m can b e uniquely determined by Eq. 21. F or the bistable region, V m is the larger ro ot of Eq. 21 if the sub ject is aw ak e and the smaller ro ot of Eq. 21 if the sub ject is asleep. When the sub ject is in sleep depriv ation, then V m = max( V m, ref ) . W e can com bine the wak e region, the bistable region of the w ake branch, and the wak e ghost to the forced wak e branc h, shown in Fig. 4 a. The sub ject will alwa ys b e aw ak e on the forced w ake branch. With the sleep branc h and the forced wak e branc h, we can determine the v alues of V m using D v . In order to apply gradien t-based optimization metho ds, suc h as the calculus of v ariation, the gradient of the optimization ob jective needs to b e differentiable. W e use t wo differentiable functions to mo del the sleep branch and the forced wak e branc h. W e use the sum of a sigmoid function and a fifth-order p olynomial to fit the sleep branch in 12/35 Fig. 4 b. V m = V 1 ( D v ) = 3 . 5702 exp( − 40( D v − 1 . 45)) − exp(40( D v − 1 . 45)) exp( − 40( D v − 1 . 45)) + exp(40( D v − 1 . 45)) + 5 X i =0 a i D i v (25) Here a i are the co efficien ts of the p olynomials and are listed in T able 3. The co efficien t v alues are obtained using the p olyfit function in Matlab. The forced w ake branch is discontin uous at D v = 2 . 46 . T o mak e it smo oth, w e con volv e the forced wak e branch with a function in Fig. 4.c defined as follows: g ( D v ) =                0 , D v ≤ 0 , 468750 D v (0 . 2 − D v ) 4 , 0 < D v < 0 . 2 , 0 , D v ≥ 0 . 2 . (26) The filter is con tinuous and the integration R ∞ −∞ g ( D v ) = 1 After the con volution, the forced wak e branch is approximated as V m = V 2 ( D v ) =                P 5 i =0 b i D i v , D v ≤ 2 . 46 , P 11 i =0 c i D i v , 2 . 46 < D v ≤ 2 . 66 , 1 . 04 , D v > 2 . 66 . (27) Here b i and c i are all co efficien ts and are listed in T able 3. T able 3. P arameters for Sleep branc h mo del and F orced W ake Branc h Mo del P arameter V alue P arameter V alue a 0 − 3 . 2369 a 1 − 3 . 9232 a 2 9 . 2384 a 3 − 7 . 3438 a 4 2 . 0482 a 5 − 0 . 1964 b 0 1 . 1236 b 1 − 0 . 3960 b 2 0 . 8783 b 3 − 1 . 0640 b 4 0 . 5328 b 5 − 0 . 0982 c 0 1 . 2155 · 10 7 c 1 − 2 . 5973 · 10 7 c 2 2 . 2560 · 10 7 c 3 − 1 . 0007 · 10 7 c 4 2 . 2781 · 10 6 c 5 − 2 . 0589 · 10 5 c 6 − 1 . 0268 · 10 4 c 7 6 . 7223 · 10 3 c 8 − 3 . 4379 · 10 3 c 9 1 . 2007 · 10 3 c 10 − 217 . 0005 c 11 16 . 6128 13/35 (a) F orced W ak e Branch (b) Sleep Branch Fitting (c) Filter g ( D v ) (d) F orced W ak e Branch Fitting Fig 4. (a) F orced wak e branch. (b) Comparison of the real sleep branch (blue) and the fitted sleep branc h in hybrid PR mo del (red) (c) The filter function g ( D v ) b et w een 0 ≤ D v ≤ 0 . 2. (d) Comparison of the real forced wak e branch (red) and the fitted forced w ake branch in hybrid PR mo del(blue). 14/35 The new h ybrid system follows the dynamics: dx dt = π 12 [ x c + µ ( 1 3 x + 4 3 x 3 − 256 105 x 7 ) + (1 − 0 . 4 x )(1 − k c x c ) u ] , (28) dx c dt = π 12 [ − ( 24 0 . 99729 τ ) 2 x + ( q x c − k x )(1 − 0 . 4 x )(1 − k c x c ) u ] , (29) u = Gα 0 ( I I 0 ) p (1 − β )(1 − n ) , (30) dn dt = 60[ α 0 ( I I 0 ) p (1 − β )(1 − n ) − γ n ] , (31) dH dt = 1 χ − H + µ H Q max 1 + exp( − V m − θ σ ) ! , (32) V m =        V 1 ( D v ) , if β = 1 , V 2 ( D v ) , if β = 0 , (33) D v = − v v c C + v v h H + A v . (34) The mo del parameters in Eq. 28-32 and Eq. 34 are same as the parameters in full PR mo del. The alertness of h ybrid PR mo del can also b e defined as: A ( t ) = [1 − β ( t )][ H + ( t ) − H ] . (35) The states of the h ybrid PR mo dels are denoted as X ( t ) = [ x ( t ) , x c ( t ) , n ( t ) , H ( t ) , β ( t )] T ∈ R 5 . (36) When the ligh t schedule is the same as Eq. 12, and the sleep state β ( t ) in the h ybrid PR mo del follows Eq. 24, we obtain the p erio dic solution X ref ( t ) = [ x ref ( t ) , x cref ( t ) , n ref ( t ) , H ref ( t ) , β ref ( t )] T . (37) T o demonstrate the accuracy of our appro ximation using the prop osed h ybrid PR mo del, we compare the state tra jectories of the hybrid PR mo del to the full PR mo del. The p eriodic solutions of circadian rh ythm x ref , x cref and n ref are plotted in Fig. 5(a) and the sleep and wak e-up time are plotted in Fig. 5(b). The phase of circadian rhythm is defined as − tan( x c x ) and the av erage phase difference b et ween the hybrid PR mo del and full PR mo del is 0.0126 radians, which is equiv alen t to 2.89 minutes. The hybrid 15/35 (a) Comparison of Circadian Rh ythms (b) Comparison of Sleep Time (c) Comparison of H (d) Comparison of V m Fig 5. These are the comparisons of three mo dels. The x axes show the time in hours. The starting times are 6 AM and the durations are 48 hours. (A) Circadian rhythm Pro cess C which includes states x, x c and n . These three states hav e the iden tical dynamics in the full PR mo del and hybrid PR mo del. (b) Sleep state β ( t ). 0 means the sub ject is aw ak e and 1 means the sub ject is asleep. F or the full PR mo del, the sub ject w akes up at 6:33 AM and sleeps at 10:25 PM. F or the hybrid PR mo del, the sub ject w akes up at 6:27 AM and sleeps at 10:17 PM. (c) Comparison of the sleep homeostasis b et w een the full PR mo del and h ybrid PR mo del (d) Comparison of the V m b et w een the full PR mo del and hybrid PR mo del. mo del wak es up 6 minutes earlier and go es to sleep 8 minutes earlier than the full h ybrid mo del. The states H and V m are plotted in Fig. 5(c) and Fig. 5(d). W e can see that the h ybrid PR mo del is a go od approximation of the full PR mo del. 2.3 V alidation T o v alidate the precision of cognitive p erformance prediction, we p erform tw o sim ulations for all three mo dels discussed and fit the predicted alertness data to the exp erimen tal data in [33, 34]. In c onstant r outine exp eriment pr oto c ol [33], sub jects were instructed to remain aw ak e for an extended duration of 60 hours under contin uous ligh t conditions (150 lux). Sub jects access their sub jective alertness b y indicating their level 16/35 T able 4. Comparison of Three-Pro cess Mo del, PR Mo del and Hybrid PR mo del Three-pro cess Mo del PR mo del Hybrid PR Mo del Num b er of States 6 6 5 System Dynamics Hybrid Con tinuous Hybrid Stiffness Non-stiff Stiff Non-stiff of alertness on a linear 100-mm visual-analog scale (V AS) [35]. The exp erimental data ha ve the sub jectiv e alertness measurement at n different time steps and the mean and co v ariance are denoted as M V AS = [ µ V AS 1 , · · · , µ V AS n ] T (38) and Σ V AS =          σ 2 V AS 1 0 · · · 0 0 σ 2 V AS 2 · · · 0 . . . . . . . . . . . . 0 0 · · · σ 2 V AS n          . (39) W e use linear fitting to find the parameters θ V AS = [ θ V AS 1 , θ V AS 2 ] that minimizes the normalized mean square error (NMSE): N M S E V AS = ( M V AS − ( θ V AS 1 A + θ V AS 2 1 )) T Σ − 1 V AS ( M V AS − ( θ V AS 1 A + θ V AS 2 1 )) (40) where A = [ A 1 , A 2 , · · · , A n ] T represen ts the predicted alertness at time step 1 to n and 1 ∈ R n is a v ector of ones. The parameters and NMSEs for each mo del are listed in T able 5. In photop erio d exp eriment [34], sub jects were instructed to stay in 16 h dayligh t in one w eek and then stay in dim light less than 1 lux for 24 hours. The Stanford Sleepiness Scale [36] are recorded during the dim light p eriod. Similarly , the exp erimen tal data ha ve mean M S S S = [ µ S S S 1 , µ S S S 2 , · · · , µ S S S n ] T (41) 17/35 and Σ S S S =          σ 2 S S S 1 0 · · · 0 0 σ 2 S S S 2 · · · 0 . . . . . . . . . . . . 0 0 · · · σ 2 S S S n          . (42) W e also find the parameters θ S S S = [ θ S S S 1 , θ S S S 2 ] that minimizes the NMSE for SSS: N M S E S S S = ( M S S S − ( θ S S S 1 A + θ S S S 2 1 )) T Σ − 1 S S S ( M S S S − ( θ S S S 1 A + θ S S S 2 1 )) (43) The parameters for all three mo dels are listed in T able 6. T able 5. Alertness Linear Fitting Parameters Three-pro cess Mo del PR mo del Hybrid PR Mo del θ V AS 1 72.92 9.17 9.52 θ V AS 2 22.19 53.25 55.19 N M S E V AS 32.99 29.24 25.37 T able 6. Sleepiness Linear Fitting Parameters Three-pro cess Mo del PR mo del Hybrid PR Mo del θ S S S 1 -0.656 3.99 -0.656 θ S S S 2 2.72 0.057 2.78 N M S E S S S 29.22 8.61 9.55 W e plot the fitted sub jective alertness and sleepiness in Fig. 6. The hybrid PR mo del has the low est N M S E V AS among three mo dels. The N M S E S S S of PR mo del is also v ery close to the full PR mo del and significan tly low er than the three-pro cess mo del. The alertness A predicted b y the hybrid PR mo del shows close agreement with the empirical data in [33, 34]. 3 Optimal Ligh t and Sleep Sc hedule F or Impro ving Alertness Extensiv e research has aimed to enhance the alertness of shift work ers by implementing structured sleep schedules [29, 37, 38]. F or shift work ers, the sub ject must remain aw ake during sc heduled night shifts. In [13], Yin et al. proposed the problem of optimizing alertness with light and sleep schedules together with the three-pro cess mo del. Inspired 18/35 (a) Alertness Fitting Comparison (b) Sleepiness Fitting Comparison. Fig 6. Alertness and Sleepiness Fitting Comparison. The x-axis shows the time that sub ject stays aw ak e since wak e up time. The y-axis shows the sub jective alertness or sleepiness. The empirical data shows the mean and plus minus standard deviation of the exp erimen tal data. b y [13], we discuss tw o problems for optimizing alertness with the h ybrid PR mo del b ecause the PR mo del can predict alertness more accurately than the three-pro cess mo del. F or the first problem, w e find the optimal light and sleep schedules during the optimization horizon to maximize the cum ulative alertness only during the shift work. F or the second problem, w e find the optimal light and sleep schedules during the optimization horizon that maximize the cum ulative alertness during the full horizon. W e optimize not only alertness during shift work but also during the time the sub ject is resting. The details of these problems will b e defined in the following subsections. In b oth problems, light is sub ject to the constrain t: 0 ≤ I ( t ) ≤ 150 lux . (44) W e adopt the ‘tunable sleep sc hedule’ in [39]. The sub ject can stay a wak e when H + ≤ H ( t ) ≤ H + max where H + max = 3 . 43 − A v + v vc C v vh corresp onds to H sta ying aw ak e for 2 hours after reac hing the sleep threshold. The introduction of this upp er limit preven ts excessiv e sleep depriv ation. The sub ject can also wak e up earlier than the sp ontaneous w ake time when H − ≤ H ( t ) ≤ H − max where H − max = 2 . 11 − A v + v vc C v vh corresp onds to the v alue of H 2 hours b efore reaching the wak e threshold. W e assume that the sub ject is alw ays aw ak e ( β (0) = 0) at the b eginning of optimization. The sleep state during 19/35 tunable sleep sc hedule can b e expressed as β ( t ) = F ′ β ( t, β ( t − )) =                1 , t ∈ { t 1 sleep , · · · , t N f sleep } , 0 , t ∈ { t 1 wak e , · · · , t M f wak e } , β ( t − ) , otherwise , (45) with the comfort constrain ts H − ≤ H ( t i sleep ) ≤ H − max , H + ≤ H ( t j wak e ) ≤ H + max (46) where t i sleep and t j wak e represen t the i -th sleeping time and j -th w ake-up time. Here, N f and M f denote the num ber of times the sub ject falls asleep and wak es up, resp ectiv ely , during the optimization p eriod. I ( t ) , t 1 sleep , · · · , t N f sleep , t 1 wak e , · · · , t M f wak e are all optimization v ariables. Note that the shift work are predetermined by user’s demand, so the sleep state during shift w ork will alwa ys b e aw ak e and do es not need to satisfy the comfort constrain ts. 3.1 Shift W ork Alertness Optimization In this subsection, we solv e the shift work alertness optimization problem. The states of the sub ject X follo w the hybrid dynamics that is piecewise contin uous and the dynamics switc h b et ween asleep mo de ( β ( t ) = 1) and aw ak e mo de ( β ( t ) = 0). Assume that the dynamics for the i -th mo de is describ ed b y the following dynamics: ˙ X = D i ( X, I ) , t ∈ [ t i − 1 , t i ) , ∀ i ∈ { 1 , 2 , · · · , N } , (47) where N is the total num ber of mo des during the optimization time horizon. All switc hing times t i are sub ject to the inequality constraints in the tunable schedule in Eq. 46, except during shift work. The shift w ork alertness optimization problem is: Giv en the initial conditions X (0) = X 0 , the dynamics of the h ybrid mo del represen ted as ˙ X = D i ( X, I , β ) for i = 1 , · · · , N , the light constraint in Eq. 44 and tunable sleep sc hedule, we maximize the cumulativ e alertness of a sub ject during shift 20/35 w ork ov er multiple days: max Z t ∈ t wor k (1 − β ( t ))[ H + ( t ) − H ( t )] dt, (48) where t wor k represen ts the set of shift work interv als. W e apply the v ariational calculus metho d in [13] to solve the alertness optimization problem. An ob jectiv e function is formulated as J = Z t f 0 L ( t, X, I ) dt, (49) and L = − 1 wor k [ H + ( t ) − H ( t )], where 1 wor k ( t ) =        1 , if t ∈ t wor k , 0 , if t ∈ t wor k , (50) is an indicator function. W e use the standard constrained optimization approach by in tro ducing Lagrange multipliers λ ( t ) ∈ R 4 . The augmented ob jectiv e function is: J a = Z t f 0 L ( τ , X , I ) dτ + N X i =1 Z t i t i − 1 λ ( τ )[ D i ( X, I ) − ˙ X ( τ )] dτ . (51) The selection of λ do es not affect the v alue of J a . By selecting the Lagrange multiplier λ ( t ) as: λ ( t N ) = 0 , (52) ˙ λ ( t ) = − δ L ( t, X, I ) δ X −  δ D i ( X, I ) δ X  T λ ( t ), when t ∈ [ t i − 1 , t i ) , (53) λ ( t − i ) = λ ( t + i ) , (54) the first v ariation of augmented cost δ J a simplify to: δ J a = N X i =1 Z t i t i − 1 ∇ I ( t ) J δ I dτ + N − 1 X i =1 ∇ t i J δ t i , (55) ∇ I ( t ) J = ∂ L ( t, X , I ) ∂ I + λ T ( t ) ∂ D i ( X, I ) ∂ I , (56) 21/35 ∇ t i J = L ( t − i , X ( t − i ) , I ( t − i )) − L ( t + i , X ( t + i ) , I ( t + i )) + λ ( t − i ) D i ( X ( t − i ) , I ( t − i )) − λ T ( t + i ) D i +1 ( X ( t + i ) , I ( t + i )) (57) where ∇ I ( t ) J and ∇ t i J are the gradient of J with resp ect to I and t i and δ I , δ t i are p erturbations in I and t i . The details of λ and gradient formulations are in the App endix. Applying gradien t descent metho d, at the j -th iteration, the light input I j ( t ) and i-th swtic hing time t j i can b e up dated by: I j +1 ( t ) = min(max(0 , I j ( t ) − η I ∇ I ( t ) J ) , I max ) , (58) t j +1 i = t j i − η t ∇ t j i J ∈ Ω t i , (59) where η I and η t are the step size for up date and Ω t i represen ts the set of feasible t i that satisfy the constrain t Eq. 46. The steps of the gradien t descent algorithm can b e summarized as follo ws. Gradien t Descent Algorithm Step 1: Start with an initial guess of ligh t I 0 ( t ) and sleep sc hedule t 0 i , i ∈ { 1 , 2 , · · · , N − 1 } iteration j = 0; Step 2: In tegrate the state dynamics with Eq. 47 forward in time to obtain state tra jectory X j ( t ) and the ob jective function J j ; Step 3: In tegrate the Lagrange multiplier with Eq. 52 - 54 backw ard in time to obtain λ j ( t ) , t ∈ [ t 0 , t f ]; Step 4: Calculate the gradien ts ∇ I ( t ) j J and ∇ j t i with Eq. 56 and Eq. 57, then p erform up dates to I j +1 and t j +1 i with Eq. 58, 59; Step 5: Rep eat Step 2 to 4 un til light and sleep schedules reach conv ergence. In Step 1 of the gradien t descent algorithm, different initial guesses can lead to differen t lo cal minima. W e construct an initial guess based on p erio dic solution in Eq. 37 with the follo wing steps: 22/35 Initial Guess Based On P erio dic Solution Step 1: Start with light and sleep schedules that follow the p erio dic solution in Eq. 37; Step 2: Change all sleep states during shift work to aw ak e; Step 3: Sim ulate the system forw ard with the up dated sleep schedule and obtain N episo des of sleep so that the i-th episo de { t i | t i sleep ≤ t i < t i wak e } satisfies β ( t i ) = 1, sleep state during shift w ork will remain aw ak e; Step 4: Calculate the up dated H ( t ) with up dated sleep in Step 3, adjust the b egin and end time of sleep episo de i , t i sleep and t i wak e so that t i sleep = min { β ( t i ) == 1 ∧ H + ≤ H ( t i ) ≤ H + max } , and t i sleep = max { β ( t i ) == 1 ∧ H − ≤ H ( t i ) ≤ H − max } , where ∧ represen ts logic ‘and’; Step 5: Increase i by 1; Step 6: Rep eat Step 3 to 5 un til i > N , the num b er of sleep episo des. This initial guess will provide a sleep schedule that satisfies the comfort constraints in Eq. 46. Note that in [13], Yin et al. assume that sleep is only a single contin uous every da y . This excludes napping, a short p eriod of sleep in addition to the usual long sleep at nigh t. Ho w ever, multiple studies hav e shown that taking a late nap b efore the night w ork can improv e the alertness during w ork [40, 41]. W e extend the metho ds in [13] assuming that late nap b egins when H ( t ) = H + and ends when H ( t ) = H − or w ork time is reac hed. Our simulations also demonstrate the same phenomenon for late naps. W e consider three consecutive night shifts that b egin at 11 PM and end at 7 AM of the next morning, and the optimization ob jective is to maximize alertness during three shifts. The optimization horizon b egins at 8 AM on the first day and ends at 7 AM on the fourth day . One optimization uses the initial guess based on the p eriodic solution as w e discussed and is optimized using the gradient descent algorithm. The sub ject will nap for 0.7 hours, 2.35 hours and 2.3 hours b efore the shift work, resp ectively . Another optimization uses the initial guess also based on p eriodic solution, but excludes the naps b efore the night shift. The optimal light and sleep schedules are plotted in Fig. 7. 23/35 When the optimal ligh t is 0 lux, the sub ject needs to wear blo c king goggles to blo c k part of the blue sp ectrum [42]. W e compare the optimized schedule of these tw o optimizations using a verage alertness A av g = J 24 . The first schedule that has naps b efore the shift work has an av erage of A av g = − 1 . 368 during the shift w ork, and the second sc hedule that do es not ha ve naps b efore the shift work has an av erage of A av g = − 1 . 4533 during the w ork shift. The improv emen t in alertness with naps b efore shifts can b e explained by the decrease in H during late napping. Ho wev er, in different scenarios, the optimal schedules will b e different. If we can adjust the circadian rh ythm ahead of the mission by adapting to different light and sleep sc hedules, alertness during shift work will b e higher [12]. W e run a series of sim ulations that start the optimization 1 to 13 days ahead of the first day of the shift w ork. The optimal schedules with 6 and 13 days of preparation are plotted in Fig. 7(c) and (d), resp ectiv ely . As the preparation time increases, the av erage alertness during shift w ork also increases, as shown in Fig. 7(e). Av erage alertness reaches a plateau after 13 days as the circadian rhythms of the sub ject are fully adjusted to align with the shift work at night. Adjusting circadian rhythms provides a more significant increase in alertness, and late napping is not required in this scenario. In summary , if the sub ject cannot adjust the light and sleep schedules multiple days in adv ance, late napping b efore shift work can increase av erage alertness during shift w ork. But if adjusting light and sleep schedules b efore shift w ork, the gradient descent algorithm can find b etter schedules that change the circadian rh ythms of the sub ject to increase the a verage alertness during shift work. 3.2 Cum ulativ e Alertness Optimization In this subsection, w e solve the cumulativ e alertness optimization problem: Giv en the initial conditions X (0) = X 0 , and the dynamics of the h ybrid mo del represen ted as ˙ X = D i ( X, I , β ) for i = 1 , · · · , N , light constraint Eq. 44 and tunable sc hedule, we maximize the cumulativ e alertness of a p erson during certain p erio d: max Z t f t 0 (1 − β ( t ))[ H + ( t ) − H ( t )] dt, (60) where t 0 is the starting time and t f is the end time of the optimization horizon. 24/35 (a) Optimal Sc hedules With Naps (b) Optimal Schedules Without Naps (c) Optimal Sc hedules With 6 Da ys Preparation (d) Optimal Schedules With 13 Da ys Preparation (e) Average Alertness During Shift W ork Fig 7. Optimal schedules for 3 consecutive night shifts. The shift work time is indicated by the blue curve. The alertness A is represented by the red curve. The green region indicates sleep time, and red region indicates when the light is turned on. (a) and (b) represen t the light and sleep schedules which the optimization b egins on the first da y of shift work. (c) represen ts the light and sleep schedules that start 6 days ahead of the shift w ork. (d) represents the light and sleep schedules that start 13 days ahead of the shift w ork. (e) represents the num b er of hours the optimization b egins b efore the first shift work versus the av erage alertness during the work shift. 25/35 The gradien t descent algorithm in Section 3.1 applies to this problem with the only exception of L in Eq. 49 is replaced by L = − (1 − β ( t ))[ H + ( t ) − H ( t )] . (61) Here, w e study tw o cases of cum ulative alertness optimization. Both cases start at 12 PM and end at 12 AM on the fifth night. In the first case, w e optimize the cumulativ e alertness of a sub ject that accommo dates the p erio dic light in Eq. 12 and sp ontaneous sleep schedule. W e find that the cum ulative alertness increases from 89.28 to 89.56 after optimization. The alertness comparison is plotted in Fig. 8(a). The difference is less than 0.3%. This sho ws that the p erio dic light schedule and sp ontaneous sleep schedule are very close to the ligh t and sleep schedules that maximize the cumulativ e alertness during the optimization horizon if there is no shift w ork. The second case has night shifts that start at 10 PM and end at 6 AM in the first three nights. After nigh t shifts, the sub ject falls asleep sp on taneously . W e start with the initial guess based on p erio dic solution. The cumulativ e alertness is − 2 . 3247 b efore optimization. After optimizing the light and sleep sc hedules, cumulativ e alertness increases to 27 . 631. The alertness comparison is plotted in Fig. 8(b). In Fig. 8(c), the circadian phase on day 4 is delay ed 3.1 hours compared to da y 1. This shows that optimal light and sleep schedules delay the circadian phase to align with activities at nigh t. In [23], Hong et al. use the PR mo del to form ulate circadian necessary sleep(CNS), whic h is the minimum sleep duration that sleep homeostasis H decreases b elow the sleep threshold H + . Natural wak e time in terv al (NWTI) is defined as the time such that { t | H ( t ) ≤ H + ( t ) } so that the sub ject is naturally aw ak e. If the sub ject’s sleep time is less than CNS, then the sleep is circadian insufficient sleep (CIS) and the sub ject will exp erience more sleepiness during the wak e-up time. The sub ject following CNS sho ws more NWTI than the sub ject follo wing CIS. W e hypothesize that optimizing cumulativ e alertness will not reduce NWTI. W e run sim ulations with our hybrid PR mo del and compare NWTI with the CNS sc hedule. Fig. 9 shows a comparison of NWTI b et w een the sub ject follo wing the CNS schedule and the sub ject following the optimal cumulativ e alertness schedule obtained from the gradient 26/35 (a) Sub ject follows p erio dic schedule (b) Sub ject follows 3 night shifts (c) Phase Comparison of the Sub ject follows 3 night shifts Fig 8. Alertness comparisons of the sub ject b efore and after optimization. (a) represen ts the case of a sub ject accommo dates to the p erio dic light and (b) represents the case that has 10 PM to 6 AM night shifts. A ref ( t ) represen ts the alertness of the sub ject b efore optimization. A ( t ) represents the alertness of the sub ject after the optimization. The green region indicate sleep time and red region indicates when the ligh t is turned on. (c) Circadian phase comparison of the reference schedule b efore optimization in red and the sc hedule after optimization in blue. descen t algorithm. The NWTI increases from 31 . 91 hours to 32 . 81 hours. W e randomly sampled 3-nigh t consecutive schedules in which every shift work b egins b etw een 3 PM 27/35 and 11:00 PM and the shift duration is randomly sampled b et ween 4 and 12 hours. The b eginning and ending time of the optimization horizon is 12 PM on the first day and 12 Am on the fifth da y . In 50 simulations, the av erage increase in cum ulative alertness is 29.97. In 42 cases, the natural wak e time in terv al was longer under the optimal cum ulative alertness schedule than under the CNS schedule. The av erage p ercen tage NWTI increase is 1.24% with a standard deviation of 1.61%. Therefore, the optimal sc hedules obtained from gradient descent algorithm can increase alertness during the full optimization horizon and not sacrifice NWTI. (a) Sub ject follows CNS schedule (b) Sub ject follows optimal cumulativ e alertness schedule Fig 9. Natural wak e time interv al comparisons of the sub ject follo wing CNS schedule and optimal cum ulative alertness schedule. This shows the plot of alertness of the sub ject b efore and after the optimization. A ref ( t ) represents the alertness of the sub ject b efore optimization. A ( t ) represen ts the alertness of the sub ject after the optimization. Green regions indicate the natural wak e time interv al and red regions indicate the shift w ork. 28/35 4 Conclusion This pap er prop oses the h ybrid PR mo del that consists of the ph ysiology underlying sleep-w ake dynamics and can accurately predict cognitive abilities b etter than the three-pro cess mo del. The hybrid mo del do es not hav e stiff equations in the full PR mo del and has contin uous gradients. W e apply the calculus of v ariation to find the light and sleep sc hedule that maximizes cumulativ e alertness during shift w ork and the full optimization horizon. W e find that if we can adjust the circadian rhythms b efore shift w ork, alertness during shift will b e higher than not adjusting the circadian rhythm. If we cannot adjust the circadian rh ythms b eforehand, taking a late nap b efore the shift will increase alertness during shift work. W e also compare the optimal cumulativ e alertness sc hedule with the CNS schedule and demonstrate that maximizing cum ulative alertness during the full optimization horizon do es not sacrifice the natural wak e time in terv al. 5 App endix 5.1 Calculus of V ariation Deriv ation Supp ose the dynamics for the i -th mo de is describ ed by the following dynamics: ˙ X = D i ( X, I ) , t ∈ [ t i − 1 , t i ) , ∀ i ∈ { 1 , 2 , · · · , N } , (62) where N is the total num ber of mo des during the optimization time horizon. t 0 = 0 and t N = t f are fixed. All the switching time t i are sub ject to the inequality constraints in tunable sc hedule in Eq. 46. An ob jective function is formulated as J = Z t f 0 L ( t, X, I ) dt. (63) W e in tro duce Lagrange m ultipliers λ ( t ) ∈ R 4 . The augmented ob jectiv e function is: J a = Z t f 0 L ( τ , X , I ) dτ + N X i =1 Z t i t i − 1 λ ( τ )[ D i ( X, I ) − ˙ X ( τ )] dτ . (64) 29/35 W e in tro duce a small scalar α and add p erturbations to the inputs αδ I and α δt i . The p erturbed augmented ob jectiv e function b ecomes: J a + αδ J a + o ( α ) = Z t f 0 L ( τ , X + αδ X + o ( α ) , I + αδ I ) dτ + N X i =1 Z t i + αδ t i t i − 1 + αδ t i − 1 λ T ( τ )[ D i ( X + α δX + o ( α ) , I + αδ I − ˙ X − α δ ˙ X − o ( α )] dτ , (65) where δ X, δ J a are the first order p erturbation in states and augmented cost, and o ( α ) represen ts the terms that has higher order than first order of α . F ormally , δ X = lim α → 0 X ( I ( t ) + αδ I ( t ) , t i + αδ t i ) − X ( I ( t ) , t i ) α , (66) δ J a = lim α → 0 J a ( I ( t ) + αδ I ( t ) , t i + αδ t i ) − J a ( I ( t ) , t i ) α . (67) The first v ariation of the augmented cost is calculated by subtracting the p erturbed ob jective function Eq. 65 by Eq. 64 and taking the limit α → 0. δ J a = Z t N t 0 (  δ L ( τ , X , I ) δ X  T δ X + δ L ( τ , X , I ) δ I δ I ) dτ + N X i =0 Z t i t i − 1 ( ˙ λ T δ X + λ T ( τ ) "  δ D i ( X, i ) δ X  T δ X + δ D i ( X, I ) δ I δ I #) dτ − N X i =1  λ T ( t − i ) δ X ( t − i ) − λ T ( t + i − 1 ) δ X ( t + i − 1 )  + N − 1 X i =1  λ T ( t − i ) D i ( X ( t − i , I ( t − i ))) − λ T ( t + i ) F i +1 ( X ( t + i ) , I ( t + i ))  δ τ i + N − 1 X i =1  L ( t − i , X ( t − i ) , I ( t − i )) − L ( t + i , X ( t + i ) , I ( t + i ))  δ τ i (68) By selecting the Lagrange m ultiplier λ ( t ) as: λ ( t N ) = 0 , (69) ˙ λ ( t ) = − δ L ( t, X, I ) δ X −  δ D i ( X, I ) δ X  T λ ( t ), when t ∈ [ t i − 1 , t i ) , (70) λ ( t − i ) = λ ( t + i ) , (71) 30/35 δ J a can b e simplified to δ J a = N X i =1 Z t i t i − 1  δ L ( τ , X , I ) δ I + λ T ( τ ) δ D i ( X, I ) δ I  T δ I dτ + N − 1 X i =1 [ L ( t − i , X ( t − i ) , I ( t − i )) − L ( t + i , X ( t + i ) , I ( t + i )) + λ ( t − i ) D i ( X ( t − i ) , I ( t − i )) − λ T ( t + i ) D i +1 ( X ( t + i ) , I ( t + i ))] δ t i . (72) F rom Eq. 72, the gradient of J with resp ect to I and t i are giv en as ∇ I ( t ) J = ∂ L ( t, X , I ) ∂ I + λ T ( t ) ∂ D i ( X, I ) ∂ I , (73) ∇ t i J = L ( t − i , X ( t − i ) , I ( t − i )) − L ( t + i , X ( t + i ) , I ( t + i )) + λ ( t − i ) D i ( X ( t − i ) , I ( t − i )) − λ T ( t + i ) D i +1 ( X ( t + i ) , I ( t + i )) (74) References 1. Iliff JJ, W ang M, Liao Y, Plogg BA, Peng W, Gundersen GA, et al. A para v ascular path wa y facilitates CSF flow through the brain parench yma and the clearance of in terstitial solutes, including amyloid β . Science translational medicine. 2012;4(147):147ra111–147ra111. 2. Iliff JJ, Lee H, Y u M, F eng T, Logan J, Nedergaard M, et al. Brain-wide pathw a y for w aste clearance captured by contrast-enhanced MRI. The Journal of clinical in vestigation. 2013;123(3):1299–1309. 3. Born J, Wilhelm I. System consolidation of memory during sleep. Psyc hological researc h. 2012;76:192–203. 4. Slater G, Steier J. Excessiv e da ytime sleepiness in sleep disorders. Journal of thoracic disease. 2012;4(6):608. 5. James SM, Honn KA, Gaddameedhi S, V an Dongen HP . Shift work: disrupted circadian rh ythms and sleep—implications for health and well-being. Curren t sleep medicine rep orts. 2017;3:104–112. 31/35 6. Phillips A, Robinson P A. A quantitativ e mo del of sleep-wak e dynamics based on the ph ysiology of the brainstem ascending arousal system. Journal of biological rh ythms. 2007;22(2):167–179. 7. Boivin DB, T rembla y GM, James F O. W orking on a t ypical schedules. Sleep medicine. 2007;8(6):578–589. 8. Do yle I II FJ, Gunaw an R, Bagheri N, Mirsky H, T o TL. Circadian rhythm: A natural, robust, m ulti-scale control system. Computers & c hemical engineering. 2006;30(10-12):1700–1711. 9. Serkh K, F orger DB. Optimal schedules of light exp osure for rapidly correcting circadian misalignmen t. PLoS computational biology . 2014;10(4):e1003523. 10. Ab el JH, Chakrabarty A, Doyle FJ. Con trolling biological time: Nonlinear mo del predictiv e control for p opulations of circadian oscillators. Emerging Applications of Con trol and Systems Theory: A F estschrift in Honor of Mathukumalli Vidy asagar. 2018; p. 123–138. 11. Booth V, Xique I, Diniz Behn CG. One-dimensional map for the circadian mo dulation of sleep in a sleep-wak e regulatory netw ork mo del for human sleep. SIAM Journal on Applied Dynamical Systems. 2017;16(2):1089–1112. 12. Zhang J, W en JT, Julius A. Optimal circadian rhythm control with light input for rapid en trainment and improv ed vigilance. In: 2012 IEEE 51st IEEE Conference on Decision and Con trol (CDC). IEEE; 2012. p. 3007–3012. 13. Yin J, Julius A, W en JT, W ang Z, He C, Kou L. Human Alertness Optimization with a Three-Pro cess Dynamic Mo del. Mathematics. 2022;10(11):1916. 14. Borb ´ ely AA, et al. A tw o pro cess mo del of sleep regulation. Hum neurobiol. 1982;1(3):195–204. 15. Ac hermann P , Borb ´ ely AA. Simulation of daytime vigilance by the additive in teraction of a homeostatic and a circadian pro cess. Biological cyb ernetics. 1994;71(2):115–121. 32/35 16. Ac hermann P , W erth E, Dijk DJ, Borb ely A. Time course of sleep inertia after nigh ttime and daytime sleep episo des. Archiv es italiennes de biologie. 1996;134(1):109–119. 17. Reifman J, Gander P . Commen tary on the three-pro cess mo del of alertness and broader mo deling issues. Aviation, space, and en vironmental medicine. 2004;75(3):A84–A88. 18. P ettersson K, M ¨ uller K, Tiet¨ av¨ ainen A, Gould K, Hæggstr¨ om E. Saccadic eye mo vemen ts estimate prolonged time aw ake. Journal of Sleep Research. 2019;28(2):e12755. 19. Diniz Behn C, Kop ell N, Brown EN, Mo chizuki T, Scammell TE. Dela yed orexin signaling consolidates wak efulness and sleep: ph ysiology and mo deling. Journal of neuroph ysiology . 2008;99(6):3090–3103. 20. Rempe MJ, Best J, T erman D. A mathematical mo del of the sleep/wak e cycle. Journal of mathematical biology . 2010;60(5):615–644. 21. Phillips AJ, Chen P , Robinson P . Probing the mechanisms of chronot yp e using quan titative mo deling. Journal of Biological Rhythms. 2010;25(3):217–227. 22. P apatsimpa C, Linnartz JP . P ersonalized office lighting for circadian health and impro ved sleep. Sensors. 2020;20(16):4569. 23. Hong J, Choi SJ, Park SH, Hong H, Bo oth V, Jo o EY, et al. P ersonalized sleep-w ake patterns aligned with circadian rhythm relieve daytime sleepiness. Iscience. 2021;24(10). 24. Khalil HK, Grizzle JW. Nonlinear systems. vol. 3. Pren tice hall Upp er Saddle Riv er, NJ; 2002. 25. Jew ett ME, F orger DB, Kronauer RE. Revised limit cycle oscillator mo del of h uman circadian pacemaker. Journal of biological rhythms. 1999;14(6):493–500. 26. Borb AA, Achermann P . Sleep homeostas is and mo dels of sleep regulation. Journal of biological rh ythms. 1999;14(6):559–570. 33/35 27. F orger DB, Jewett ME, Kronauer RE. A simpler mo del of the human circadian pacemak er. Journal of biological rhythms. 1999;14(6):533–538. 28. P ostnov a S, Lo c kley SW, Robinson P A. Prediction of cognitiv e p erformance and sub jective sleepiness using a mo del of arousal dynamics. Journal of biological rh ythms. 2018;33(2):203–218. 29. Song YM, Choi SJ, P ark SH, Lee SJ, Jo o EY, Kim JK. A real-time, p ersonalized sleep in terven tion using mathematical mo deling and wearable devices. Sleep. 2023;46(9):zsad179. 30. F ulcher B, Phillips AJ, Robinson P A. Modeling the impact of impulsive stimuli on sleep-w ake dynamics. Physical Review E—Statistical, Nonlinear, and Soft Matter Ph ysics. 2008;78(5):051920. 31. F ulcher B, Phillips AJ, Robinson P A. Quan titative physiologically based mo deling of sub jectiv e fatigue during sleep depriv ation. Journal of theoretical biology . 2010;264(2):407–419. 32. P ostno v a S, Postno v DD, Seneviratne M, Robinson P A. Effects of rotation in terv al on sleepiness and circadian dynamics on forward rotating 3-shift systems. Journal of biological rh ythms. 2014;29(1):60–70. 33. Johnson M, Duffy J, Dijk D, Ronda J, Dyal C, Czeisler C. Short-term memory , alertness and p erformance: a reappraisal of their relationship to b ody temp erature. Journal of sleep research. 1992;1(1):24–29. 34. W ehr T A. In short photop erio ds, human sleep is biphasic. Journal of Sleep Researc h. 1992;1(2):103–107. 35. Monk TH. Sub jective ratings of sleepiness—the underlying circadian mechanisms. Sleep. 1987;10(4):343–353. 36. Shahid A, Wilkinson K, Marcu S, Shapiro CM. Stanford sleepiness scale (SSS). In: STOP , THA T and one hundred other sleep scales. Springer; 2011. p. 369–370. 37. ˚ Ak erstedt T. Is there an optimal sleep-wak e pattern in shift work? Scandinavian Journal of W ork, En vironment & Health. 1998; p. 18–27. 34/35 38. Vital-Lopez FG, Doty TJ, Reifman J. Optimal sleep and work schedules to maximize alertness. Sleep. 2021;44(11):zsab144. 39. Yin J, Julius AA, W en JT. Optimization of light exp osure and sleep schedule for circadian rh ythm entrainmen t. Plos one. 2021;16(6):e0251478. 40. Drisk ell JE, Mullen B. The efficacy of naps as a fatigue countermeasure: a meta-analytic in tegration. Human F actors. 2005;47(2):360–377. 41. Sc h weitzer PK, Randazzo AC, Stone K, Erman M, W alsh JK. Lab oratory and field studies of naps and caffeine as practical countermeasures for sleep-w ake problems asso ciated with night work. Sleep. 2006;29(1):39–50. 42. Rahman SA, Marcu S, Shapiro CM, Brown TJ, Casp er RF. Spectral mo dulation atten uates molecular, endo crine, and neurob eha vioral disruption induced b y no cturnal light exp osure. American Journal of Physiology-Endocrinology and Metab olism. 2011;300(3):E518–E527. 35/35