A Mathematical Model of the Circadian Phase-Shifting Effects of Exogenous Melatonin

Model of Light’s Effects on the Circadian Pacemaker

The Kronauer et al. (1999) model of the effects of light on the human circadian pacemaker includes a light preprocessor in the retina, Process L, and the circadian pacemaker in the SCN, Process P. When light reaches the retina, it activates photoreceptors that project via the retinohypothalamic tract to SCN neurons.

In the model, photoreceptors exist as “activated” or “ready” (since transition time between states is negligible). Photoreceptors are activated at rate

α = α 0 ( I I 0 ) p I I + I 1 ,

where I₀ and I₁ are constants, and I represents environmental light intensity in lux. The value for p was empirically determined by Kronauer et al. (1999) and updated by St. Hilaire et al. (2007a, 2007b). Activated photoreceptors are converted back to the ready state at constant rate β (Kronauer et al., 1999; St. Hilaire et al., 2007a, 2007b). The fraction of activated photoreceptors, n, is modeled by

d n d t = α ( 1 − n ) − β n .

The resultant photic drive, B, to the circadian pacemaker is assumed to be proportional to the rate, α(1 – n), at which photoreceptors are activated:

B = G ( 1 − b x ) ( 1 − b x C ) α ( 1 − n ) .

In these equations, x represents circadian pacemaker activity, and x_C is a complementary variable. G is a constant that modifies photic drive strength, and (1 – bx)(1 – bx_C) modulates circadian phase dependence of light sensitivity through parameter b (Kronauer et al., 1999). The value b = 0.4 was chosen by Kronauer et al. (1999) to achieve a photic sensitivity maximum-to-minimum ratio of ~3:1 with maximum sensitivity approximately 3.8 h before the core body temperature nadir.

Two coupled differential equations model Process P (St. Hilaire et al., 2007b):

κ d x d t = x c + γ ( x 3 + 4 x 3 3 + 256 x 7 105 ) + B − η M ,

κ d x c d t = B x c 3 − x ( 24 f τ C ) 2 + K B − ζ M .

Parameter κ modifies photic drive B’s strength, τ _C represents intrinsic period, γ determines the oscillator’s stiffness (i.e., how quickly it returns to its limit cycle after perturbation) (Kronauer, 1990), and f is a correction term included to fix the period at precisely τ _C (Forger et al., 1999).

Constants η and ξ in Equations 4 and 5 modify the strength of the melatonin drive, M. Their values are experimentally determined (see the Results section) to achieve correct PRC shape.

Oscillator phase, φ, uses the standard polar definition

φ = tan − 1 ( x x C ) .

Endogenous and Exogenous Melatonin

The multicompartmental melatonin concentration model includes the pineal gland (H₁), blood plasma (H₂), and oral absorption site (H₃). The first 2 are adapted from St. Hilaire et al. (2007b), whereas absorption site is our addition. These compartments are modeled by the following differential equations:

d H 1 d t = − β I P H 1 + A ( ϕ ) ( 1 − m B ) S ( H 1 , B ) ,

d H 2 d t = β I P H 1 − β C P H 2 + β A P H 3 ,

d H 3 d t = − β A P H 3 .

During biological night, the SCN induces the pineal gland to produce melatonin, which infuses into plasma at rate β _IP . Melatonin is cleared from plasma by the kidneys at rate β _CP . Absorption of exogenous melatonin into plasma happens at rate β _AP ; this new parameter’s value is determined by experimental measures of time to peak (details below). Parameter values for a, β _IP , and β _CP , used in Equations 7, 8, 10, and 11, were previously estimated by St. Hilaire et al. (2007b). Minor updates to these parameters’ values were derived for our new model as described below.

A(φ) in Equation 7 represents pineal melatonin production and is defined below. (1 – mB)S(H₁,B) models melatonin suppression by light input, B, where m is a constant modifying B’s strength, originally derived by St. Hilaire et al. (2007b). S is a step function that prevents H₁ from falling below 0 in response to strong light stimuli, since a negative concentration is unphysiological. When (1 – mB) ≥ 0, S = 1, but when H₁ < 0.001 pmol/L (threshold chosen to approximate 0) and (1 – mB) < 0, S = 0.

At the time of a dose delivery, the value of H₃ is increased by an amount proportional to the exogenous dose at a rate determined by the pharmacokinetic parameters in Equations 7, 8, and 9. To convert dose mass to an effective blood concentration, we assume equal distribution across 5 L of blood volume (Nadler et al., 1962), the accepted “average” adult’s blood volume. To account for dose loss by first-pass hepatic metabolism (Claustrat et al., 2005), we multiply the effective blood concentration by a parameter 0 ≤ ε ≤ 1, which represents the fraction of the dose that is absorbed into the bloodstream. This parameter may be different for different dose preparations. An alternative approach would be to explicitly model first-pass metabolism by adding an additional compartment; this is not necessary for the results presented here.

Serotonin to melatonin conversion is catalyzed by arylalkylamine N-acetyltransferase (AA-NAT). This enzyme-catalyzed melatonin precursor is then O-methylated by hydroxyindole O-methyltransferase (HIOMT), producing melatonin. This circadian phase dependent process is modeled by A(φ) (St. Hilaire et al., 2007b):

A ( ϕ ) = a 1 − e − δ M 2 π ( ϕ on − ϕ ) 1 − e − δ M 2 π ( ϕ on − ϕ off ) , for ϕ > ϕ off and ϕ < ϕ on ,

A ( ϕ ) = a e − r M 2 π ( ϕ on − ϕ off ) , f o r ϕ on < ϕ < ϕ off .

The phases φ_on and φ_off are those at which melatonin synthesis is turned on and off. The function M_2π(x) returns x modulo 2π (cycle period). The constant r converts from time to phase, for consistency with previous work in which melatonin onset and offset times were used rather than phases (St. Hilaire et al., 2007b).

The melatonin drive to the SCN represents input to melatonin receptors from both endogenous and exogenous melatonin blood plasma. It is modeled by the sigmoid function,

M = M max 1 + e H s a t − H 2 σ ,

where M decreases to a negligible value as H₂ decreases and approaches M_max as H₂ increases. The midpoint of this rising portion of the function occurs at H₂ = H_sat and M = M_max/2. The curve becomes steeper for smaller values of σ.

A sigmoid function was chosen to simulate receptor saturation. The value of H_sat was estimated from dose-responses of several studies (Zhdanova, 2005; Dollins et al., 1994; Stone et al., 2000; Burgess et al., 2010). Stone et al. (2000) compared the effects of 0.1-, 0.5-, 1.0-, 5.0-, and 10.0-mg doses on sleep onset and observed no dose-dependence over the range 0.5 to 10.0 mg. Studies by Dollins et al. (1994) and Zhdanova et al. (1995) found that doses between 0.3 and 10.0 mg produce effects of similar magnitude. Dollins et al. (1994) found that the physiological concentration of melatonin during the biological night was roughly 200 pg/mL, a concentration that could be duplicated during the biological day by administering a dose between 0.1 and 0.3 mg. Zhdanova et al. (1995) found that above a 0.3-mg dose there was no significant increase in response. Together, these findings suggest that melatonin receptors are saturated at levels above ~200 pg/mL. These findings are consistent with a more recent study (Burgess et al., 2010) in which only minimal dose-dependence in phase shifting was observed between 0.5 mg and 3.0 mg. We therefore take 200 pg/mL (861 pmol/L) to be the value of H_sat. Other studies have also reported saturation of melatonin receptors; results from Deacon and Arendt (1995), who looked at acute effects on core body temperature, suggest that the half-maximal response occurs between 0.05 mg and 0.5 mg, which after digestion and metabolism produce plasma concentrations of ~430 pmol/L and ~4300 pmol/L, respectively, consistent with the value of H_sat chosen here.

After fixing a nominal value for H_sat, we constrained σ based on dose responses observed experimentally. Burgess et al. (2010) derived phase response curves for 3 consecutive daily doses of 0.5 mg or 3.0 mg, finding a 4% increase in PRC amplitude for the higher dose. These findings place an upper bound of ~10³ for σ. Larger values of σ would result in a flatter sigmoid curve, and therefore a more gradual dose response, which would not be consistent with experimental data. Therefore, a nominal value of σ = 50 pmol/L was chosen. To achieve tighter constraints on σ, PRC data are needed from doses well below H_sat. Those PRCs have not been produced, however, and doses used clinically (including those modeled below) are large enough to saturate melatonin receptors. At these concentrations, the model is insensitive to the precise value of σ.

The values of M_max and other parameters were determined by fitting our model to the PRC of Burgess et al. (2008) for a 3.0-mg dose of exogenous melatonin, described below in the “PRC Methods” section

Pharmacokinetics Methods

Melatonin’s pharmacokinetics are determined by 3 parameters: ε (preabsorption metabolism), β _CP (effective clearance rate), and β _AP (effective absorption rate). To fit values to β _AP and β _CP , we compared the model to blood plasma melatonin concentration data from Wyatt et al. (2006). Specifically, we required that the model accurately reproduce both time to peak (i.e., time from ingestion to peak concentration) and physiological half-life, for both 0.3-mg and 5.0-mg doses. Their data showed a half-life of ~1.0 h and a time to peak of ~1.2 h. Consistent with this, Waldhauser et al. (1984) found the half-life of melatonin to be ~50 min. Wyatt et al. (2006) analyzed melatonin’s pharmacokinetic profile at 2 different circadian phases, one in which endogenous melatonin was present and one in which endogenous melatonin was absent. Although concentrations were slightly higher during the exogenous melatonin-present trials, the time to peak and clearance rate were not affected by circadian phase. Therefore, we did not model difference in pharmacokinetics with circadian phase. Parameter ε modifies the amplitude of melatonin’s concentration profile over time. Its value was set to 0.1 based on the findings of Claustrat et al. (2005) and confirmed after fitting values for β _AP and β _CP to the Wyatt et al. (2006) data. The value of ε and the values of β _CP and β _AP were used to simulate the Burgess et al. (2008, 2010) protocols below.

Exogenous melatonin pharmacokinetics are preparation dependent; therefore, a single parameter set will never fit all data. Aldhous et al. (1985) showed that for an equal dose of 2 mg in gelatin capsules, corn oil, and a slow-release pill, different plasma concentration profiles resulted, presumably because of different degrees of preabsorption metabolism and rates of absorption. Comparison of gelatin preparations used by Aldhous et al. (1985) and Wyatt et al. (2006) reveals further differences. Wyatt et al. (2006) used a gelatin capsule in which 5.0 mg of melatonin was suspended in cellulose, whereas Aldhous et al. (1985) used a gelatin capsule in which 2.0 mg of melatonin was suspended in lactose. The 2 studies show different times to peak and half-lives. These differences cannot be accounted for by changing dose size and instead must be accounted for by changing pharmacokinetic parameters, because the model’s solution to a single dose is a biexponential equation, in which time to peak and half-life are independent of dose size. To illustrate the model’s ability to reproduce different plasma melatonin profiles, we modeled a 2-mg dose for both the Aldhous and Wyatt gelatin preparations by fitting values of ε and β _AP for each.

PRC Methods

We simulated the Burgess et al. (2008) PRC protocol to fit the model’s sensitivity to exogenous melatonin. The protocol began with 7 days of habitual sleep schedules, followed by 5 days on a 4.0-h ultradian sleep-wake and light-dark cycle alternating between 1.5 h of permitted sleep with lights off and 2.5 h kept awake in room light (>150 lux). The ultradian cycle minimized phase-shifting effects of light by evenly distributing light exposure across circadian phases. For the first three 24-h days of this ultradian cycle, participants were given a 3.0-mg gelatin-based capsule of oral melatonin (Burgess et al., 2008) or a placebo at the start of a wake episode at the same clock time each 24-h day. A PRC was constructed based on the last day’s phase assessment. Values of parameters M_max, η, and ξ in Equations 4, 5, and 12, which together determine the strength of the phase resetting effect, were explored over a range of 1 order of magnitude above and below final values and optimized using a least-square data-fitting method to fit the 3.0 mg PRC. The 3-dimensional parameter space spanned by M_max, η, and ξ was explored using a fixed grid loop algorithm written in MatLab 7.11.0 (The MathWorks, Natick, MA) that simulated a PRC at each coordinate of the parameter space. Residuals were calculated as the vertical distance between experimental data points and the simulated PRC. The final set of model parameters is listed in Table 2.

OPEN IN VIEWER

Table 2.

Model parameters.

Parameter	Value and Units (Previous Value)	Parameter	Value and Units (Previous Value)
I0	9500 lux	η	0.04
	(9500 lux)
I1	100 lux	ζ	0.54
α0	0.1 sec	β IP	7.83 × 10−4 sec−1
	(0.05 sec)		(7.83 × 10−4 sec−1)
β	0.007 sec	β CP	3.35 × 10−4 sec−1
	(0.0075 sec)		(4.17 × 10−4 sec−1)
p	0.5	β AP	1.62 × 10−4 sec−1
r	15.36 sec	a	1.0442 × 10−3 sec−1
G	37	ϕon	6.113
b	0.4	ϕoff	4.352
γ	0.13	δ	600 sec
κ	12/π h	M max	0.019513
τ c	24.1 h	H sat	861 pmol/L
f	0.99729	σ	50 pmol/L
m	7 sec−1	ε	0.10

Best fit parameter values are given in the order that they appear in the Methods section. Previously used parameter values from St. Hilaire et al. (2007b) are given when applicable, even if the values are unchanged.

Using these parameter values, we simulated 2 independent data sets. First, we simulated an experiment by Burgess et al. (2010), which is identical to the 2008 experiment except for the use of 0.5-mg instead of 3.0-mg doses. Then, to investigate the effects of a light-dark cycle on the melatonin PRC, we simulated the Lewy et al. (1992, 1998) protocol in which baseline DLMO was assessed on day 1, followed by a week of protocol acclimation with placebo. On day 7, pretreatment DLMO was obtained followed by 2 days of placebo (days 8 and 9), 4 days of 0.5 mg melatonin (days 10-13), and posttreatment DLMO measurement on day 14. With the exception of DLMO measurement days (1, 7, and 14), subjects were at home in uncontrolled and unrecorded light exposure environments. For individuals previously adapted to dim light (e.g., in inpatient experiments), the dynamic range of the circadian phase shifting response to light falls between 0 and 1000 lux (Zeitzer et al., 2000); we therefore simulate light intensities from this range. Light-dark cycles simulated had 16 h of 5, 50, or 500 lux with 8 h of 0 lux in each simulated 24-h day. The parameter set was estimated using the least-squares data fitting method described above.

Pharmacokinetics Results

Simulated and experimental (Wyatt et al., 2006) pharmacokinetics results for 0.3-mg and 5.0-mg data are shown in Figure 2. For the gelatin capsule preparation modeled here, we use ε = 0.1, β _AP = 1.62 × 10⁻⁴ sec⁻¹, and β _CP = 3.35 × 10⁻⁴ sec⁻¹ to achieve a 100-min time to peak and 100-min half-life with peak melatonin concentration following a 0.3-mg and 5.0-mg dose of 5 × 10² and 10⁴ pmol/L, respectively. Simulations of the Wyatt et al. (2006) and the Aldhous et al. (1985) oral doses of gelatin capsule melatonin are shown in Figure 3. The Aldhous et al. preparation has a faster time to peak, faster clearance rate, and greater preabsorption metabolism than the Wyatt et al. preparation, with ε = 0.03, β _AP = 6.67 × 10⁻⁴ sec⁻¹, and β _CP = 3.35 × 10⁻⁴ sec⁻¹, representing a time to peak of ~50 min, a half-life of ~40 min, and a peak melatonin concentration following a 2.0-mg dose of 1.15 × 10³ pmol/L. We assessed model sensitivity to these parameters. Increasing β _CP by 20% decreased time to peak by 10% and half-life by 18%. Increasing β _AP by 20% decreased time to peak by 7% and half-life by 9%.

OPEN IN VIEWER

Figure 2.

Experimental (symbols, from Wyatt et al. 2006) and modeled (lines) plasma melatonin concentrations over 7 h following 5.0-mg and 0.3-mg doses of oral melatonin.

OPEN IN VIEWER

Figure 3.

Simulated pharmacokinetic profiles following the Aldhous (1985) (dashed line) and Wyatt (2006) (solid line) 2.0-mg gelatin-capsule melatonin doses. Experimental data from Aldhous et al. (1985) indicated by symbols.

PRC Results

To verify that a phase response is indeed present, we simulated zero phase response (i.e., a horizontal line) at all phases. The adjusted R² value was 0.616 when phase shifts were simulated using the model and only 0.321 when zero phase shifts were simulated. We therefore proceeded with the model for all simulations. We note that the uncertainty in the original data is not reported, thus limiting our data fitting comparison to that of the residual sums for different simulations.

Simulated and experimental PRC data for the 3.0-mg dose are shown in Figure 4. Exogenous melatonin administered around the time of DLMO causes negligible phase resetting. Doses given prior to DLMO (left of 0 h in the figure) cause phase advances, whereas doses given after DLMO (right of 0 h in the figure) cause phase delays. The model thus predicts type 1 resetting to 3.0 mg of exogenous melatonin.

OPEN IN VIEWER

Figure 4.

Melatonin PRC for 3 consecutive daily doses of 3.0 mg of oral melatonin. Data extracted from Burgess et al. (2008) are represented as symbols; model output is represented by the solid curve. The vertical line marks the stimulus time required for maximum phase advance (horizontal line). The thin, solid line across the figure demarks a phase shift of 0 h.

In a more recent experiment, Burgess et al. (2010) showed that compared with a 3.0-mg dose, a 0.5-mg dose of the same preparation results in a slightly lower peak that occurs at a slightly later phase. The model correctly predicted these findings without modification to model parameters (Figure 5).

OPEN IN VIEWER

Figure 5.

Simulated PRCs for 3 consecutive daily doses of 3.0 mg (solid line) and 0.5 mg (dashed line) of melatonin. The vertical lines mark timing required for the peak phase advance and phase delay (horizontal lines) for each PRC. The horizontal straight, dashed line at zero shows zero phase shift.

Since Lewy et al. (1998) did not control light exposure, and light is a stronger phase-shifting stimulus than is melatonin, the simulations of that protocol were performed with 4 different light levels during the wake episodes (Figure 6). As light intensity decreases, the PRC peak occurs earlier and at higher amplitude, similar to the effects observed when increasing exogenous melatonin dose (Figure 5). Higher light intensities reduce melatonin-induced phase shifts due to light’s entraining effects, which may lock the rhythm to a stable phase.

OPEN IN VIEWER

Figure 6.

Simulated PRC for 4 consecutive doses of 0.5 mg of exogenous melatonin. Symbols represent experimental data (Lewy et al. 1998). The 4 curves show results from constant darkness (DD; 0 lux) and 3 light levels (5, 50, and 500 lux) within a light-dark (LD) cycle of 16 h light and 8 h dark. Note that PRCs for DD and 5 lux conditions overlap. The horizontal straight, solid line at zero shows zero phase shift.