Application of a Limit-Cycle Oscillator Model for Prediction of Circadian Phase in Rotating Night Shift Workers

Practical alternatives to gold-standard measures of circadian timing in shift workers are needed. We assessed the feasibility of applying a limit-cycle oscillator model of the human circadian pacemaker to estimate circadian phase in 25 nursing and medical staff in a field setting during a transition from day/evening shifts (diurnal schedule) to 3–5 consecutive night shifts (night schedule). Ambulatory measurements of light and activity recorded with wrist actigraphs were used as inputs into the model. Model estimations were compared to urinary 6-sulphatoxymelatonin (aMT6s) acrophase measured on the diurnal schedule and last consecutive night shift. The model predicted aMT6s acrophase with an absolute mean error of 0.69 h on the diurnal schedule (SD = 0.94 h, 80% within ±1 hour), and 0.95 h on the night schedule (SD = 1.24 h, 68% within ±1 hour). The aMT6s phase shift from diurnal to night schedule was predicted to within ±1 hour in 56% of individuals. Our findings indicate the model can be generalized to a shift work setting, although prediction of inter-individual variability in circadian phase shift during night shifts was limited. This study provides the basis for further adaptation and validation of models for predicting circadian phase in rotating shift workers.

Working and sleeping out of phase with the endogenous circadian pacemaker is a major factor contributing to adverse health and safety outcomes 1,2 . Biological adaptation of the circadian clock to a night shift schedule, marked by a shift in the timing of melatonin secretion, into the daytime rest period, can improve sleep quality and reduce performance impairments [3][4][5] . Even partial circadian adaptation is, however, rarely observed in most real-world settings 6 .
Multiple field studies have reported large variability between individuals in the direction and magnitude of circadian phase shift in response to night shift work [7][8][9][10] . This inter-individual variability has been linked to differences in the timing of light exposure 7,9,11 , which is the key timing signal for the circadian clock [12][13][14] . Night workers with delayed circadian phase after consecutive night shifts are exposed to more light during the evening compared to those with advanced phase 7 . Inter-individual variability in the magnitude and direction of circadian phase response to night shift can be largely explained by the relative difference in the amount of light exposure in the delay and advance portions of the phase response curve 11 .
Individual differences in circadian response to night work translate to variation between workers in the time course of alertness and performance across the night shift 5,10,[15][16][17] . Knowledge of individual circadian phase in shift workers could identify times of impaired alertness and thereby inform individualized countermeasures for improving workplace safety, overall health, and wellbeing. Direct assessment of circadian phase via measurement of the melatonin rhythm is impractical and prohibitively expensive for ongoing use in occupational settings.

Results
Characteristics of the 25 participating doctors and nurses are summarized in Table 1. Activity and light data, measured using wrist-worn actigraphs, were available for 7.16 ± 3.01 (M ± SD; range 2-14) consecutive 24-h intervals over a work schedule consisting of day/evening shifts and days off (the diurnal schedule), and 3.60 ± 1.08 (M ± SD; range 3-5) consecutive 24-h intervals over a schedule consisting of 3-5 consecutive night shifts (the night schedule). The number of days of actigraphy data (light and activity counts) available varied between participants depending on when they were enrolled into the study relative to their night schedule.
Model estimations of circadian phase. Photic only model. Using only the photic input, B(t), the model estimated an average aMT6s acrophase for all participants on the diurnal schedule of 3:55 ± 0:43 h (range: 2:58-5:39 h). After 3 to 5 consecutive night shifts the photic model predicted an average aMT6s acrophase of 4:34 ± 0:46 h (range: 3:16-6:03 h). The average predicted phase shift was a delay of 0:39 ± 0:42 h (range 2:24 h phase delay to 0:33 h phase advance). The photic-only model accurately predicted aMT6s acrophase to within ±30 minutes in 52% and 32% of individuals, and to within ±1 hour in 64% and 56% on the diurnal and night schedules, respectively ( Table 2). The direction of phase shift was predicted accurately in 80% of participants (12% incorrect prediction of phase advance; 8% incorrect prediction of phase delay). The model predicted a total of 5 phase advances, of which only 2 accurately reflected the direction of aMT6s phase shift. Model phase estimates generated using the photic model had a significant positive relationship with measured aMT6s acrophase on the diurnal schedule (r = 0.63, r 2 = 0.39, p = 0.001; Fig. 2A), and on final night shift (r = 0.60, r 2 = 0.33, p = 0.001; Fig. 2C), and with the degree of phase shift (r = 0.48, r 2 = 0.23, p = 0.015; Fig. 2E). Thus, 39% and 33% of the variance in aMT6s acrophase was explained by the photic model on the diurnal schedule and final night shift, respectively, and 23% of the variance in aMT6s phase shift. At the group level, predictions of phase on the diurnal schedule were not significantly different from aMT6s acrophase times (t(24) = 0.30, p = 0.768). Predicted phase on final night shift, however, was significantly different to measured aMT6s acrophase times (t(24) = 2.52, p = 0.019), indicating an under-prediction of the delay in aMT6s acrophase times. Correspondingly, predicted phase shift differed significantly from measured shift in aMT6s acrophase over the night shift schedule (t(24) = −2.15, p = 0.04), indicating an under-prediction of the magnitude of phase shifts, particularly when there were phase shifts of 2 hours or more.  Acrophase times on third, fourth or fifth consecutive night shift for n = 25 participants. (C) Phase shift from diurnal schedule to final night shift for n = 25 participants. Positive phase shift indicates a phase advance (i.e., final night shift phase occurring at an earlier clock time than diurnal phase); negative phase shift indicates a phase delay (i.e., final night shift phase occurring at a later clock time than diurnal phase). Circles represent measured aMT6s acrophase time (reference phase); squares represent predicted phase using the photic only mode; diamonds represent predicted phase using the photic and non-photic (PNP) model. in 56% and 28% of individuals, and to within ±1 hour in 80% and 68% of individuals on the diurnal and night schedules, respectively ( Table 2). The predicted phase shift based on PNP model predictions was an average phase delay of 1:40 ± 1:01 h (range: 3:40 h phase delay to 0:03 h phase advance). The PNP model correctly predicted the direction of phase shift in 88% of participants (12% incorrect prediction of phase delay). Only one phase advance was predicted, which accurately reflected the direction of the measured aMT6s phase shift for that participant.
Prediction error relative to measured aMT6s acrophase. The prediction error for each model is summarized in Table 2. Model performance using the photic-only and PNP models were comparable when predicting diurnal phase. On the night shift schedule the PNP model showed slightly improved performance compared to the photic-only model, across all metrics of prediction error. There were no significant differences in mean absolute prediction error between photic and PNP models on either diurnal (p = 0.66) or night shift (p = 0.23) schedules, or in phase shift (p = 0.21). When comparing raw error values mean prediction error was smaller when predicting diurnal phase using the photic only model compared to the PNP model (t(24) = −3.678, p = 0.001; Fig. 3). Conversely, mean prediction error was smaller when using the PNP model to predict phase on final night shift (t(24) = 5.757, p < 0.001), or to predict phase shift between the shift schedules (t(24) = −7.889, p < 0.001). There was an offset in mean prediction error (raw error values) using the photic model to predict acrophase on the night schedule such that predicted phase was systematically earlier than aMT6s acrophase.
For both models, there were significant relationships between prediction error and (1) the measured diurnal acrophase, (2) final night shift acrophase, and (3) phase shift (Fig. 4). This relationship was strongest when using the photic model to estimate acrophase on final night shift (r 2 = 0.81, p < 0.001; Fig. 4C) or phase shift from day to final night shift (r 2 = 0.85, p < 0.001; Fig. 4E). There were no significant relationships between prediction error (for predicted diurnal or night phase, or predicted phase shift) and age, BMI, average mid-sleep time on the diurnal schedule, or MEQ score.
Removal of outlier on diurnal schedule. One participant (81) had a relatively late aMT6s acrophase time on the diurnal schedule (6:45 am; more than 2 SD from the mean of 3:58 am ± 1:06 h), with large prediction error when using both photic and PNP models. To examine a potential bias caused by this participant, analyses were also conducted with this participant removed. This resulted in a smaller average prediction error for the photic model (absolute mean = 0.61 ± 0.49 h), but did not improve the strength of the relationship between predicted phase and measured aMT6s acrophase (r 2 = 0.34, p = 0.003). There was also a reduction in prediction error for the PNP model (absolute mean = 0.60 ± 0.54 h), with an improved relationship between predicted and measured phase (r 2 = 0.33, p = 0.003).  www.nature.com/scientificreports www.nature.com/scientificreports/ Impact of initial values on model estimations. To test whether model performance for predicting final night shift phase could be improved with the application of individual starting acrophase, the models were initialized using diurnal aMT6s acrophase (the "oracle" method). Initializing with the known initial acrophase made only a small improvement in the predictions of acrophase on final night shift, and in the predictions of phase shift (see Table S1). When applying the "oracle" model initialization using the PNP model, the mean absolute prediction error in the night schedule reduced from 0.95 h to 0.94 h, and the mean absolute prediction error in the phase shift reduced from 1.11 h to 1.04 h. Model predictions using the PNP model improved from accounting for 52% to 53% of variance in aMT6s night shift acrophase; and from 42% to 49% of variance in phase shift (Fig. S2).

Discussion
This study tested the capability of a mathematical limit-cycle oscillator model to predict circadian timing in the context of rotating night shift work, using easily collected ambulatory light and activity data. The main findings are: (1) the photic and PNP models both performed well on diurnal schedules, with absolute mean prediction errors <0.7 h, and >50% of aMT6s acrophases predicted to within ±30 minutes; (2) the addition of a non-photic input (i.e., incorporating the effects of sleep/wake patterns on circadian phase) improved prediction of aMT6s acrophase after night shifts from an absolute mean error of 1.19 h to 0.95 h; (3) predicted phase shifts were of smaller magnitude than the measured aMT6s phase shifts; (4) knowledge of aMT6s acrophase prior to night shifts only modestly improved model phase estimations compared to initializing with a starting point derived from mid-sleep time.
The photic and PNP models that we tested here have been extensively developed from measurements of human circadian responses under laboratory conditions in highly-screened, healthy individuals 18,19,25 . Although these models have been used to assess real-world schedules, including rotating night shift schedules 30 , prediction www.nature.com/scientificreports www.nature.com/scientificreports/ accuracy under such conditions had not been previously assessed. Recent field studies of individuals living on diurnal schedules demonstrated that these models can predict timing of melatonin phase markers to within ±1 hour in about half of individuals using wrist-worn light and activity data [26][27][28] . We replicated these findings and for the first time evaluated prediction accuracy on night schedules, including considerable variations in individual shift rosters. We conclude that the PNP model is a suitable tool for predicting phase on night schedules, albeit with slightly less accuracy than on diurnal schedules. The importance of including non-photic inputs is an important finding for operational use, and likely reflects the fact that abnormal phase angles between sleep timing and circadian timing often arise on rotating shift schedules. Specifically, we observed improved prediction of later circadian phases, potentially due to the impact of the delayed sleep-wake timing during the night schedule on model accuracy.
Whenever predicting circadian phase, an important but rarely addressed question is: what level of accuracy is required? For example, is a 1 hour error tolerable? The answer will naturally depend upon the outcome being predicted. A key application of predictive models in shift work settings is determining times of likely poor cognitive performance (e.g., 30 ). We can estimate the consequences of different circadian phase prediction errors from laboratory data. Using median reaction time on a psychomotor vigilance task (PVT) under forced desynchrony conditions 31 , we calculated the absolute errors in reaction time and their effect sizes for a 1 or 2 hour error in phase (see Fig. 5). This analysis indicates that effect size varies by circadian phase (i.e., the consequences of a mis-estimate are greater at some circadian phases than others), and that a 1 hour error ensures a medium or smaller effect size, whereas a 2 hour error can result in large effect sizes. These findings suggest that, at least for a cognitive performance outcome, the mean absolute error of 0.95 h obtained for the PNP model under night shift conditions is tolerable. We note, however, that other, non-circadian factors could contribute to prediction errors for cognitive performance, such as the sleep homeostatic process.
Our findings identify areas of potential improvement for the models. First, the models did not replicate the amount of inter-individual variability in circadian phase observed, either on diurnal or night shift schedules. This discrepancy is indicated by the regression lines for predicted versus actual phase which have slopes consistently <1, similar to previous findings 28 . This is unsurprising given that the model uses the same parameter values for all individuals, whereas we know there are individual differences in circadian physiology, including intrinsic circadian period 12,32 and individual light sensitivity 33,[34][35][36][37] . Tuning the τ c parameter for males and females based on published sex differences in circadian period 12 did not improve model performance in this sample. This issue could potentially be addressed by the development of tools to individualize estimation of model parameters [38][39][40] . We did not attempt to perform model parameter individualization on this dataset, since we did not have direct measures of τ c or other model parameters. Second, the models systematically underestimated the amount of phase delay that occurred, particularly for individuals with phase shifts of 2 hours or more. This finding is consistent with a previous application of the model, where the model underestimated phase shifts over a simulated night shift protocol 41 .
A possible reason for the underestimation of phase shifts is limited validation of the model in real-world environments with low light levels such as during night shifts. Participants in this study were exposed to lower average intensity light, and spent a smaller proportion of time in bright light whilst working night shifts, compared to on the diurnal schedule. Participants on night shift spent 24% of waking hours at <10 lux and 61% of waking hours at <100 lux. Laboratory experiments have demonstrated a nonlinear relationship between light intensity and circadian phase response, with half the maximal phase resetting observed at ~100 lux 42 . While the light parameters in the model have been refined to reflect these studies 19 , our findings suggest that the model may still underestimate the biological response to the lower intensity light levels experienced on the night shift schedule.
A limitation of the current study is the measurement of light exposure using a wrist-worn actigraph device, which is an imperfect estimate of retinal light exposure 43 . Future work using a sensor worn at eye-level, or better estimates of overall environmental lighting conditions, would provide a more accurate measure and potentially improve model estimations. Validation of models using wrist-worn sensors is important, however, given their current ubiquity and ease of use in field studies. Additionally, it is now recognized that photopic lux is limited in its ability to quantify the effect of light on the circadian system 44 . The circadian pacemaker is particularly sensitive to the phase-resetting effects of short-wavelength light 45,46 ; however, the model's light input is still based on photopic lux 18,19 . Modification of the model to account for other dimensions of light (e.g., melanopic lux 44 ), along with improved light sensors, such as in 47 , would likely provide a more accurate estimate of the circadian phase response to the changed light-dark cycle in shift workers.
In conclusion, this study indicates the potential utility of a mathematical modeling approach to estimating circadian phase in a real-world shift work context, as an alternative to costly and/or burdensome measurements of circadian phase markers. The study also identified important caveats or limitations to the approach, including lower prediction accuracy across night shifts than on diurnal shifts, and failure to capture the full range of inter-individual differences, meaning while the model performs well for most individuals, the model may perform poorly in individuals with extreme early or late phase angles of entrainment. With further model development, there is potential for use of this approach to predict times of worst performance impairment, which would inform both personal and occupational countermeasures. www.nature.com/scientificreports www.nature.com/scientificreports/ Methods Study protocol. Participants were 20 nurses (15 female) and 5 doctors (3 female), aged 33.3 ± 9.2 y (mean ± SD) with a BMI of 23.8 ± 3.7 kg/m 2 , working rotating shifts in an Intensive Care Unit (ICU) at Austin Health, Melbourne, Australia. Doctors worked a consistent rotating roster of 7 consecutive day shifts (08:00-20:30 h) and 7 consecutive night shifts (20:00-08:30 h), with 7 days off in between each block of shifts (i.e., 7 day, 7 off, 7 night, 7 off). Nurses worked variable shift schedules consisting of day (07:00-15:30 h), evening (15:00-21:30 h), and night (20:00-08:30 h) shifts. Participants were recruited via scheduled in-service presentations and targeted recruitment during shifts. Recruitment sought workers with an upcoming roster consisting of at least 4 day or evening shifts or days off (the 'diurnal schedule') followed by a minimum of 3 consecutive night shifts (the 'night schedule').
Study approval. All participants provided written informed consent prior to enrolment into the study, and received financial compensation for their time. The study protocol was approved by the Austin Health and Monash University Human Research Ethics Committees and was conducted in compliance with standards set by the latest revision of the Declaration of Helsinki.
Light exposure and activity measurements. Light (photopic lux) was recorded in one-minute epochs across the entire shift pattern (diurnal schedule followed by night schedule) using a wrist actigraph device (Actiwatch Spectrum or Spectrum Plus, Philips Respironics, Bend, OR, USA; serviced and calibrated prior to data collection by Philips Respironics), worn on the non-dominant wrist. The same device also recorded wrist activity (activity counts/minute). Participants were asked to wear the device at all times, with the light sensor uncovered, except while showering or when the device would interfere with hospital operational requirements.
Questionnaires. Daily work diaries were completed in which participants recorded the start and end times of each shift worked, including breaks. All shift times were confirmed using payroll information provided by hospital administration. Participants also completed an online sleep-health questionnaire, which included demographic information such as age, sex, body mass index (BMI), and shift work history. Diurnal preference was measured using the Horne & Ostberg Composite Morningness Eveningness Questionnaire (MEQ) 48 . www.nature.com/scientificreports www.nature.com/scientificreports/ Reference circadian phase measurements. Circadian phase was assessed by measuring the rhythm of the urinary melatonin metabolite 6-sulphatoxymelatonin (aMT6s) over 24-48 hours on both the diurnal schedule (day/evening shifts or days off) and end of the night schedule (night 3, 4 or 5). Participants collected all urine voids over the collection period, in 4-hourly sequential blocks (8-hourly during sleep). For each block the total urine volume and sampling time were recorded, and a 5-mL sample retained and stored at −20 °C. Samples were analysed for aMT6s concentration using radioimmunoassay 49 at the Adelaide Research Assay Facility (University of Adelaide, Australia) with reagents provided by Stockgrand Ltd (University of Surrey, Guildford, UK). Assay details are reported elsewhere 11 .
Cosinor analysis was conducted to determine the acrophase (peak) time in the aMT6s rhythm 50 at two time points: once during a sequence of confirmed diurnal schedules, and once on the final night shift. All included urine collections passed visual inspection of the aMT6s rhythm, and had a significant cosinor fit (98% had α = 0.10, 2% (n = 2) p < 0.12).

Sleep assessment.
Rest intervals were determined based on activity counts using Actiware software (Actiware 6 software, Philips Respironics, Bend, OR, USA). Main rest intervals were examined for the following criteria: minimum duration of 3 hours, mean activity below a proprietary threshold, and mean light levels below a proprietary threshold. Rest intervals that did not meet these criteria were visually inspected: time boundaries were adjusted at the start and end of each episode where there was a discrepancy between the rest interval and the activity and light distributions of ≥ 60 minutes. Rest intervals shorter than 3 hours that occurred outside the expected sleep window were marked as naps and excluded from the non-photic input into the model. Five rest intervals from four participants were missing due to the participant removing the actigraph overnight. For these cases, where off-wrist intervals aligned approximately with the time and duration of the expected sleep episode, activity and light levels were set to zero, and treated as a main rest interval.
The mid-sleep time was calculated as the mid-point between sleep onset and sleep offset for all main sleep intervals, for use in model initialization (described below). Main rest intervals only were used to calculate the square waveform non-photic input for the model (details below).
Model structure. In this study we used an existing limit-cycle oscillator model 19 implemented with two approaches: (a) using the photic drive B(t) only (photic model, N(t) set to zero), and (b) using both the photic B(t) and non-photic N(t) driving terms (PNP model), see Fig. 6.
The main pacemaker equations are listed below, which describe the modeling framework used to calculate circadian phase estimations.   6. Schematic of the model of the central pacemaker driven by light via B (t) and rest derived from activity via N (t). The model is a limit-cycle oscillator of Van der Pol type in the (X, X ) c plane (see Eqs 1 and 2). The light B(t)) and activity (N(t)) driving terms are obtained from the corresponding external stimuli through specific sensitivity modulators depending on the oscillator state (see Eqs 3 and 4). B (t) is obtained from the light intensity measurements and N (t) from the activity counts. Both are collected by a wrist-worn device with a oneminute sampling rate and result from non-linear processing stages described in 19 . The output of the model consists of time-dependent trajectories in the phase plane, from which daily core body temperature minimum times are estimated.
variable based on a Lienard transformation. This enables the pacemaker to be modeled in a two-dimensional plane. The resulting state vector generates elliptic-like trajectories over time, the shapes of which are modulated by external stimuli. Equations (3) and (4) are circadian modulators, respectively, for the photic B(t) and non-photic N(t) stimuli. In Eq. (3), B t ( ) is the output of a dynamic light processor fed with the measurements of white light intensity (photopic lux) over successive one minute epochs. Similarly, N t ( ) results from the processing of activity counts from wrist actigraphy, enabling the epochs to be tagged in terms of wake or rest status. In the above equations, there are 6 additional model parameters, namely: 1. μ: stiffness parameter, set to 0.13 2. τ c : intrinsic period, set to 24.2 h 3. L q : model parameter, set to 1/3 4. L k : model parameter, set to 0.55 5. l s : light sensitivity parameter, set to 0.4 6. a s : activity sensitivity parameter, set to 10.0 All details relevant to the computation of B t ( ) and N t ( ) from the light and activity measurements have been described previously 19 , including the model parameter values that have been used in the present work. All data pre-processing and phase estimations were implemented in Matlab (R2016b, Mathworks).
The model reports core body temperature minimum (CBTmin) as the circadian phase marker. Based on published relationships between CBTmin and plasma melatonin midpoint 51 , and between plasma melatonin midpoint and urinary aMT6s acrophase [52][53][54] , it was assumed that the timing of aMT6s acrophase in the current sample may be used as a proxy for CBTmin for the purpose of comparing model phase estimations with measured aMT6s acrophase times. initialization procedure. Equations (1) and (2) require initial values {x(t 0 ), x c (t 0 )} at the start time, t 0 , of each recording. To derive these initial values, average mid-sleep time was used as an estimate of the first CBT nadir, which is converted to an initial position in phase space, based on assuming a uniform angular velocity and constant amplitude of the rhythm. Average mid-sleep time was calculated for each individual over the diurnal shift schedule, prior to the subsequent night shift schedule. This approach relies on previously reported relationships between average mid-sleep time and circadian phase (dim light melatonin onset (DLMO)) on a regular diurnal sleep schedule 55,56 , such that mid-sleep time was used as a proxy for initial circadian phase.
To test whether having prior phase measurements could improve model predictions, an alternate initialization procedure was also tested for the night schedules, where the model was initialized using the aMT6s acrophase times measured on the diurnal schedule. We refer to this alternate initialization as the "oracle" method to emphasize that use is made of a reference acrophase time to derive initial parameter values for Eqs (1) and (2), as opposed to the applied heuristic (i.e., average mid-sleep time on diurnal schedule), which does not require a prior phase measurement. Data analysis. Statistical analyses were performed using SPSS version 23 (SPSS Inc., Chicago, IL, USA).
Measured aMT6s phase shift from diurnal schedule to the end of the night schedule was calculated as diurnal acrophase minus night acrophase, such that negative values indicate a phase delay and positive values indicate a phase advance. Predicted acrophase times on the diurnal and night schedules on the same dates that measurements were taken were used to calculate an estimated phase shift (i.e., predicted phase shift). Pearson's correlations were used to examine the relationship between measured aMT6s acrophase and predicted acrophase times on both the diurnal and night schedules, and between measured and predicted phase shifts between diurnal and night schedules. Planned paired samples t-tests were used to compare estimated circadian phase to measured aMT6s phase, and to compare the prediction error between photic and PNP models at each time point. Pearson's correlations were used to examine whether age, MEQ, BMI, and measured aMT6s acrophase were related to individual differences in prediction error on the diurnal schedule or night shift schedule, or to predicted phase shift.