Melissa A St

Hilaire
email: msthilaire@rics.bwh.harvard.edu

Claude Gronfier

Jamie M Zeitzer

Elizabeth B Klerman

A physiologically-based mathematical model of melatonin including ocular light suppression and interactions with the circadian pacemaker

Keywords: melatonin, mathematical models, sleep, circadian, suppression, plasma, saliva, ocular

The rhythm of plasma melatonin concentration is currently the most accurate marker of the endogenous human circadian pacemaker. A number of methods exist to estimate circadian phase and amplitude from the observed melatonin rhythm. However, almost all of these methods are limited because they depend on the shape and amplitude of the melatonin pulse, which varies among individuals and can be affected by environmental influences, especially light. Furthermore, these methods are not based on the underlying known physiology of melatonin secretion and clearance, and therefore cannot accurately quantify changes in secretion and clearance observed under different experimental conditions. A published physiologically-based mathematical model of plasma melatonin can estimate synthesis onset and offset of melatonin under dim light conditions. We amended this model to include the known effect of melatonin suppression by ocular light exposure and to include a new compartment to model salivary melatonin concentration, which is widely used in clinical settings to determine circadian phase. This updated model has been incorporated into an existing mathematical model of the human circadian pacemaker and can be used to simulate experimental protocols under a number of conditions.

Introduction

The activity of the endogenous circadian pacemaker located in the suprachiasmatic nucleus (SCN) within the hypothalamus cannot be measured directly in humans. In order to extract information about circadian phase or amplitude or other characteristics of the circadian system, overt marker rhythms must be used. A comparison of the plasma melatonin rhythm to two other marker rhythms, core body temperature and cortisol, revealed that phase markers measured by melatonin were the least variable under relatively constant conditions, and therefore are more precise markers of the endogenous circadian rhythm [START_REF] Klerman | Comparisons of the variability of three markers of the human circadian pacemaker[END_REF]. Melatonin is also commonly used as a marker in circadian rhythms research because it can be assayed in blood, saliva or urine and because it is relatively unaffected by sleep-wake state or changes in posture [START_REF] Arendt | Melatonin: characteristics, concerns, and prospects[END_REF].

Multiple mathematical methods have been published (for review, see [START_REF] Arendt | Melatonin: characteristics, concerns, and prospects[END_REF] and ( 1)) by which circadian phase can be extracted from the plasma melatonin rhythm. Most of these methods depend on curve-fitting of the melatonin profile and/or the crossing of a threshold to determine phase. Thus, these methods assume certain characteristics about the shape and amplitude of the melatonin profile. For example, to calculate phase as the "quarter crossing of the fit" (DLMO 25% ), the amplitude of a melatonin pulse is calculated from a 3-harmonic fit of the data, which assumes that the shape of the melatonin pulse can be captured by a fundamental plus 2harmonic (e.g., approximately 24, 12, and 8 hr) sinusoidal curve. Then, 25% of this fit amplitude is calculated and linear interpolation is used to estimate the time at which the plasma melatonin concentration crosses this value. This time is used as the phase marker (see Fig. 3

of (3)).

This or similar phase estimates may not be accurate for several reasons: a) There is large intraindividual variability in melatonin amplitude and shape across several days, especially during an experimental intervention such as a light pulse under different background lighting conditions (4); b) There is large inter-individual variability in melatonin amplitude [START_REF] Van Coevorden | Neuroendocrine rhythms and sleep in aging men[END_REF], [START_REF] Zeitzer | Do plasma melatonin concentrations decline with age?[END_REF]. Therefore, some threshold methods may not be accurate in individuals with low amplitude of melatonin concentration; c) Linear interpolation assumes that the change in melatonin concentration across the interpolation interval is linear, which may not be true, especially when sampling rate is infrequent.

There are two other problems with most current measures of melatonin phase and amplitude.

Most investigators agree that the increase of melatonin secretion is a marker of circadian phase [need reference]. However, the different definitions of this increase -which depend on the threshold or fit used -make it difficult to compare results across studies, or even to associate the time of the change with the physiology. A second problem is that the error in the calculation is not reported. Some circadian rhythms experiments report phase shifts of as little as 30 minutes as significant. However, without knowledge of the error in the phase marker, the statistical significance of the results and the statistical power of the experiments may not be known.

The suppression of melatonin by ocular light exposure has been used to probe the sensitivity of the circadian pacemaker to light in blind persons [START_REF] Klerman | Photic resetting of the human circadian pacemaker in the absence of conscious vision[END_REF], to varying intensities and durations of polychromatic light exposure, and to long-term changes in light exposure [START_REF] Smith | Adaptation of human pineal melatonin suppression by recent photic history[END_REF]. This method is appropriate because the pathway for ocular light suppression of melatonin begins in the retina, travels via the single-synaptic retino-hypothalamic tract (RHT) to the SCN of the hypothalamus, through the paraventricular nuclei (PVN) and the superior cervical ganglion (SCG) and ends in the pineal [START_REF] Klein | Pineal N-acetyltransferase and hydroxyindole-Omethyltransferase: Control by the retinohypothalamic tract and the suprachiasmatic nucleus[END_REF]. Intensity response curves [START_REF] Zeitzer | Sensitivity of the human circadian pacemaker to nocturnal light: Melatonin phase resetting and suppression[END_REF] demonstrate that the half-maximum sensitivity to light estimated from plasma melatonin suppression occurs at approximately 100 lux, with maximum melatonin suppression occurring between 1000 and 10000 lux. For comparison, ordinary room light levels are approximately 100 lux; studies conducted in room light may inadvertently include melatonin suppression in addition to other specified experimental effects.

More recent studies have also reported melatonin suppression by monochromatic light with peak sensitivity in the short-wavelength bands [START_REF] Brainard | Action spectrum for melatonin regulation in humans: Evidence for a novel circadian photoreceptor[END_REF] [START_REF] Thapan | An action spectrum for melatonin suppression: Evidence for a novel non-rod, non-cone photoreceptor system in humans[END_REF] [START_REF] Wright | Effect of light wavelength on suppression and phase delay of the melatonin rhythm[END_REF] [START_REF] Cajochen | High sensitivity of human melatonin, alertness, thermoregulation, and heart rate to short wavelength light[END_REF].

In 1997, Brown et al. [START_REF] Brown | A mathematical model of diurnal variations in human plasma melatonin levels[END_REF] introduced a method for calculating melatonin phase and amplitude using a physiologically-based differential equation model of the diurnal variation of plasma melatonin concentrations. Details of the model are included below, in the Methods section. This model can be used to estimate the onset (SynOn) and offset (SynOff) of plasma melatonin production, as well as the amount of melatonin secreted (amplitude). The onsets and offsets estimated by the model were comparable to assessments of melatonin phase based on threshold methods. This method also yields information about the infusion rates of melatonin into and clearance rates of melatonin out of the plasma, information that cannot be extracted using the other methods but that may be useful in understanding the physiology. A method to estimate melatonin phase that is based on physiology may be a much more useful tool for melatonin analysis because there are no intrinsic assumptions about the shape and amplitude of the melatonin pulse and because the parameters (e.g., synthesis onset and offset, infusion and clearance rates, amplitude) can be directly related to secretory processes. Such information may provide physiological insight into the differences in melatonin levels seen across individuals, and in different conditions. Salivary melatonin rhythms have also been used as markers of the endogenous circadian pacemaker. Saliva samples for melatonin assay are often collected when it is impossible or impractical to collect plasma samples. Phase estimates are usually calculated by the same curvefitting and threshold methods used for plasma melatonin. However, in addition to the problems detailed above for analyzing plasma samples for circadian information, there are at least two additional problems with using saliva for this purpose: (a) Salivary melatonin concentrations are ~30% lower than plasma melatonin concentrations [START_REF] Voultsios | Salivary melatonin as a circadian phase marker: Validation and comparison to plasma melatonin[END_REF] which may make threshold methods even more difficult to apply and (b) often only the rising portion of the salivary melatonin pulse is collected, because most melatonin secretion occurs when an individual is asleep and therefore not available for producing saliva samples. The lack of data from the entire time of melatonin secretion makes it impossible to use any threshold calculation that depends on the overall amplitude of the pulse.

We developed two additions to the existing physiologically-based model of melatonin: (a) a photic effect to predict melatonin suppression by light and (b) a saliva compartment to predict salivary melatonin concentration. This model was then integrated with an existing mathematical model of the effect of light on the human circadian pacemaker [START_REF] Kronauer | Quantifying human circadian pacemaker response to brief, extended, and repeated light stimuli over the photopic range[END_REF] and used to estimate melatonin synthesis onset and offset, amplitude and infusion and clearance rates from data collected under a variety of conditions, including (a) a constant routine procedure (CR, described below); (b) exposures to continuous bright light, continuous dim light and intermittent bright light; and (c) conditions in which melatonin sampling rate was low or only included the rising portion of the melatonin curve. We compare these predictions to other known methods of melatonin phase estimation and discuss which methods are most accurate under a variety of experimental conditions.

Methods

Model Equations

The physiologically-based model of melatonin developed by Brown and colleagues [START_REF] Brown | A mathematical model of diurnal variations in human plasma melatonin levels[END_REF] was based on the physiological data of melatonin synthesis, secretion and clearance from rats; it is assumed that the same or similar physiology exists for humans. The pathway for pineal melatonin synthesis includes the following steps: tryptophan is hydroxylated and then decarboxylated to form serotonin. Arylalkylamine N-acetyltransferase (AANAT) converts serotonin to N-acetylserotonin [START_REF] Weissbach | Biosynthesis of melatonin: enzymic conversion of serotonin to N-acetylserotonin[END_REF], which is then O-methylated by hydroxindole-Omethyltransferase to produce melatonin [START_REF] Axelrod | Enzymatic O-methylation of N-acetylserotonin to melatonin[END_REF]. Melatonin is not stored in the pineal, which is on the body side of the blood brain barrier, and readily diffuses into the bloodstream [START_REF] Rollag | Route of melatonin secretion in sheep[END_REF], where it is transported through the blood, partially bound to albumin [START_REF] Cardinali | Binding of melatonin to human and rat plasma proteins[END_REF]. The hormone is cleared by the liver and excreted through the urine [START_REF] Kopin | The fate of melatonin in animals[END_REF]. The mathematical model of this synthesis is a firstorder kinetic process with two compartments: the pineal compartment and the plasma compartment. These two compartments are represented by two first-order differential equations:

) ( ) ( ) ( 1 1 t A t H dt t dH IP + - = β (1) ) ( ) ( ) ( 2 1 2 t H t H dt t dH CP IP β β - = ( 2
)
where β I is the rate of infusion of melatonin into the plasma, β C is the rate of clearance of melatonin from the plasma, H 1 (t) is the concentration of melatonin in the pineal compartment and H 2 (t) is the concentration of melatonin in the plasma. A(t) is the AANAT activity at time t and is calculated by the following function:

( ) ( ) ( )            - - = - - - - - - off

Data Analysis

Several datasets were used to develop and test the additions to the physiologically-based model of melatonin. Fluorescent white light was used in all studies (see [START_REF] Gronfier | Entrainment of the human circadian pacemaker to longer-than-24h days[END_REF] for details, including light spectrum) . For all protocols, circadian phase was assessed during constant routines (CRs) using a 2-harmonic plus noise fit of Core Body Temperature (CBT min ) (23) and melatonin as phase markers. During a CR, individuals remain awake in bed under dim light with frequent small meals; the CR is designed to minimize changes known to affect markers of the circadian pacemaker (24) [START_REF] Duffy | Entrainment of the human circadian system by light[END_REF]. CBT min was used as the phase marker during the experiment to time events because CBT data are available in real time and can be analyzed immediately, unlike melatonin which must first be assayed. A description of each dataset follows: In each condition, the stimulus was preceded and followed by 30 minutes of <1 lux. The stimulus day was followed by a 64-hr CR. The light level during wakefulness was ~1.5 lux from the start of the initial CR until the end of the second CR, except during the light exposure session. Although the IBL stimulus represented only 23% of the total duration of the continuous light exposure, the phase shift of DLMO 25% , their marker of melatonin phase, during IBL exposure was 74% as great as during a CBL stimulus. Melatonin suppression was observed in the CBL and IBL conditions ( 26), but was not reported in [START_REF] Gronfier | Efficacy of a single sequence of intermittent bright light pulses for delaying circadian phase in humans[END_REF]. The data from this study were used to introduce the effect of light in the pineal melatonin concentration equation (Equation 1). The plasma melatonin concentrations from the light exposure days were used to develop and test the addition of a light suppression effect to the melatonin model. The plasma melatonin concentrations from the initial CR were used to determine the entrained phase angles between CBT min and SynOn and SynOff. Baseline AUC was calculated from the 4-hr before the melatonin midpoint during the initial CR and the AUC during the light exposure was calculated using the same clock hours 24 hours later. AUC was calculated using the trapezoidal method. These data were used to validate the addition of a light suppression effect to the melatonin model.

Intermittent

3-Pulse

Salivary Melatonin Concentration dataset

Both salivary and plasma samples were collected for melatonin assay from totally blind subjects scheduled to various inpatient protocols. The details of these protocols were described in Klerman et al. 1998 (28). Data from these totally blind subjects were not confounded by ocular light suppression of melatonin concentrations. The salivary melatonin data were used to test and validate the addition of a salivary compartment to the melatonin model.

Incorporation of an effect of light on melatonin

Previous Model

In 1999, Kronauer et al. ( 16) described a mathematical model of the effect of light on the human circadian pacemaker with a dynamic stimulus processor to represent the effect of light as a chemical reaction in which photons entering the retina activate "ready elements" (Process L) that are converted into a light drive that sends a signal to the circadian pacemaker in the SCN. Once the ready elements have been used, they are recycled back to a ready state to await reactivation.

The dynamic stimulus processor assumes a prompt response occurring at the time the light source is turned on, with pseudosaturation occurring within 10 -15 minutes before the drive returns to a lower steady state response until the light stimulus ends.

In this model, this light drive acts on the circadian pacemaker to produce phase shifts. The equations of Process L represent the physiological process by which light initiates a chemical reaction within the photo-pigments of the retinal photoreceptors that transmit a photic signal from the retina through the RHT to the circadian pacemaker within the SCN. This chemical reaction assumes a forward rate constant

p I I       = 0 0 α α (4)
where α 0 = 0.05, I 0 = 9500 and p = 0.5. Equation 4 converts ready elements into a drive B ˆ onto the pacemaker such that

( )α n G B - = 1 ˆ ( 5
)
where G is a scaling constant and n is the fraction of elements in the system that are "used".

Used elements are recycled back into the ready state at a rate of β. In general, the rate at which elements are activated (processed from "ready" to "used") in Process L is given at any time by the formula

( ) [ ] n n n β α - - = 1 60 & . ( 6
)
The drive generated in Equation 5 acts onto the dynamic circadian pacemaker, Process P.

Process P is divided into two components: the circadian pacemaker and the circadian sensitivity modulator. The coupled pacemaker equations are

      +       + + + = B x x x x x c 7 3 105 256 3 4 3 1 12 µ π & (7)                   +       - = kB x qBx x x c c 2 99729 . 0 24 12 τ π & (8)
Equations 7 and 8 represent a higher-order limit cycle oscillator, and the coefficients in Equation 7were chosen so that the amplitude of the limit cycle would equal 1.00, as in previous iterations of the model [START_REF] Kronauer | A quantitative model for the effects of light on the amplitude and phase of the deep circadian pacemaker, based on human data[END_REF].

B ˆ enters Process P via the circadian sensitivity modulator, where it is converted into a direct drive, B, onto the pacemaker that is dependent on the state variables x and x c such that

( )( ) c x x B B 4 . 0 1 4 . 0 1 ˆ-- = ( 9
)
to characterize the feature that the human circadian pacemaker has varying sensitivity to light throughout the circadian day [START_REF] Jewett | Human circadian pacemaker is sensitive to light throughout subjective day without evidence of transients[END_REF].

To compare the results of the model to core body temperature (CBT) data in humans, a phase relationship was derived [START_REF] May | A new mathematical definition of CBTmin improves model predictions of the effect of light on the circadian pacemaker in amplitude suppression protocols[END_REF] such that Time of CBT min = Time of φ xcx + φ ref [START_REF] Brainard | Action spectrum for melatonin regulation in humans: Evidence for a novel circadian photoreceptor[END_REF] where CBT min is as described above, φ ref = 0.97, and φ xcx is defined as the polar phase angle between the state variables x and x c such that

°- = 7 . 170 arctan x x c ( 11
)

Incorporation of light into melatonin model

To incorporate an effect of light into the melatonin model, the AANAT activity, A(t), in Equation 3was redefined as:

) 1 )( ( ) ( B m t A t A - = ′ ( 12
)
where m is a constant whose value was determined by parameter fitting to the Intermittent Bright Light dataset and validated on the 1-Pulse Intensity Response Curve dataset.

Parameter Fitting

Parameter fits to the data were derived by a similar method to the one described in the original publication [START_REF] Brown | A mathematical model of diurnal variations in human plasma melatonin levels[END_REF]. The following equation is used to relate the model predictions to plasma melatonin data:

i i t t H y i ε + = ) ( 2 ( 13
)
where H 2 (t i ) is the plasma melatonin level at time t i and ε i is an approximate Gaussian random variable with mean 0 and variance 2 The -2 log-likelihood of the model is

[ ] [ ] { } ∑ = + - = - N i t i t K t H t H y t H y f i i 2 2 2 1 2 1 ) ( , | ) ( , | log 2 σ θ θ ( 15
)
where

∑ = + - = N i t i N K 2 2 2 log ) 1 ( σ π (16)
The maximum likelihood estimate of θ, , θ is computed by minimizing the likelihood. However, in the original model, a Bayesian analysis was used to calculate 2 i t σ . This analysis was based on the counts per minute recorded from the immunoassay. These data were not available for this analysis, nor are they available for many clinical applications. Therefore, we derived 2 i t σ from the intra-assay coefficient of variance obtained from the laboratory at which the melatonin immunoassays were run.

To integrate this model with an existing mathematical model of the effect of light on the human circadian pacemaker [START_REF] Kronauer | Quantifying human circadian pacemaker response to brief, extended, and repeated light stimuli over the photopic range[END_REF] so that it can predict synthesis onset and offset under a variety of simulated protocols, it is necessary to assign average parameter values to β IP , β CP , α, λ, and A. We used the CR data from the Intermittent Bright Light dataset and the 3-Pulse Phase-Resetting Protocol. The individual fits of these 51 subjects were calculated using the parameter fitting method described above and each parameter was averaged across subjects to determine the average parameter value. The final average parameter estimates (mean ± s.e.m.) are: β IP = 0.047 ± 0.004, β CP = 0.025 ± 0.002, α = 0.067 ± 0.365, λ = 0.588 ± 0.113, and A = 4.834 ± 0.038. The parameter value of m in Equation 12was determined by a best-fit to the Intermittent Bright Light dataset. A value of m = 7 was found to provide the best fit for all subjects.

Saliva Compartment

The design of the saliva compartment was modeled by the same principles as the pineal and plasma compartments were created in the original mathematical model of melatonin. It is assumed that all melatonin is released from the pineal gland into the plasma [START_REF] Iguchi | Melatonin serum levels and metabolic clearance rate in patients with liver cirrhosis[END_REF]. In the plasma, melatonin partially binds to albumin and is eventually cleared through the liver. Melatonin that does not become bound to albumin readily diffuses into the saliva [START_REF] Kennaway | Circadian rhythm of free melatonin in human plasma[END_REF]. Melatonin in the saliva is eventually cleared through unknown processes. A third equation is added to account for the infusion of melatonin from the plasma into the saliva compartment and for the clearance from saliva. This equation is:

) ( ) ( ) ( 3 2 3 t H t H dt t dH CS IS β β - = ( 17
)
where β IS is the rate of infusion from the plasma to the saliva compartment and β CS is the rate of clearance from the saliva. The plasma compartment equation (Equation 1) has also been modified to include the infusion from the plasma to the saliva. This equation is:

) ( ) ( ) ( ) ( 2 2 1 2 t H t H t H dt t dH IS CP IP β β β - - = (18)
A schematic of the model including the saliva compartment is in Fig. 1.

Comparison of Model Simulations to Data

To determine which model fit the experimental data better, both the adjusted R 2 (R 2 adj ) and the Akaike Information Criterion (AIC) were calculated. R equals the covariance between the experimental data and the model divided by the product of the standard deviation of the experimental data and the standard deviation of the model predictions; R 2 therefore has the same units as variance. R 2 adj is calculated as ( )

1 1 2 2 - - - - k n k R R ( 19
)
where n is the number of samples and k is the number of parameters estimated for each model.

The inclusion of "k" in the calculation imposes a penalty for adding new parameters to the model. A higher R 2 adj indicates a better fit to the experimental data. When the degrees of freedom n -k -1 was less than or equal to zero, R 2 was reported in place of R 2 adj .The AIC was calculated as n log ( 2 σ ) + 2k, where 2 σ is the residual sum of squares, and n and k are defined as in the R 2 adj calculation [START_REF] Anderson | Null hypothesis testing: problems, prevalence, and an alternative[END_REF]. Similarly to the R 2 adj calculation, the "2k" term in the AIC formula imposes a penalty for adding new parameters to a model. A lower AIC indicates a better fit to the experimental data. For the original model (Equations 1 -3), k = 7. With the addition of the light suppression term (Equation 12), the number of parameters fit increases to 8. For the salivary melatonin model (Equation 18), two additional parameters were added to the model; therefore, k = 10.

Results

For the 51 subjects from the Intermittent Bright Light dataset and the 3-Pulse Phase-Resetting Protocol, SynOn occurred 8.26 ± 1.58 hr before CBT min while SynOff occurred 1.28 ± 3.33 hr before CBT min on average. Thus the time of SynOn is defined as:

Time of SynOn = Time of φξ cxon + φ ref (19)
where SynOn is the synthesis onset, φ ref = 0.97, and φ xcxon is defined as the polar phase angle between the state variables x and x c of the circadian model such that

°- = 3 . 65 arctan x x c (20)
Similarly, the time of SynOff is defined as:

Time of SynOff = Time of φ xcxoff + φ ref ( 21
)
where SynOff is the synthesis offset, φ ref = 0.97, and φ xcxoff is defined as the polar phase angle between the state variables x and x c such that

°- = 2 . 166 arctan x x c ( 22
)
The model was fit to data from the Intermittent Bright Light dataset with CBL, IBL and VDL conditions. While both the original melatonin model and that with the incorporated light effect can predict the data from the continuous VDL condition (Fig. 2, bottom panel), only the latter model can simulate the melatonin suppression of the CBL or IBL conditions (Fig. 2, top and middle panels). Table 1 compares the R 2 and AIC of both models with respect to the experimental data for each subject in each condition. SynOn indicates the time when melatonin synthesis onset occurred and SynOff indicates the time when melatonin synthesis offset occurred.

A indicates the maximum N-acetyltransferase activity. R 2 and AIC are two goodness of fit measures; a higher R 2 and a lower AIC indicate a better fit to the data. For subjects in the CBL and IBL conditions, the R 2 values increase and the AIC values decrease from the original melatonin model to the model with the incorporated light effect, indicating that the model with the incorporated light effect is a better fit to the data. For the VDL condition, the R 2 values are approximately equal for both models, while the AIC increases for the model with the incorporated light effect. Because the light level in this condition is low, the light suppression is not activated; the AIC values increase due to the "2k" parameter penalty in the AIC equation.

This model was validated by fitting the model to data from the 1-Pulse Intensity Response Curve

Protocol with m = 7 as found by fits to the Intermittent Bright Light dataset as described above. The data from one subject are plotted in Fig. 5 against model predictions. Table 3 reports the R 2 , AIC and SynOn for each subject in which the rising portion of the salivary melatonin concentration was available. In some cases, the SynOn value predicted was not a reasonable estimate of the true clock time at which SynOn would be estimated to occur based on other methods. The possible reasons for this include a poor fitting procedure, a bad choice of initial parameters, data that are too noisy, or not enough data points.

Another way to assess the model's ability to calculate SynOn when only the rising portion of a curve is available is to truncate full plasma and salivary melatonin curves and fit the model. A total of 16 full salivary melatonin curves were available for this analysis from the Salivary Melatonin Concentration dataset. The curves were truncated at the peak value of melatonin obtained. Again, an arbitrary SynOff, equal to the time at which the final available data point was measured, was used. Table 4 reports the R 2 and the SynOn values estimated from the truncated curves compared to the SynOn values estimated from the full curves. There is larger variability between the two SynOn values, which suggests that a full curve may be necessary to obtain an accurate SynOn.

The existing mathematical model of the human circadian pacemaker uses CBT min as a phase marker of the circadian system, although the CBT rhythm itself is not an output of the model. By incorporating melatonin as both a phase marker and an output rhythm of the circadian pacemaker, the model can predict both phase and melatonin suppression. Fig. 6 is an example of the output of the melatonin model generated by simulating the 3-Pulse Phase-Resetting Protocol, as described above. The phase response to the light stimulus is determined by subtracting initial phase found during the initial pre-stimulus CR from the final phase found during the poststimulus CR. Additionally, melatonin suppression is observed on the three light exposure days, which affects the estimate of melatonin phase. By incorporating melatonin into our existing model of the circadian pacemaker, the model can predict both melatonin phase and melatonin suppression under a number of experimental protocols.

Discussion

The original physiologically-based mathematical model of melatonin was able to estimate synthesis onset, offset and amplitude from plasma melatonin samples collected under dim light conditions. We have updated this model to include an effect of ocular light exposure. This model can be used to estimate synthesis onset, offset and amplitude from plasma melatonin samples collected under a variety of lighting conditions, including exposure to intermittent bright light pulses. Phase estimates are often not calculated on light exposure days because threshold methods cannot account for the change in shape and amplitude of the melatonin profile that is observed during the light exposure. However, the model with an incorporated light effect can account for the changes in melatonin concentration observed in response to a light pulse. Furthermore, the additional parameters of the model, including infusion rate of melatonin into the plasma and clearance rate from the plasma, may provide useful information about the dynamics of melatonin suppression by light.

We also added a salivary melatonin compartment to the existing model. The structure of the salivary melatonin compartment is similar to the plasma melatonin compartment and includes rates of infusion from the plasma to saliva and clearance from the saliva. Thus, this model can be used to estimate synthesis onset, offset and amplitude from salivary melatonin samples when plasma melatonin samples have not been collected. Because full salivary melatonin profiles are not always available, we additionally tested whether the model was able to predict accurately the synthesis onset from the rising portion of salivary melatonin only. As can be seen in Table 4, the results of this analysis were inconclusive: for some melatonin pulses analyzed, the SynOn estimates obtained from the rising portion of the curve were within 2 minutes of the SynOn estimates obtained from the full curve. For other melatonin pulses, this was not the case. Further analysis with a larger set of data is necessary. The updated melatonin model was incorporated into an existing model of the circadian pacemaker. The original model of the circadian pacemaker reported CBT min as the phase marker; however, melatonin is the preferred marker rhythm and most experimental results are reported using melatonin. Therefore, the addition of this melatonin model to the circadian pacemaker model facilitates comparisons of experimental results to model predictions. Additionally, melatonin suppression by light is often used to evaluate the sensitivity of the pacemaker to varying durations, timings and intensities of light. This model of the circadian pacemaker can predict melatonin suppression as well as phase shifts of the melatonin rhythm.

In the original publication of the melatonin model [START_REF] Brown | A mathematical model of diurnal variations in human plasma melatonin levels[END_REF], it was reported that estimates of synthesis onset and offset were comparable to melatonin phase markers estimated by threshold methods. Because existing threshold methods depend on the shape and amplitude of the melatonin pulse, they are not adequate measures to estimate phase under experimental conditions in which the shape and amplitude of the melatonin profiles changes. We propose that a physiologically-based method to estimate melatonin phase is more appropriate because these changes in shape and amplitude are considered in the parameters of the model. Furthermore, this method of analysis can provide additional information about the melatonin profile that cannot be measured with threshold-based methods.

Many studies have shown that melatonin suppression by light (10) (11) (12) (13), as well as circadian phase-shifting measured from the melatonin rhythm (12) (35) (36) [START_REF] Wright | Differential effects of light wavelength in phase advancing the melatonin rhythm[END_REF], is differentially sensitive to light wavelength. In particular, action spectra that have been constructed from melatonin suppression data reveal peak sensitivity around 460 nm [START_REF] Brainard | Action spectrum for melatonin regulation in humans: Evidence for a novel circadian photoreceptor[END_REF] [START_REF] Thapan | An action spectrum for melatonin suppression: Evidence for a novel non-rod, non-cone photoreceptor system in humans[END_REF] Future work on the circadian model will include the incorporation of wavelength of light information to predict the effect of narrowband and broadband monochromatic light stimuli on melatonin suppression and melatonin phase shifts.

A user-friendly version of the model detailed in this paper is available for download at: http://dsm.bwh.harvard.edu/bmu/melatonin/ or can be obtained directly from the author.

Figures and

  Bright Light (IBL) dataset Gronfier et al. (3) performed a series of experiments to determine the response of the circadian pacemaker to a sequence of intermittent bright light pulses compared to the response of the pacemaker to a continuous bright light or very dim light pulse. In these experiments, 20 subjects were scheduled to two 24.0-hr days in ~90 lux followed by a 26.2-hr CR. After the initial CR, subjects were scheduled to 8 hr of sleep in 0 lux. For the next wake episode, subjects were randomly assigned to a 6.5-hr stimulus, centered ~3.5hr before CBT min , in one of three lighting conditions: (1) continuous bright light (CBL) at ~9500 lux (n = 6); (2) intermittent bright light (IBL) consisting of six 15-minute bright light pulses of ~9500 lux separated by 60 minutes of very dim light of < 1 lux (n = 7); (3) continuous very dim light (VDL) of <1 lux (n = 7).

  Phase-Resetting Protocol Duffy et al. (27) performed a series of experiments in which 31 subjects were scheduled to two 24.0-hr days in ~90 lux followed by a 26 to 33-hr CR, depending on the desired center of the bright light stimulus. After the initial CR, subjects were scheduled to 8 hr of sleep in 0 lux. During the following three wake episodes, subjects were exposed to a 5-hr ~9500 lux pulse centered in the middle of the 16-hr wake episode, which occurred at different circadian phases to either phase-advance or phase-delay the pacemaker. The 3-cycle bright light stimulus was followed by a 40-hr CR. CBT and plasma melatonin rhythms were assessed throughout the protocol. The light level during wakefulness was ~10 lux from the start of the initial CR until the end of the second CR, except during the bright light stimuli. The plasma melatonin concentrations from the initial CR were used to determine the entrained phase angles between CBT min and SynOn and SynOff. 1-Pulse Intensity Response Curve Protocol Zeitzer et al. (9) performed a series of experiments to determine the response of the pacemaker to light pulses of varying intensities. In these experiments, 21 subjects were scheduled to two 24.0-hr days in 150 lux followed by a ~50-hr CR. After the initial CR, subjects were scheduled to 8 hr of sleep in 0 lux. During the following wake episode, subjects were randomly assigned to receive a single 6.5-hr light pulse centered ~3.5 hr before CBT min with intensities ranging from 3 lux to 9100 lux. The stimulus day was followed by a 30-hr CR. The light level during wakefulness was ~10 lux from the start of the initial CR until the end of the second CR, except during the stimulus. Melatonin suppression was calculated as the difference between the area under the curve (AUC) of the baseline melatonin pulse and the melatonin pulse during the light exposure divided by the baseline AUC.

A

  phase relationship was derived to compare the results of the model to melatonin and CBT data in humans. CBT min and melatonin SynOn and SynOff were calculated from the initial CR data from the 3-Pulse Phase-Resetting Protocol and the Intermittent Bright Light dataset.

σ

  determined by the melatonin immunoassay.The objective is to estimate θ,

Fig. 3

 3 Fig. 3 plots the melatonin suppression calculated from the model fits as a function of light

Fig. 1 :

 1 Fig. 1: Schematic diagram of the melatonin model adapted from Brown et al. 1997.

Fig. 2 :Fig. 3 :

 23 Fig. 2: The melatonin concentration of a subject exposed to a 6.5-hr continuous bright light

Fig. 4 :

 4 Fig. 4: The 24-hr salivary melatonin profile of a blind subject fit by the salivary melatonin model. For these data, the fit had R 2 = 0.975 and AIC = 37.009. Arrows pointing down indicate

Fig. 5 :

 5 Fig. 5: The salivary melatonin concentration in a blind subject in which data from the rising

Fig. 6 :

 6 Fig. 6: Predictions of the plasma melatonin concentration (solid line, based on average parameter

Fig. 1 :Fig

 1 Fig. 1:

Fig. 3 :Fig. 4 :Fig. 5 :Fig. 6 :

 3456 Fig. 3:

  where λ is the AANAT activity growth, α is the AANAT activity decline, t on is the onset time of increased AANAT activity and t off is the offset time of growth in AANAT activity and A is the maximum AANAT activity (pmol•l -1 •min 1 ). Although the original model was based on the assumption that AANAT activity is the rate limiting step of melatonin synthesis, more recent studies indicate that this may not be the case [Liu and Borjigin 2005 J. Pineal Research -need to add to refman].

	t ( A	)	A	1 1		e e	t λ λ	off t	t	on on t	t	on	≤	t	<	t	off	(3)
			Ae	α	t	t					t	off	≤	t

  Table 2 reports the R 2 , SynOn and SynOff for each subject in whom a full salivary melatonin profile was available. Table 2 also reports the SynOn and SynOffvalues for each subject calculated from plasma melatonin, when available. In general, the AIC values are lower for the salivary than for plasma melatonin data. This does not indicate that the fits by the salivary melatonin model are better; the difference in values may be due to difference in sampling rates of saliva and plasma melatonin concentrations. The salivary melatonin data were also more variable than the plasma data, which may also explain why the R 2 values are lower in the salivary melatonin fits.When only a portion of a melatonin pulse is available for analysis, relative threshold methods such as DLMO 25% cannot be used to estimate circadian phase. Often under clinical settings, only the rising portion of a salivary melatonin pulse is obtained, because collecting samples throughout the night would require interrupting the sleep episode. Only the rising portion of the melatonin profile was analyzed for 14 melatonin profiles from 4 subjects in the Salivary Melatonin Concentration dataset. These data were fit by the model to determine whether an accurate estimate of SynOn could be extracted. Since SynOff cannot be estimated from the data, an arbitrary SynOff, equal to the time at which the final available data point was measured, was used in the model, since solving Equation 3 requires t off . Therefore, only

	θ =	(	β	IP	,	β	CP	,	β	IS	,	β	CS	,	α	,	λ	,	, A	Syn	on	)	was estimated.

Table 1 :

 1 Estimates of the original melatonin model and the model with incorporated light effect to plasma melatonin data from subjects scheduled to a 6.5-hr continuous bright light (CBL), intermittent bright light (IBL) or very dim light (VDL) exposure. "--" indicates that data were missing or that an accurate value could not be obtained.

			Original melatonin model		Model with incorporated light effect
	Subject	SynOn (clock hour)	SynOff (clock hour)	A(t) (pmol•L 1 • min -1 ) -	R	2	adj	AIC	SynOn (clock hour)	SynOff (clock hour)	A(t) (pmol•L 1 • min -1 ) -	R	2	adj	AIC
	Condition
	BL
	19a1v	20.87	8.11	1.19	0.13	3892.12	20.2	6.10	1.29	0.96	1844.38
	2006v	20.09	21.91	7.70	0.38	5420.69	20.01	4.17	11.28	0.98	3238.42
	2012v	5.44	7.87	3.71	0.99	3536.74	22.96	8.03	3.71	0.99	3541.13
	2036v	4.39	8.30	4.37	0.99	3279.80	4.43	5.78	2.06	0.99	2912.51
	2041v	--	20.51	1003.11	-0.30	6266.01	18.55	5.69	5.40	0.99	3172.92
	2043v	0.76	3.16	4.05	0.82	2789.14	--	6.55	0.72	0.96	1958.00
	Condition
	IBL
	19b9v	21.42	7.45	0.89	0.94	2376.36	22.23	7.59	1.07	0.97	1852.74
	19c4v	20.8	8.09	2.83	0.71	5178.28	20.95	7.29	8.53	0.90	4480.78
	2004v	23.6	8.01	1.01	0.82	5752.39	23.16	8.69	11.52	0.97	4589.86
	2024v	21.3	8.37	4.00	0.94	4631.57	21.79	7.37	3.98	0.98	3893.72
	2035v	20.54	5.78	3.84	0.91	4035.70	20.42	5.42	3.06	0.98	3019.08
	19c9v0t2	20.62	6.72	5.66	0.81	4562.74	20.58	6.17	4.01	0.90	4151.10
	2042v	21.76	--	1.14	0.93	4213.36	21.7	--	2.15	0.96	3783.21
	Condition
	VDL
	19b1v	22.17	6.38	6.04	0.90	4507.15	22.08	7.27	3.42	0.88	4636.21
	19c3v	22.13	5.88	4.05	0.81	5032.53	21.84	7.25	4.26	0.82	4995.12
	2008v	21.41	2.71	2.86	0.92	3942.86	21.31	2.17	2.68	0.92	3913.54
	2018v	22.50	7.26	4.44	0.94	4410.99	22.3	7.83	3.85	0.93	4547.55
	2038v	19.68	20.66	11.77	0.96	3601.56	20.33	6.45	10.81	0.87	4474.42
	2039v	20.09	23.39	1.59	0.94	3665.99	20.21	5.16	4.42	0.90	3988.82
	2056v	23.10	6.55	3.49	0.96	3710.22	22.97	7.14	2.94	0.96	3697.57

Table 2 :

 2 Estimates of the model with an incorporated salivary melatonin component to full salivary melatonin profiles compared to estimates of the original plasma melatonin model to full plasma melatonin profiles in the same subjects. "--" indicates that data were missing or that an accurate value could not be obtained.

	Salivary

Table 3 :

 3 Estimates of the model with an incorporated salivary melatonin component to the rising portion of the salivary melatonin concentration. "--" indicates that data were missing or that an accurate value could not be obtained.

			SynOn
		Day of	(clock
	Subject	Protocol	hour)	R 2	AIC
	1415U7T6	7	4.30	1.00	11.56
		8	5.60	1.00	17.42
		9	3.02	0.80	40.10
		15	5.98	0.92	39.18
		16	7.05	0.97	32.48
		21	8.03	0.98	15.90
		30	--	0.97	23.60
		31	--	1.00	14.00
		33	--	1.00	6.76
		38	3.15	1.00	9.82
	1459U5T6	5	--	0.79	24.23
	1459U5T7	10	20.87	0.98	22.91
		11	--	1.00	10.12
	1531U6T6	8	18.97	0.99	24.49

Table 4 :

 4 Estimates of the model with an incorporated salivary melatonin component to the rising portion of full salivary melatonin profiles. "--" indicates that data were missing or that an accurate value could not be obtained.

			SynOn (clock			SynOn (clock
			hour)			hour)
			estimated			estimated
		Day of	from rise			from full
	Subject	Protocol	only	R 2	AIC	curve
	1415U7T6	4	5.15	1.00	14.04	5.14
		10	6.27	0.97	24.08	4.98
		11	--	0.99	22.49	6.14
		17	6.68	0.91	32.12	6.96
		18	--	0.91	26.78	7.91
		20	7.00	1.00	16.46	6.98
		46	8.50	1.00	20.06	--
	1451U4T4	3	22.98	0.96	18.73	0.06
		13	0.19	0.96	20.05	1.96
	1459U4T4	6	17.08	0.98	22.04	20.09
	1459U5T6	2	20.03	0.99	19.08	20.02
		6	21.00	0.95	23.01	22.49
		7	20.88	0.80	32.28	--
	1459U5T7	5	--	0.99	21.96	--
	1531U6T6	2	11.42	1.00	14.24	--
		7	17.02	0.98	33.47	--

Acknowledgements

We thank Drs. Jeanne F. Duffy and Charles A. Czeisler for the use of their data. This work was supported by National Aeronautics and Space Administration (NASA) Cooperative Agreement NCC9-58 with the National Space Biomedical Research Institute HPF 00203 and HPF 00405 and by NASA Grant NAG69-1035. The experimental studies were performed in a General Clinical Research Center supported by National Institutes of Health (NIH) Grant MO1-RR02635, NIA-1-R01-AG06072 from the National Institute on Aging and NIMH-1-R01-MH45130 from the National Institute of Mental Health. EBK is supported by NIH K02-HD045459. CG is supported by FP6-EUCLOCK


Keywords: Adaptation, Physiological, physiology, radiation effects, Adult, Body Temperature Regulation, physiology, radiation effects, Circadian Rhythm, physiology, radiation effects, Female, Humans, Light, Male, Melatonin, blood, Photic Stimulation, methods, Photoperiod, Sleep, physiology, radiation effects intermittent light, melatonin, circadian pacemaker, photoreception, phase shift

It has been shown in animal studies that exposure to brief pulses of bright light can phase shift the circadian pacemaker, and that the resetting action of light is most efficient during the first minutes of light exposure. In humans, multiple consecutive days of exposure to brief bright light pulses have been shown to phase shift the circadian pacemaker. The aim of the present study was to determine if a single sequence of brief bright light pulses administered during the early biological night would phase delay the human circadian pacemaker. Twenty-one healthy young subjects underwent a 6.5 h light exposure session, in one of three randomly assigned conditions: 1) continuous bright light of ~9,500 lux; 2) intermittent bright light (six 15-min bright light pulses of ~9,500 lux separated by 60 min of very dim light of < 1 lux); 3) continuous very dim light of < 1 lux. Twenty subjects were included in the analysis. Core body temperature (CBT) and melatonin were used as phase markers of the circadian pacemaker. Phase delays of CBT and melatonin rhythms in response to intermittent bright light pulses were comparable to those measured after continuous bright light exposure, even though the total exposure to the intermittent bright light represented only 23 of the 6.5 h continuous exposure. These results % demonstrate that a single sequence of intermittent bright light pulses can phase delay the human circadian pacemaker and shows that intermittent pulses have a greater resetting efficacy on a per minute basis than does continuous exposure.

Efficacy of a single sequence of intermittent bright light pulses for delaying circadian phase in humans

Gronfier Claude1 * , Wright Kenneth P. 1 2 , Kronauer Richard E. 3 , Jewett Megan E. 1 , Czeisler Charles A. 1

Division of Sleep Medicine 1 Brigham and Women s Hospital -Harvard Medical School

INTRODUCTION

The daily alternation of light and darkness is the primary synchronizer of endogenous circadian oscillators to the 24-hour day. The role of light as a phase resetting agent of the human circadian pacemaker has been demonstrated in a variety of experiments using single or multiple sequences of continuous bright light exposure (for review see ref.

( )). The phase-shifting effect of bright light on the human 17 pacemaker has been shown to be phase-dependent ( ). Exposure to several hours of bright light for up to 3 consecutive days have been 33 found to reset the endogenous circadian rhythms of core body temperature, cortisol, melatonin, alertness, performance and TSH by up to 12 hours ( , , , , , , , , , , , ). Moreover, it has been shown that even room light can significantly reset the 1 9 12 14-16 [START_REF] Dawson | Improving adaptation to simulated night shift: Timed exposure to bright light versus daytime melatonin administration[END_REF] [START_REF] Dawson | Phase resetting of the human circadian pacemaker with use of a single pulse of bright light[END_REF] 21 [START_REF] Duffy | Phase-shifting human circadian rhythms: Influence of sleep timing, social contact and light exposure[END_REF] [START_REF] Eastman | Circadian rhythm adaptation to simulated night shift work: Effect of nocturnal bright-light duration[END_REF] 30 48-50 55 phase of the human circadian clock ( , ). In fact, Zeitzer ( ) found that half of the maximum phase delay resetting response 5 6 et al. [START_REF] Zeitzer | Sensitivity of the human circadian pacemaker to nocturnal light: Melatonin phase resetting and suppression[END_REF] achieved by a single exposure to bright light (~9,000 lux) during the early biological night can be obtained by ~1 of this light intensity, % i.e., ordinary room light of ~100 lux.

It has been reported in studies of nocturnal animals that the mammalian circadian system is responsive to brief exposures of light.

Nelson and Takahashi ( ) demonstrated that pulses of bright light as short as 3 seconds produced measurable phase shifts in the golden 43 hamster. They also showed that extending the duration of the stimulus beyond 5-min produced little additional phase shift. Furthermore, they found that the shift obtained by two 5-min light pulses separated in time by up to one hour produced a phase resetting response similar to that of a single 5-min pulse ( ). Others have reported that the circadian system can be reset by a chain of very brief light pulses 42 of 2 milliseconds in mice ( ), and 10 microseconds in duration in rat ( ). The structure of the physiological system responsible for such a 51 2 dynamic resetting of the circadian pacemaker is not yet known. They have been hypothesized by Kronauer and colleagues ( , ) to 38 39 involve light sensitive elements that are set to a used state after being exposed to light, and slowly recover their initial ready state after a " " " " period of recycling after light removal.

To date there has been few systematic evaluation of the relative effectiveness of different durations of light exposure in the human circadian system ( ), although the results of studies on the resetting response to intermittent bright light exposure suggest that humans are 25 also responsive to brief pulses of light. Indeed, three consecutive days of exposure to six 40-min pulses of bright light (~5,000 lux) have been shown to phase delay circadian rhythms in humans ( ). In other studies, three to five consecutive days of exposure to intermittent 3 pulses of bright light (3000 to 11,000 lux) as brief as 5 min in duration have been shown to induce phase advances of the endogenous circadian pacemaker in humans ( , , ). These findings suggest that humans, like animals, are sensitive to brief bright light pulses.

However, interpretation of the temporal dynamics of the resetting response was complicated by the use of multiple, consecutive days of light exposure in those studies.

In the present study, we therefore quantified the phase-resetting efficacy of a single sequence of intermittent bright light pulses given during the biological night, as compared to subjects exposed to continuous bright light or to very dim light at the same circadian phase.

Based on prior results, we anticipated a robust phase delay in the groups of subjects exposed to both continuous ( ) and intermittent ( ) 33 45 bright light vs. small drifts in phase in the group of control subjects exposed to continuous very dim light due to the near-24-h period of the human circadian pacemaker ( , , , ). In addition, based on a mathematical model of the effect of light on the human circadian 13 33 34 54 pacemaker developed by Kronauer and co-workers ( , , ), we predicted that the resetting effect of the intermittent light exposure 32 38 39 would be more efficient than continuous bright light exposure, as measured by the amount of shift induced per minute of bright light exposure.

MATERIAL AND METHODS

Subjects

Twenty-one healthy subjects participated in the study; 20 subjects 24.3 3.9 (SD) years old, 15 males, 6 females, BMI 22.8 2.3(SD)

[ ± ±
were included in the analysis due to non-compliance to constant routine procedures by one subject. The subjects had no medical, ] psychiatric, or sleep disorders as determined by history, physical examination, electrocardiogram, blood and urine biochemical screening tests, and psychological screening questionnaires (Minnesota Mutliphasic Personality Inventory and Beck Depression Inventory). A staff psychologist interviewed subjects and those with a history of or a current psychiatric pathology were excluded. Subjects reported that they were not taking any medication and were instructed to abstain from the use of alcohol, nicotine, recreational drugs, and foods or beverages containing caffeine for three weeks prior to the study. All participants were drug free at the time of study as verified by a comprehensive toxic analysis conducted upon admission to the laboratory. All experimental procedures were carried out in accordance with the principles of the Declaration of Helsinki, and the protocol was approved by the Human Research Committee at the Brigham and Women s Hospital.

'

Prior to beginning the protocol, all participants gave informed, written consent.

Overall study design

Subjects were required to maintain a regular 8:16 h sleep: wake schedule at home for at least 3 weeks prior to admission to the laboratory. In order to ensure compliance to this protocol, subjects were required to call into a date/time-stamped answering machine just prior to going to bed and immediately upon awakening, and the times were compared to sleep-wake logs on the day of admission. In addition, wrist activity and light exposure were monitored for 1 week immediately prior to admission to the laboratory (Actiwatch-L, Mini

Mitter, Sunriver, OR, USA) and were used to verify the stability of their sleep-wake cycle during that last week and throughout the entire impatient protocol.

Upon admission to the study on experimental day 1, subjects were maintained in an environment free from external time cues, including clocks, radios, television, visitors, and sunlight. Subjects maintained contact with staff members specifically trained to avoid communicating time of day or the nature of the experimental conditions to the subjects. Subjects were adapted to the laboratory with three baseline days (Days 1 3), during which time they continued to sleep and wake at their habitual times ( ). In order to assess their -Fig. 1 endogenous circadian phase before the light stimulus and to appropriately center the light stimulus during the following scheduled episodes of wakefulness, subjects underwent a 26.2-h initial constant routine (CR1; procedure described below) from day 4 to day 5. The duration of the constant routine and the associated shift of the sleep episodes 4 and 5 (see ) was designed so that the 6.5-h light

- Figure 1 -
exposure session (centered in the middle of the wake period) would coincide with the delay portion of the human phase response curve ( 33), and was started 9.05 h before habitual bedtime and ended 2.55 h before habitual wake time (corresponding to a center at ~3.5 h before the subjects core body temperature minimum CBTmin ).

' [ ]

On day 5, subjects were randomly assigned to one of the three light exposure conditions (described below). Then, they were scheduled to a final constant routine of 64-h duration (Days 6 9) to re-assess their endogenous circadian phase and thereby estimate the phase shift from the initial constant routine. Subjects were then discharged on day 10, after a ~ 22-h recovery sleep episode.

During wake episodes, subjects were free to move about the suite as desired, except that they were asked not to lie down, nap, or exercise. Subjects compliance with the protocol was monitored by means of closed-circuit cameras and frequent interaction with depending on the experimental condition (described below), <1 lux (~0.5 lux measured horizontally at a height of ~1.37 m) or ~9,500 lux (in the direction of gaze at a wall-mounted target) during the light exposure session. Sleep episodes were conducted in darkness. The experimental suites were equipped with hand-held terminals for on-line event recording, a porthole for 24-h blood sample collection with minimal sleep disturbance, a video camera and a voice-activated audio system.

Constant Routine Procedure

The constant routine (CR) procedure was used to assess the phase of two markers of the endogenous circadian pacemaker: the core body temperature rhythm and the melatonin rhythm. The methodology we used involved a refinement of a technique first proposed by Mills ( ), according to which the subjects are studied under constant environmental and behavioral conditions in order to eliminate et al. [START_REF] Mills | Adaptation to abrupt time shifts of the oscillator s controlling human circadian rhythms[END_REF] or distribute across the circadian cycle the physiologic responses evoked by environmental or behavioral stimuli, such as sleeping, eating, and changes in posture ( ). The CR consisted of a regimen of enforced wakefulness in a semi-recumbent posture in constant dim 22 illumination of ~1.5 lux. Subjects were required to maintain a very low level of physical activity and were not permitted to change their posture throughout the constant routine. This posture was also maintained for urine samples and bowel movements. Nutritional intake was divided into hourly isocaloric snacks to maintain an equal caloric intake across the circadian cycle. Caloric requirements were calculated with use of the Wilmore nomogram ( ) to determine the basal metabolic rate, and adjusted upward by a 10 activity factor. Fluid intake 53

% was calculated for each subject to account for the sedentary nature of the CR. A staff member monitored the subject to help maintain wakefulness and to ensure compliance with the posture and activity level requirements.

Light exposure conditions

Subjects were randomly assigned to one of the three light exposure conditions ( ). The light exposure session of 6.5 hours Fig. 2 duration was centered 5.8 h before habitual wake time (i.e., it started ~1.1 h before habitual bedtime and ended ~2.6 h before habitual waketime) and therefore ~3.5 hours before CBT ( ). The three conditions were: continuous bright light (BL) of ~9,500 lux, defined as min 23 100 bright light; intermittent bright light (IBL) consisting of six 15-minute bright light pulses of ~9,500 lux separated by 60 minutes of % very dim light (<1 lux), defined as 23 bright light; and continuous very dim light (VDL) of <1 lux, defined as 0 bright light. In all three % % conditions, the 30 minutes prior to and the 30 minutes following the light exposure session were conducted in <1 lux. Seven subjects were randomized to each condition. Subjects were seated in a chair from 2 hours before until 2 hours after the end of the light exposure session.

A technician was present at all times during the light exposure session to ensure compliance with the protocol. During the light exposure sessions, subjects were instructed to fix their gaze on a target mounted on the wall of their suite for 5 of every 10 minutes and then allowed free gaze for 5 of every 10 minutes to ensure consistency of light exposure between subjects. This alternating fixed and free gaze episodes continued during the entire duration of the light exposure sessions. Light intensity measurements were taken in the direction of gaze every five minutes to ensure that subjects were exposed to the appropriate light intensity for the appropriate duration. Clear polycarbonate lenses filtered 99.9 of the light in the UV range from the light source. In addition, all subjects wore clear goggles (model S379; luminous % # transmittance 90 , UV absorption >99 , Uvex Safety, Smithfield, RI, USA) during the exposure sessions to bright light for additional % % protection against UV exposure.

Data Collection

Temperature data

A real time, on-line data acquisition system utilizing IBM-PC-compatible computers was employed to monitor and collect temperature data. Core body temperature was continuously recorded via a disposable rectal thermistor (Yellow Springs Instrument, Yellow Spring, OH, USA), and room temperature was recorded by means of a permanent air temperature thermistor; both were recorded every minute.

Hormonal data

Blood samples were collected every 40 minutes during constant routines, through an indwelling intravenous catheter that was inserted into a forearm vein on day 2 of the study, and every 5 to 10 minutes during the light exposure session (data will be presented elsewhere). A solution of heparinized saline (0.45 sodium chloride, 10 units of heparin/ml) was infused at a slow rate (20 ml/h, 200 IU heparin/h) % between blood samples. We provided participants with Ferrous Gluconate (324 mg) pills to be taken at breakfast and dinner: 1) for a minimum of one week prior to participation in the inpatient research protocol; 2) during the inpatient portion of the protocol; and 3) for three weeks after completion of the protocol. Subjects hemoglobin levels were tested every 1 2 day(s) to ensure appropriate levels (>11.0 Osceol WI, USA). The assay sensitivity was 2.5 pg/ml. The average intra-assay coefficients of variation (CVs) were 6.4 below 50 pmol/L, and 4.9 above.

% %

Estimation of circadian phase

For temperature data collected during constant routines, a non-linear least-squares dual-harmonic regression model with correlated noise ( ) was used. The time of the fitted core body temperature minimum (CBT ) was defined as the marker of the temperature phase.

sleep episode, or from changes in posture from flat to the semi-recumbent CR posture. In order to minimize uncertainty in the phase estimate, subjects for whom amplitude of the temperature rhythm was below 0.14 C were not included in the analysis ( ) (see ).

°24

Table 1 In order to minimize the effect of inter-individual differences in circadian period on the phases shift assessments only the first circadian oscillation of melatonin and CBT of the final CR was included in the analysis, and phase was estimated on the same time window length as for CR1.

For melatonin data, a three-harmonics least-square regression analysis was applied to data collected during CR1 ( , thick line) to Fig. 3 estimate the amplitude of each subject s melatonin rhythm ( , dashed line). A threshold of 25 of the peak-to-trough amplitude of the ' Fig. 3 % fitted curve ( ) was calculated ( , dotted line) and applied to data of CR1 and CR2 (first circadian oscillation considered). Three 54 Fig. 3 markers of melatonin phase were computed for each CR data: 1) the midpoint (MP ) between 2) the upward (DLMOn ), and 3) the 25% 25% downward (DLMOff ) crossings of this 25 threshold.

25%

%

Phase shifts of temperature and melatonin rhythms were calculated as the difference between the phase in CR1 and the phase in CR2. exposure, the continuous and intermittent bright light groups demonstrated phase delays of their endogenous circadian rhythms of core body temperature and melatonin ( ). Although the phase shifts observed in the IBL group were slightly smaller than those observed in Fig. 4 the BL group, no statistical differences were found for any marker (

Statistical analysis

Repeated

). Fig. 4 The resetting responses measured with CBT were more variable than those measured with melatonin ( ), however, we found min Table 1 highly significant and positive linear relationships between the markers of the two rhythms (r 0.896, r

CBTmin/DLMOn25% = CBTmin/MP25% = 0.871, r 0.779; p 0.0001 for all three correlation coefficients).

CBTmin/DLMOff25% = ≤

We found no significant linear correlation between observed melatonin phase shifts and the initial phases of light exposure ( )

φ init
within the limited range of initial phases to which subjects were exposed within each condition (r 0.182, r φ init/DLMOn25% = φ init/MP25% = 0.071, r 0.031; p non significant for all three correlation coefficients)

φ init/DLMff25% =
In appraising the variability of the phase shift data presented here ( ) we recognize three independent variability sources whose Table 1 variances should sum: 1) noise in the individual phase assessments whose difference comprises the phases shift, 2) inter-subject differences in the endogenous circadian period that leads to different phase drifts over the two cycles of shift assessments, 3) inter-subject differences in photic response arising either from sensitivity differences or from differences in the timing of the stimulus. Klerman et al. (

) found that the method of single-phase assessment yielding the lowest variance was a normalized level-crossing on the rising front of 36 the melatonin pulse. This variance could be as low as 0.04 h (50 percentile) so that for phase shifts the variance might be 0.08 h . For deviation of circadian period for young adults to be 0.13 h so that a two-cycle phase drift would show a variance of 0.07 h . We recognize 2 that variance estimates are extremely unreliable for data groups of 6 or 7 members. However, the value of 0.15 h for the VDL condition is 2 consistent with the sum of 0.08 h and 0.07 h for sources (1) and ( 2) described above. The variance of 0.20 h for the BL condition suggests that perhaps as little as 0.05 h might be ascribed to source [START_REF] Baehr | Eastman CI Intermittent bright light and exercise to entrain human circadian rhythms to night work[END_REF]. However, by this reasoning the variance of 0.12 h for the IBL

2 2
condition is inexplicably small and may represent a remarkably homogeneous subject group with little difference in endogenous circadian period.

DISCUSSION

The current results demonstrate that a single sequence of intermittent bright light exposure is effective in inducing phase delays of the human circadian timing system. The phase delays observed after intermittent bright light exposure were comparable to those measured after continuous bright light exposure ( ), even though the bright light represented only 23 of the 6.5 h light exposure session. Given The small change in phase observed in the control group may be attributed to the longer-than-24 h period of the circadian pacemaker (

), and the potentially the weak entraining effect of non-photic synchronizers ( ). We conclude that the additional 2.06 h average change 13 37 in phase (DLMOn) observed in the intermittent bright light group as compared to the control group represents the direct effect of the exposure to intermittent bright light on the circadian pacemaker.

Based on previous studies ( , , ) we expected a robust phase delay for both BL and IBL exposure. Indeed, the dynamic resetting significantly different from those after BL exposure, and phase delay shifts after IBL and BL exposures were significantly greater than those after VDL exposure.

The mathematical model developed by Kronauer and colleagues ( , , ) predicts that there is a greater activation rate (a measure 32 38 39 of resetting effectiveness of light) at the very beginning of a light pulse. As activator elements are depleted, the activation rate becomes limited by the rate at which elements are recycled. This model is consequently characterized by two rate parameters: a forward rate, , that α depends on light intensity, I, but is fixed at any given intensity, and a backward rate, , that is independent of intensity. For any given I, at β least two different temporal patterns are needed to estimate and from the response data ( , ). For I 10,000 lux, and were For a periodic stimulus pattern consisting of alternating light and dark episodes, we denote that the fraction of cycle time the light is " on as the duty cycle . If a window of time is available for stimulation, the maximum response will be generated by using light the entire " " " time (duty cycle 1). However, most of this light is used inefficiently since after a time, t , the incremental response is limited by the

= ≅ α 1 -
recycle rate and consequently is low. The maximally efficient use of a stimulus (i.e., maximum response per minute of stimulus) is found when the duty cycle is very low, but then the total response is also low. Intermediate duty cycles effect a compromise: less than maximum efficiency but total response not a great deal below maximum. The present experiment was designed with several objectives. The first was to expose subjects to a single sequence of bright light pulses, as opposed to multiple sequence of exposure over several days. A second was to extend the range of the duration of near darkness from the former range of 20 45 minutes to 60 minutes. A third was to study a short duty cycle of 0.23. Both these changes in intermittent light patterns, besides improving the accuracy of rate estimate, increase the practical usefulness of intermittent light intervention. Our fourth goal was to validate Kronauer s mathematical model of resetting efficacy of ' intermittent bright light ( , , ); simulations based on existing rate estimates indicated that the pattern of light exposure used in this 32 38 39 experiment would produce ~70 of the resetting drive onto the circadian pacemaker achieved by continuous bright light.

%

These results are in agreement with animal studies that reported that the resetting action of light was most efficient at the beginning of the light exposure, and that minimal additional phase shift was produced by further extension of the light stimulus ( , , , ) due to a 2 42 43 51 reduction in photic responsiveness of the circadian pacemaker ( ). 35

It is now recognized that the visual photopic and scotopic systems are not the only photosensitive structures in the retina. Recent studies in mammals have identified a subset of retinal ganglion cells (RGCs) that are intrinsically photosensitive and are involved in conveying photic input to the circadian system ( , , , , ). These cells contain a photopigment melanopsin, that appears to have a 4 26 28 44 46 broad role in non-visual photoreception ( ). Although the involvement of cryptochromes in phototransduction to the circadian system has 27 been suggested ( ), recent evidence suggest that these non-visual (melanopsin-containing RGCs) along with the classical visual (rod and 52 cones) photoreceptive systems account for all non-image forming functions involved in pupillary light reflex and circadian photoentrainment ( , , ).

Given the unusual light responses of the melanopsin-containing RGCs ( ), it is possible that these cells embody the Process L 4

proposed in the model of Kronauer et al. ( , ) that intervenes between light and the pacemaker oscillator (process P). Process L can 38 39 itself be parsed into two functions. One is a dynamical system in which activation elements are used and subsequently recycled. The other function is a simple power-law relation between the strength of light and the utilization rate parameter. The exponent for the power law inferred for human phase-shift data is proportional to the 0.6 power of light intensity. In their report, Berson ( ) include data on the et al 4 acute (1 sec average) intracellular depolarization of a ganglion cell as a function of irradiance at several wavelengths of light. A plausible representation of the irradiance dependence is a power law with an exponent of 0.6 over a range of 1.5 to 2 log units (from threshold to saturation). It is possible that these ganglion cells might also embody the dynamical system portion of process L, but the data presented thus far do not cover the temporal range required either to support or to reject this possibility. The response dynamics of rods and cones to light make them unlikely to embody process L. Further investigations, however, including exposures to intermittent light of different duty cycles, and different wavelengths are required to clarify this point.

These findings have important implication as they provide a greater understanding of the effects of brief episodes of light on the human circadian timing system, and suggest that the intermittent exposure to bright light that occurs in everyday life ( ) may have a 47 much greater impact on circadian entrainment than previously recognized. The present results also indicate that brief episodes of light exposure may be a cost-and time-effective way of resetting the circadian timing system in humans. Such strategies could be used not only

'

  technicians. The experimental suites were equipped with ceiling-mounted cool-white (4,100 K) fluorescent lamps Phillips (Eindhoven,[The Netherlands) T8 and T80 . A computer system automatically turned the lighting to the required pre-set intensity at the scheduled ] times. Maximum light intensities, as measured vertically at a height of ~1.83 m using an IL1400 photometer (International Light, Newburyport, MA), were ~190 lux during the waking hours of the three baseline days (~90 lux measured horizontally at a height of ~1.37 m); <8 lux (~1.5 lux measured horizontally at a height of ~1.37 m) during the constant routines and the day of light exposure; and

  >10.3 for women). Blood samples were collected in ethylenediaminetetraacetate (EDTA)-K tubes, rapidly centrifuged at 4 C, 2 °and the plasma was stored at or below 25 C until assayed. Samples were assayed for melatonin using radioimmunoassay techniques -°(Diagnostech/Pharmasan,

2 th 2 our

 22 results the DLMOn does indeed give consistently the lowest variance estimates. Czeisler et al. ( ) estimated the standard 25% 13

Fig. 4 %

 4 that the median phase delay measured with DLMOn was 3.03 h under continuous bright light condition (6.5 hours of bright light), 25% and 2.34 h in the intermittent bright light (1.5 h of bright light) condition, the phase resetting ability of intermittent bright light exposure, taking into account the reference 0.31 h phase shift measured in the control very dim light condition, was 3.2 times greater than that of continuous bright light exposure on a per minute basis of bright light exposure.

54 -7 10 31 Figure 1 Figure 2

 543112 Figure 1 Raster plot of the 10-day experimental protocol for subject 19a1. Scheduled sleep episodes (8 hours in darkness) are illustrated as black bars. During the baseline days (days 1 3), subjects were exposed to ~90 lux during wakefulness (16 hours). For the remainder of the study, except during the light exposure session, subjects were exposed to ~1.5 lux during wakefulness (gray bars). A ~ 26.2-h constant routine (CR1) was scheduled on days 4 5. The light exposure session (white bar) was scheduled on days 5 6, and consisted of 6.5 hours of exposure centered 5.8 -hours before habitual wake time (i.e., it started ~1.1 h before habitual bedtime and ended ~2.6 h before habitual waketime). Following the light exposure day, subjects underwent a 64-h CR (CR2) and were discharged after an ~22-h recovery sleep episode. Dim light melatonin onsets (DLMOn ) in CR1 and CR2 are indicated by . For this subject, the light exposure (BL condition) generated a phase delay of 3.15 h 25%

]

Figure 3

 3 Figure 3 Method used for melatonin phase assessment using Dim Light Melatonin Onset (DLMOn ), Dim Light Melatonin Offset (DLMOff ) 25% 25% and midpoint (MP ) phase markers. The thick line represents the three-harmonics fitted curve; the dashed line represents the amplitude of 25% the fitted curve; the dotted line represents 25 of the amplitude. The figure also illustrates the bright light exposure session (L.E.S.)-induced % melatonin suppression. The study days corresponding to the profile are given at the top of the figure. CR1: constant routine 1, CR2: constant routine 2.

Figure 4

 4 Figure 4 Individual resetting responses to the very dim light (0 ), intermittent bright light (23 ) and continuous bright light (100 ) exposure % % % conditions, measured with core body temperature minimum (CBT ) and melatonin phase markers: DLMOn (dim light melatonin onset), min 25% MP (midpoint), and DLMOff (dim light melatonin offset). Phase-advance shifts are plotted as positive values, whereas phase-delay 25% 25% shifts are plotted as negative values. Horizontal bars represent the median phase shift observed in each condition.

  -measure ANOVA (condition phase marker) with Greenhouse-Geisser correction (only original degrees of freedom are × reported), and Tukey tests for post-hoc multiple comparisons were used to compare phase shifts between phase markers, with H :The light stimulus was on average centered 3.15 0.72 h after DLMOn and 3.39 1.28 h prior to CBT (i.e., it started on Table 1 significant differences between the three groups in the initial phases of the temperature and melatonin rhythms at which the light intervention occurred. There were minimal changes in phase observed in the control group of subjects exposed to very dim light (median

								0 μ
	DLMOn	25% = μ	MP	25% = μ	DLMOff	25% = μ	CBT . Pearson correlations were calculated between phase shifts and initial phases of light min
	exposure, between CBT	min	and melatonin phases shifts, and between phase shifts of the melatonin phase markers (DLMOn	25%	and
	DLMOff ).
								±	25%	±	min
	average 5 minutes (0.09 0.72 h) before DLMOn ± min) phase change of 0.44 h for CBT and 0.31 h for DLMOn and ended 5 minutes (0.09 1.56 h) after CBT 25% ± ; ). A significantly greater phase delay was observed in subjects ( ). There were no + min -25% Fig. 4
	exposed either to continuous bright light median phase delay [	= -	2.70 h for CBT	min	(p<0.0002 vs. controls), 3.03 h for DLMOn	25%
	(p<0.0002 vs. controls), 2.71 h for MP -	25%	(p<0.0002 vs. controls), and 2.38 h for DLMOff -	25%	(p 0.0003 vs. controls) or intermittent = ]
	bright light median phase delay [	= -	1.73 h for CBT (p<0.0003 vs. controls), 2.34 h for DLMOn min -	25%	(p<0.0002 vs. controls), 2.30 h -
	for MP25 (p<0.0002 vs. controls), and 1.63 h for DLMOff % -	(p 0.001 vs. controls) . As expected for the selected phase of light

25%

RESULTS

25%

=

]

Table 1

 1 Timing of Light Exposure Sessions (L.E.S.) and individual phase shifts of the endogenous circadian rhythms of temperature (CBT ) and melatonin (DLMOn BL), intermittent light (IBL) and very dim light (VDL) exposure groups. By convention, phase delays are negative ( ) shifts, phase advances are positive shifts. Note that due to low amplitude of their CBT rhythms, phase shifts of CBT are not given for subjects 2008, and 2041, and the center of the light exposure sessions relative to CBT is not given for subject 2008. Due to missing blood min min samples, melatonin MP and DLMOn could not be calculated for subject 2056. Centers of the light exposure sessions and phase shifts are given in hours.

	min	25%	, MP	25%	, DLMOff	25%	) in the bright

± 4:28 1:40 ± 22:19 1:16 ± 22:42 0:46 ± 5:12 0:46 ± -2.51 1.22

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

[START_REF] Crowley | Combinations of bright light, scheduled dark, sunglasses, and melatonin to facilitate circadian entrainment to night shift work[END_REF]

Ackowledgements:

We thank the subject volunteers, the subject recruiters, the laboratory technicians involved in data collection, and Dr S.W. Lockley for his comments on an early version of the manuscript. This work described in this article was supported by grant NAG 5-3952 from the National Aeronautics and Space Agency, by, and by General Clinical Research Center grant GCRC-M01-RR02635 from the National Center for Research Resources.