Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Indirect Effect Model with Surge Function for Describing Melatonin Circadian Rhythm: Quantitative Comparison and Application between Normal Sleepers and Patients with Delayed Sleep-Wake Phase Disorder

Indirect Effect Model with Surge Function for Describing Melatonin Circadian Rhythm: Quantitative... IntroductionDelayed sleep-wake phase disorder (DSWPD) is a circadian rhythm disorder characterized by delayed sleep phases that cause difficulty falling asleep or waking up on time, leading to daytime dysfunction [1]. The prevalence of DSWPD in adults ranges from 1.51% to 8.90% [2]. In individuals with DSWPD, the endogenous melatonin rhythm is delayed and no longer aligns with desired sleep times [3]. Patients with this disorder are usually unable to fall asleep at the desired bedtime, which can result in sleep deprivation, daytime sleepiness, and fatigue, which can have a serious impact on the patient’s life [4, 5]. Consequently, it is of paramount importance to implement timely intervention and treatment for DSWPD.The principal therapeutic modalities for DSWPD encompass light therapy, cognitive-behavioral therapy, and melatonin treatment, all of which are chronobiological in nature [6‒8]. Melatonin is regarded as the optimal pharmacological intervention for the treatment of DSWPD in patients who have not responded to non-pharmacological treatments [3]. Exogenous melatonin can alter the timing of the biological clock, and depending on the time of administration, melatonin can advance or delay the phase of the biological clock [9]. The advancement of the phase occurs when melatonin is administered between 2 and 7 h prior to the dim light melatonin onset (DLMO). Conversely, the phase is delayed when melatonin is administered in the early morning [10]. For patients with DSWPD, the goal of treatment is to advance the circadian time so that patients can sleep at a more regular time. Due to the individual differences in the circadian rhythm of melatonin, treatment with melatonin may not only be ineffective, but may even be counterproductive or cause adverse effects such as drowsiness and disruption of the biological clock if it is not administered in accordance with the individual patient’s circadian rhythm, which is usually assessed using a DLMO [3, 11]. Therefore, selecting the appropriate time of administration based on the circadian rhythm characteristics of melatonin is crucial for individualized treatment of DSWPD.Melatonin is a neuroendocrine hormone primarily synthesized in the pineal gland, though it is also produced in peripheral organs such as the skin, heart, liver, uterus, placenta, kidneys, immune cells, and thymus. However, only the secretion from the pineal gland and the visual system is regulated by the environmental light/dark cycle, as other tissues lack this regulatory mechanism [12, 13]. Melatonin secretion is regulated by circadian rhythms, usually increasing at night and decreasing during the day, which helps regulate sleep and wake cycles [14]. By comparing the circadian rhythm characteristics of DSWPD patients and normal sleepers, the etiology of DSWPD patients can be more accurately assessed, thereby providing more references for the development of individualized treatment plans. Micic et al. [15] observed that the average melatonin concentration-time curve is steeper in the ascending phase in normal sleepers than in DSWPD patients, suggesting that the reduced rate of melatonin production may contribute to the pathogenesis of DSWPD patients. Shibui et al. [16] discovered that the melatonin circadian rhythm in DSWPD patients is significantly delayed compared to normal sleepers by analyzing their average melatonin concentration-time curves.Population pharmacokinetic/pharmacodynamic (PK/PD) modeling is a powerful tool for characterizing the variability in drug exposure and response across individuals, integrating both PK (drug concentration) and PD (drug effect) data to describe the dose-response relationship [17, 18]. This approach is particularly useful for understanding inter-individual variability and optimizing dosing regimens in diverse populations. In contrast, indirect effect models focus primarily on PD effects, often used to describe the time-dependent changes in endogenous substances or biomarkers without directly modeling drug concentrations [19].Indirect PD models are particularly valuable for characterizing endogenous substances, as they account for the production and elimination processes of these compounds. Nagaraja et al. [20] utilized an indirect effect model including a surge function to characterize the delay in the luteinizing hormone surge following single or multiple administrations of cetrorelix. Similarly, an indirect effect model with a surge function has been successfully employed to describe the circadian rhythms of ACTH, cortisol, and melatonin in healthy males [21, 22].In this study, we aimed to apply an indirect effect model combined with a surge function to characterize the properties of endogenous melatonin in normal sleepers and DSWPD patients. By leveraging this approach, we seek to elucidate the circadian rhythm of melatonin production and elimination, providing insights into its role in sleep regulation and potential therapeutic applications.MethodsData SourceThe modeling dataset comprised data from 36 subjects, including data from 10 normal sleepers and 26 DSWPD patients. Additional data from three normal sleepers and 3 DSWPD patients were used for model application (Fig. 1).Fig. 1.Flowchart of model building process.Modeling data for 10 normal sleepers (mean age 33.60 ± 11.22 years) were obtained from a clinical study (registration number: NCT06467851). Subjects were admitted to the study ward before 16:00, and saliva samples were collected between 16:00 and 11:00 the following day. Samples were collected at 1-h intervals between 16:00 and 7:00 the following day, and at 2-h intervals between 7:00 and 11:00. Saliva samples were collected in the study ward in dim light. Furthermore, the use of light-producing electronic devices, such as mobile phones, was prohibited. Products that could affect salivary melatonin concentrations, such as alcohol, coffee, and strong tea, are prohibited throughout the study period. Melatonin concentrations were determined using the salivary melatonin ELISA (IBL International GmbH, Hamburg, Germany, RE54041). In daytime samples (<8 pg/mL), the mean coefficients of variation within and between groups were 17.0% and 20.5%, respectively. In nighttime samples (>10 pg/mL), the mean coefficients of variation within and between groups were 13.9% and 18.4%, respectively.Modeling data for 26 DSWPD patients were digitally extracted from individual patient profiles reported by Micic et al. [15]. The original study by Micic et al. [15] met the following screening criteria: (1) patients with DSWPD were included; (2) patients were not treated with any medication that affects endogenous melatonin levels or with phototherapy interventions; (3) salivary melatonin was measured as the study endpoint; (4) the characteristics of a complete 24-h melatonin secretion were assessed; (5) subjects were sampled in dim light; (6) strictly control the intake of products that affect salivary melatonin level, such as alcohol, coffee, and cigarettes.Data collection for both normal sleepers and DSWPD patients was confined to autumn and winter, limiting the assessment of seasonal effects on the model. Complementing our findings, a comprehensive study by Titone et al. [23] examined salivary melatonin circadian rhythms across diverse racial groups, including Caucasians, African Americans, Asians, Pacific Islanders, and other ethnicities (e.g., individuals of Mexican, Native American, and Hispanic descent). Their results revealed no statistically significant racial differences in melatonin rhythms. Consequently, racial variation was not incorporated into the final model in this study. The data used for model application, employing the Bayesian estimation method, included data from three normal sleepers sourced from Burgess et al. [24] and data from 3 DSWPD patients sourced from Burgess et al. [25].Model DevelopmentThis study describes the circadian rhythm of endogenous melatonin using an indirect effect model that includes surge functions. Based on the final model, Monte Carlo simulations were conducted to quantitatively analyze the circadian rhythm characteristics of normal sleepers and patients with DSWPD. Lastly, Bayesian feedback was used to predict DLMO and melatonin concentration at different time points based on sparse data. The process of model building and application is shown in Figure 1.Structural ModelIn normal sleepers, melatonin synthesis begins at night, between 20:00 p.m. and 22:00 p.m., and reaches a maximum at midnight or between 2:00 a.m. and 3:00 a.m. After that, melatonin production gradually decreases and remains very low during the day [26]. Melatonin production and elimination can be described by an indirect effect model, as shown in formula (1) [20, 22, 27]:(1)dRdt=kin−kout×Rwhere R is the concentration of melatonin, Kin is the zero-order release rate constant of melatonin, and Kout is the first-order elimination rate constant of melatonin. Melatonin is a hormone secreted by the pineal gland. In mammals, the pineal gland is inhibited by light during the day and activated at night through a polysynaptic pathway in the SCN, which promotes melatonin secretion [28]. Therefore, the surge in melatonin at night may be attributed to an increase in its production rate. Charles et al. [22] have successfully applied a surge function to characterize the melatonin surge at night, so we chose an indirect effect model that included a surge function to represent the changes in melatonin over time (shown in Fig. 2), as shown in Equation (2-4) [20]:(2)dRdt=kin×Surget−Kout×R(3)Surget=1+AMPt−T0WIDN+1(4)Kin=Kout×BaselineIn Equations (2-4), Surge(t) denotes the surge function; Baseline refers the baseline level of melatonin; AMP indicates the amplitude of the peak, which is a dimensionless parameter related to the amplitude of the peak; t represents time; T0 is the time to peak; WID denotes the width of the peak; and N signifies the exponential of the slope of the peak, where N is an even number (e.g., 2, 4, 6, etc.) to ensure that the melatonin concentration remains positive both before and after the peak occurs. If t ≫ T0 or t ≪ T0, melatonin concentration is close to baseline, and if t ⟶ T0, there is a sharp increase in melatonin concentration.Fig. 2.An indirect effect model with surge function. R represents the PD response variable. Kin(t) denotes periodic production rates. Kout denotes the elimination rate. Dashed red arrow indicates stimulation. The surge function illustrated in this figure is adapted from the study by Olta Tarko [29].Random Effects ModelThe nonlinear mixed-effects model was established using Phoenix NLME (version 8.3.5, Certara, Co., Princeton, NJ, USA). The model parameters were estimated using first-order conditional estimation-extended least squares. The inter-individual variation (IIV) was expressed using an exponential model (Equation 5):(5)Pij=Pj×eηijPj represents the typical value of the jth parameter in the population, while Pij denotes the true value of the jth parameter for the ith subject. ηij denotes the IIV in the jth parameter for the ith subject. The IIV (η) conforms to a normal distribution with mean 0 and variance ω2. The residual variances were expressed using additive, proportional, and combined additive and proportional models (Equations 6, 7, and 8):(6)Cij=IPREDij+εij(7)Cij=IPREDij×1+εij(8)Cij=IPREDij×1+εij,1+εij,2Cij represents the observed concentration of the ith subject at the jth sampling point, whereas IPREDij denotes the model-predicted value for the subject. εij denotes the residual variance at the jth sampling point for the ith subject. The residual variance (ε) conforms to a normal distribution with mean 0 and variance σ2.Covariate and Covariance ModelsThe model examined the effect of disease state on melatonin circadian rhythms, so disease state was included in the model as a categorical variable. Assess the correlation between the IIV distributions of T0, Kout, AMP, Baseline, and WID through graphical analysis, and examine this relationship by incorporating the covariance between IIV parameters within the model.Model EvaluationThe degree of fit of the model to the observed data and the accuracy of the model were evaluated by diagnostic graphs (e.g., observations vs. population predictions [PRED], observations vs. individual predictions [IPRED], conditional weighted residuals [CWRES] vs. PRED, and CWRES vs. time). The predictive power of the final model was evaluated using the visual predictive check. A total of 1,000 simulations were performed for the final model, calculating the 5%, 50%, and 95% percentiles of both the observed data set and the simulated data set at various time points to compare their distribution characteristics. The bootstrap resampling method was employed 1,000 times to calculate the median values and 95% confidence intervals (2.5–97.5%) of the parameters, which were then compared to the original parameter estimates from the model to assess the robustness of the final model.SimulationMonte Carlo simulation (n = 1,000) was performed based on the final model, comparing the circadian rhythm characteristics of melatonin in normal sleepers and DSWPD patients, and DLMO was calculated from the simulation results. The absolute threshold method was used to calculate DLMO. DLMO is considered as the time when the salivary melatonin concentration reaches a threshold of 7 pg/mL [30].The external data used for Bayesian estimation consisted of data from three normal sleepers and 3 DSWPD patients, which were obtained from the references Burgess et al. [24] and Burgess et al. [25], respectively. The sparse data of three normal sleepers and 3 DSWPD patients were imported into Phoenix, respectively. Based on the established final model, the maximum a posteriori Bayesian method was employed to estimate melatonin concentrations at different time points and DLMO.ResultsSubjectsThe study included 1,740 concentrations from 36 subjects. A total of 180 concentrations were obtained from 10 normal sleepers, while 1,560 concentrations were obtained from 26 patients with DSWPD. The demographic characteristics and sleep quality scores of the 36 subjects are summarized in Table 1. The 10 normal sleepers consisted of 4 males and 6 females with a mean age of 33.6 ± 11.22 years, mean weight of 66.43 ± 7.07 kg, and a mean BMI of 24.57 ± 3.16 kg/m2. The DSWPD patients included 18 males and 8 females with an average age of 21.73 ± 4.98 years [15].Table 1.Summary of demographic characteristics and sleep quality scoresNormal sleeper (N = 10)DSWPD patients (N = 26)aMale, n (%)4 (40%)18 (69%)Female, n (%)6 (60%)8 (31%)Age, years33.6±11.2221.73±4.98Weight, kg66.43±7.07NABMI, kg/m224.57±3.16NAMEQ score49.4±7.2933.09±5.76PSQI score3.4±1.906.15±3.32Values are expressed as means ± SD.BMI, body mass index; NA, not applicable.aThe demographic characteristics and sleep quality scores of DSWPD patients were obtained from Micic et al. [18].Melatonin Population ModelAn indirect effect model incorporating a surge function was employed as the structural model, which effectively characterized the circadian rhythm distribution property of endogenous melatonin. Compared to additive models and additive-proportional error models, the proportional error models provided a better description of residual variance. An exponential error model was used to describe IIV. The inter-individual variants were included in the T0, Kout, AMP, Baseline, and WID. The final model estimation results are shown in Table 2. The model parameters were estimated with high precision, and their RES% values were low. Typical values of T0, Kout, AMP, Baseline, and WID estimated by the final model were 23:59, 1.23 h−1, 7.95, 3.21 pg/mL, and 4.12 respectively. The proportional-type residual variance of the model estimate was 0.21.Table 2.Parameter estimates of the final model and bootstrap evaluationParameterFinal modelBootstrapestimateRES, %95% CImedian95% CItv AMP7.9515.475.54∼10.368.005.79∼10.86tv T0, hh:mm23:594.1322:37∼01:2200:0222:44∼01:23tv WID4.125.783.65∼4.584.163.63∼4.52tv Kout, h−11.2321.820.70∼1.751.210.80∼1.93tv Baseline, pg/mL3.2123.271.75∼4.683.221.77∼6.45Inter-individual variabilityω2 T00.0119.950.006∼0.010.010.00∼0.01ω2 Kout1.0133.020.36∼1.660.910.32∼1.63ω2 AMP0.8526.670.41∼1.300.820.38∼1.28ω2 Baseline0.2153.38−0.01∼0.430.18−0.02∼0.42ω2 WID0.0672.33−0.03∼0.150.06−0.02∼0.16Residual variability (σ)Stev0a0.216.720.18∼0.240.210.18∼0.24AMP, amplitude; T0, time of peak melatonin concentration; WID, a peak width parameter; Kout, first-order elimination rate constant of melatonin; Baseline, baseline melatonin concentration.aStev0 represents residual variability (or intra-individual variability).The disease status was included as a categorical variable affecting T0 and Baseline in the model, resulting in a reduction of −2LL by approximately 146 (3,769 vs. 3,623). The results of the graphical evaluation indicate a correlation among the IIVs of AMP, Baseline, and WID. Incorporating their covariance structure enhanced model fit, reducing −2LL by 17 (3,623 vs. 3,606). Therefore, the disease status was included as a categorical variable in the final model, and the covariances among AMP, Baseline, and WID were examined.Model EvaluationThe final model was diagnosed by goodness-of-fit plots, including observations versus PRED, observations versus IPRED, CWRES versus PRED, and CWRES versus time (shown in Fig. 3). Although the trend lines of the observations versus PRED scatterplot deviate slightly from the reference line, the trend line of the observations versus IPRED scatterplots is highly coincident with the reference line, indicating that the observations are well fitted to the model predictions. The CWRES was distributed symmetrically on both sides of the reference line (y = 0) and did not exhibit a significant trend change with group predicted values and time. This indicates that the model was unbiased and well fitted. The results of the visual predictive check indicate that observations at the 5th, 50th, and 95th percentiles exhibit similar distributional characteristics to those of the modeled data (shown in Fig. 4). The results of the nonparametric bootstrap are shown in Table 2. The parameter values estimated by the final model are close to the median of the bootstrap results and are within the bootstrap 95% confidence interval, indicating that the final model is stable.Fig. 3.Goodness-of-fit plots for the final model. a Observations (DV) versus PRED. b DV versus PRED logarithm. c DV versus IPRED. d DV versus IPRED logarithm. e CWRES versus time. f CWRES versus PRED.Fig. 4.VPC plot of the final model. The blue lines represent the 5th, 50th, and 95th percentiles of observed melatonin concentration over time. The black lines represent the 5th, 50th, and 95th percentiles of the predicted value over time. The shaded region describes the 95% prediction interval of the 5th, 50th, and 95th percentiles for predicting melatonin concentration (n = 1,000 simulations). VPC, visual predictive check.SimulationAfter adding disease state as a covariate to the model, the −2LL of the model decreased (3,769 vs. 3,623), suggesting that disease state affects the circadian rhythm of melatonin. Therefore, the resampling of normal sleepers and DSWPD patients produced two virtual groups. Figure 5 illustrates that endogenous melatonin secretion in both normal sleepers and DSWPD patients exhibits a distinctive circadian rhythm feature. The area under the curve (AUC) in normal sleepers was significantly greater than the AUC in patients with DSWPD, which was 3.1 times greater than the AUC in patients with DSWPD. The results of the final model simulation, as presented in Table 3, indicated that patients with DSWPD had a delayed T0 of approximately 4.06 h, their Baseline was approximately 2.44 pg/mL lower than that of normal sleepers, and the AMP and WID were similar to that of normal sleepers. Additionally, the DLMO for DSWPD patients was delayed by approximately 7 h compared to normal sleepers. The rate of melatonin elimination in DSWPD patients was similar to that of normal sleepers, but their production rate constant was significantly lower. One potential explanation for the delayed melatonin circadian rhythm observed in patients with DSWPD is a reduced rate of melatonin production.Fig. 5.Average melatonin concentration-time curve after 1,000 simulations of the final model. Red line represents normal sleeping subjects, and blue line represents DSWPD patients.Table 3.The parameter values outputted after 1,000 simulations of the final modelParametersNormal sleepersDSWPD patientsAMP12.13±14.1012.20±14.07T0, hh:mm00:05±1.6904:08±2.11WID4.26±1.104.25±1.09Baseline, pg/mL3.58±1.741.14±0.56Kout, h−12.04±2.582.03±2.73Kin, pg/h7.35±10.912.33±3.64t1/2, ha0.93±1.200.94±1.30DLMO, hh:mmb17:4200:54AUC, pg/mL*hb422.30134.50Values for characteristics are presented as mean ± SD unless otherwise indicated.AMP, amplitude; T0, time to peak; WID, the peak width; Baseline, baseline levels of melatonin; Kout, first-order elimination rate constant of melatonin; Kin, zero-order release rate constant of melatonin; t1/2, half-life; DLMO, dim light melatonin onset; AUC, area under the melatonin concentration-time curve.aThe half-life is calculated based on the Kout, using the formula: t1/2 = 0.693/Kout.bDLMO and AUC were calculated from the average melatonin concentration-time curves generated using simulated data for normal sleepers and DSWPD patients.Sparse data from 3 normal sleepers and 3 DSWPD patients were imported into Phoenix, respectively. Based on the established final model, the maximum posteriori Bayesian method was used to estimate melatonin concentration at different time points and DLMO. The final results indicated that the DLMO estimated by the model for DSWPD patients and normal sleepers was similar to the DLMO calculated from the original data, with an error margin within ±10% (Table 4). Additionally, the circadian rhythm characteristics of melatonin for both DSWPD patients and normal sleepers derived from the model show a high degree of similarity to the literature data, with a significant overlap in the melatonin concentration-time profiles (shown in Fig. 6).Table 4.The DLMO of the final model estimate and the DLMO of the original data calculationClassSubjectsDLMOapredicted data, hh:mmObserved data, hh:mmPEbDSWPD patients101:2601:28−0.11%201:2801:250.33%322:0221:570.54%Normal sleepers423:4723:390.78%523:4523:380.72%622:1321:265.47%aDLMO denotes the dim light melatonin onset.bPE denotes the prediction error, which is given by the formula: PE = [(prediction data−observed data)/(observed data)] × 100%.Fig. 6.Comparison of simulated and observed concentration-time profile of melatonin in saliva. The blue line represents the simulated profile, and the red line represents the observed data. a–c DSWPD patients. d–f Normal sleeping subjects.DiscussionThe population PK/PD model successfully characterized the PK and PD profile of drug, capturing the variability in drug exposure and response across the study population [17]. In this study, a comprehensive population PD model was developed to characterize the properties of endogenous melatonin using the population PK/PD approach, leveraging pooled data from 10 normal sleepers and 26 DSWPD patients. Unlike exogenous drugs, endogenous melatonin cannot be studied using traditional PK non-compartment analysis or standard compartment models, as these approaches rely on the administration of an exogenous drug and are not suitable for capturing the unique dynamics of endogenous substances. To address this, an indirect response PD model combined with a surge function was employed to accurately describe the production and elimination of melatonin. The surge function, Surge(t), effectively captures the time-dependent surge in melatonin production, reflecting the natural rise and fall of melatonin levels during the night. Charles et al. [22] were the first to describe the circadian rhythm characteristics of melatonin in the blood of healthy males using an indirect effect model that incorporated a surge function. It was reported that the concentration of melatonin detected in saliva is 31.6 ± 2.2% of that in blood, with a significant correlation observed between saliva and blood melatonin levels during both the production and elimination phases [31]. Therefore, the indirect effect model incorporating a surge function was equally applicable to the fitting of melatonin data from saliva. Compared with study by Charles et al. [22], our study had several advantages. We collected saliva samples from the subjects instead of blood samples, increasing patient compliance. We used an indirect effect model with a surge function not only to characterize the circadian profile of melatonin, but also to enable quantitative comparative analyses between the circadian profiles of normal sleepers and DSWPD patients, as well as advanced our understanding of DSWPD etiology and enabling prediction of melatonin concentrations at different time points based on sparsely sampled data.Previous research indicated that melatonin production began at 3 to 4 months of age and increased throughout childhood, typically peaking between 8 and 10 years of age, followed by a sudden decline during puberty [26]. Waldhauser et al. [32] and Arendt et al. [33] found no significant age dependency of melatonin concentration in adults aged 20 to 39 years. In the present study, although there was some difference in age between normal sleepers (mean age: 33.6 ± 11.22 years) and DSWPD patients (mean age: 21.73 ± 4.98 years), the effect of this difference on melatonin levels was limited. In the indirect effect model incorporating the surge function (Equation 3), N represented the exponent of the surge slope. The shape of the surge function curve could be modified by adjusting the value of N to better fit the observed data and improve the accuracy of the model. In our study, when N = 12, our model exhibited the optimal fit, with the CV% of the typical values of the model estimates all being less than 30%. Furthermore, the low values of −2LL, AIC, and BIC indicated a high degree of precision and reliability. In addition, we investigated the correlations between T0, Kout, AMP, Baseline, and WID using a covariance model. The results indicated that considering the correlations among AMP, Baseline, and WID led to a significant optimization of our model, with the −2LL value decreasing by 17 compared to the model that did not account for parameter IIV correlations. The inter-individual variability correlation between AMP, Baseline, and WID suggested the individual’s secretion and metabolism of melatonin during which these parameters did not exist independently but rather interacted with each other. For example, higher Baseline concentrations were associated with greater amplitudes, suggesting that an individual’s physiological state in circadian regulation may have influenced their melatonin secretion pattern. Conversely, alterations in WID may also have interacted with changes in amplitude and baseline concentration to affect overall melatonin levels.It had been shown that DLMO is delayed by approximately 2–6 h in DSWPD patients compared to normal sleepers [15, 16, 34‒36]. Our findings indicated that, compared to normal sleeper, individuals with DSWPD exhibited a significant delay in their melatonin circadian rhythm, with the DLMO delayed by approximately 7 h. It was currently postulated that a delay in the timing of the melatonin circadian rhythm was responsible for DSWPD. The International Classification of Sleep Disorders-Third Edition (ICSD-3) recommended the measurement of a patient’s DLMO for the diagnosis of circadian rhythm sleep-wake disorders [37]. In addition, it was reported that sleep latency in DSWPD patients might have been related to a reduced rate of melatonin secretion [15, 34]. However, no study was reported to compare the rates of melatonin production and elimination between normal sleepers and DSWPD patients. In contrast, we applied an indirect effect model including a surge function to more directly relate endogenous melatonin production and elimination to melatonin concentration. Our findings suggested that the rate of melatonin production was significantly lower in DSWPD patients compared to normal sleepers (2.33 ± 3.64 pg/h vs. 7.35 ± 10.91 pg/h). The lower rate of melatonin production in DSWPD patients might have been due to a variety of factors including dysregulation of the biological clock, the effects of light, genetic factors, and psychological stress. Therefore, it was possible that the cause of DSWPD, in addition to the delay in the timing of the melatonin circadian rhythm, might have been due to the reduced rate of melatonin production.The endogenous melatonin profile was a reliable marker of circadian rhythm phase [38]. DSWPD treatment was based on the principle of aligning the endogenous circadian rhythms with the timing of sleep. Exogenous melatonin affected circadian rhythms in a phase-dependent manner, making oral melatonin a common treatment for DSWPD. The administration of melatonin doses of approximately 1–5 mg had been demonstrated to be the most efficacious in the treatment of DSWPD [15]. Due to the chronopharmacological properties of melatonin, the timing of melatonin administration was critical in the treatment of DSWPD. The optimal timing of melatonin administration was contingent upon DLMO. Oral administration of melatonin between 1.5 and 6.5 h prior to DLMO had been demonstrated to advance circadian rhythm phases, with the magnitude of phase advancement exhibiting a linear correlation with the timing of melatonin administration [39]. It had been demonstrated that a discrepancy of merely a few hours in the time of initial administration could result in a threefold difference in the phase-advancing effect of melatonin on circadian rhythms. This variability represents the distinction between therapeutic success and failure [40]. Therefore, individualized treatment for DSWPD was essential. However, current studies on determining DLMO and melatonin secretion characteristics predominantly relied on collecting saliva or plasma samples from subjects every hour or 30 min. This method required repeated sampling (usually lasting at least 6 h) to assess the time at which an individual’s melatonin secretion exceeds a certain threshold [41]. The intensive sampling of subjects at night not only reduces patient compliance, but also requires stringent environmental and personnel requirements, and expensive testing, making clinical diagnosis as well as individualization of patients’ medication very difficult. Our approach incorporated sparse data into model, and the estimated DLMO of DSWPD patients calculated through the established final model showed less than a 10% error compared to original DLMO observed values (Table 4; Fig. 6).These results suggested that our final indirect effect model can effectively predict melatonin circadian rhythm characteristics and DLMO in patients with DSWPD using sparse sample data. This approach offers a valuable tool for personalizing medication strategies for DSWPD patients, with the potential to enhance treatment outcomes. The model was developed using intensive data from 36 subjects, comprising a total of 1,740 sample points, which allowed for robust population modeling and meaningful results. However, the relatively small number of subjects may limit the generalizability of our findings. Furthermore, the lack of individual covariate information, such as age, weight, ethnicity, and living conditions which were not available, might constrain the broader applicability of the results. Future studies should aim to include larger and more diverse populations, incorporate comprehensive demographic data to validate and extend the applicability of our model, and strictly control dietary intake, as melatonin is present in natural products such as honey and plant-derived substances [42], in addition to alcohol, coffee, and cigarettes.ConclusionIn this study, a population PK/PD approach was employed to quantitatively compare normal sleepers’ melatonin circadian rhythm characteristics with those of DSWPD patients. In addition to delayed DLMO, patients with DSWPD had significantly lower rate of melatonin production compared to normal sleepers. Based on the final model, the Bayesian method can be used to estimate melatonin concentration at different time points and DLMO in DSWPD patients using sparsely sampled data, which helps in individualizing medication for DSWPD patients.AcknowledgmentsWe thank the doctors, nurses, and sample collection personnel involved in this study. Special thanks to the subjects who consented to take part in this study.Statement of EthicsThis study protocol was reviewed and approved by Ethics Committee of The First Affiliated Hospital of Nanjing Medical University, Approval No. 2023-SR-431. The age of the subjects ranges from 18 to 60 years. All subjects have signed an informed consent form prior to their enrollment in the study.Conflict of Interest StatementThe authors have no conflicts of interest to declare.Funding SourcesThis study was not supported by any sponsor or funder.Author ContributionsZhanhui Lv: writing – original draft, data curation, and methodology. Zimo Zhang: writing – original draft, data curation, and data analysis. Lu Wang: methodology. Sufeng Zhou: project administration. Jianguo Sun: writing – review. Feng Shao: writing – review and editing, supervision, methodology, resources, and project administration. Xiuqin Wang: writing – review, resources, and project administration. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Neuroendocrinology Karger

Indirect Effect Model with Surge Function for Describing Melatonin Circadian Rhythm: Quantitative Comparison and Application between Normal Sleepers and Patients with Delayed Sleep-Wake Phase Disorder

Download PDF
14 pages

Loading next page...
 
/lp/karger/indirect-effect-model-with-surge-function-for-describing-melatonin-3RdnBgpDM4

References (41)

Publisher
Karger
Copyright
© 2025 The Author(s). Published by S. Karger AG, Basel
ISSN
0028-3835
eISSN
1423-0194
DOI
10.1159/000546125
Publisher site
See Article on Publisher Site

Abstract

IntroductionDelayed sleep-wake phase disorder (DSWPD) is a circadian rhythm disorder characterized by delayed sleep phases that cause difficulty falling asleep or waking up on time, leading to daytime dysfunction [1]. The prevalence of DSWPD in adults ranges from 1.51% to 8.90% [2]. In individuals with DSWPD, the endogenous melatonin rhythm is delayed and no longer aligns with desired sleep times [3]. Patients with this disorder are usually unable to fall asleep at the desired bedtime, which can result in sleep deprivation, daytime sleepiness, and fatigue, which can have a serious impact on the patient’s life [4, 5]. Consequently, it is of paramount importance to implement timely intervention and treatment for DSWPD.The principal therapeutic modalities for DSWPD encompass light therapy, cognitive-behavioral therapy, and melatonin treatment, all of which are chronobiological in nature [6‒8]. Melatonin is regarded as the optimal pharmacological intervention for the treatment of DSWPD in patients who have not responded to non-pharmacological treatments [3]. Exogenous melatonin can alter the timing of the biological clock, and depending on the time of administration, melatonin can advance or delay the phase of the biological clock [9]. The advancement of the phase occurs when melatonin is administered between 2 and 7 h prior to the dim light melatonin onset (DLMO). Conversely, the phase is delayed when melatonin is administered in the early morning [10]. For patients with DSWPD, the goal of treatment is to advance the circadian time so that patients can sleep at a more regular time. Due to the individual differences in the circadian rhythm of melatonin, treatment with melatonin may not only be ineffective, but may even be counterproductive or cause adverse effects such as drowsiness and disruption of the biological clock if it is not administered in accordance with the individual patient’s circadian rhythm, which is usually assessed using a DLMO [3, 11]. Therefore, selecting the appropriate time of administration based on the circadian rhythm characteristics of melatonin is crucial for individualized treatment of DSWPD.Melatonin is a neuroendocrine hormone primarily synthesized in the pineal gland, though it is also produced in peripheral organs such as the skin, heart, liver, uterus, placenta, kidneys, immune cells, and thymus. However, only the secretion from the pineal gland and the visual system is regulated by the environmental light/dark cycle, as other tissues lack this regulatory mechanism [12, 13]. Melatonin secretion is regulated by circadian rhythms, usually increasing at night and decreasing during the day, which helps regulate sleep and wake cycles [14]. By comparing the circadian rhythm characteristics of DSWPD patients and normal sleepers, the etiology of DSWPD patients can be more accurately assessed, thereby providing more references for the development of individualized treatment plans. Micic et al. [15] observed that the average melatonin concentration-time curve is steeper in the ascending phase in normal sleepers than in DSWPD patients, suggesting that the reduced rate of melatonin production may contribute to the pathogenesis of DSWPD patients. Shibui et al. [16] discovered that the melatonin circadian rhythm in DSWPD patients is significantly delayed compared to normal sleepers by analyzing their average melatonin concentration-time curves.Population pharmacokinetic/pharmacodynamic (PK/PD) modeling is a powerful tool for characterizing the variability in drug exposure and response across individuals, integrating both PK (drug concentration) and PD (drug effect) data to describe the dose-response relationship [17, 18]. This approach is particularly useful for understanding inter-individual variability and optimizing dosing regimens in diverse populations. In contrast, indirect effect models focus primarily on PD effects, often used to describe the time-dependent changes in endogenous substances or biomarkers without directly modeling drug concentrations [19].Indirect PD models are particularly valuable for characterizing endogenous substances, as they account for the production and elimination processes of these compounds. Nagaraja et al. [20] utilized an indirect effect model including a surge function to characterize the delay in the luteinizing hormone surge following single or multiple administrations of cetrorelix. Similarly, an indirect effect model with a surge function has been successfully employed to describe the circadian rhythms of ACTH, cortisol, and melatonin in healthy males [21, 22].In this study, we aimed to apply an indirect effect model combined with a surge function to characterize the properties of endogenous melatonin in normal sleepers and DSWPD patients. By leveraging this approach, we seek to elucidate the circadian rhythm of melatonin production and elimination, providing insights into its role in sleep regulation and potential therapeutic applications.MethodsData SourceThe modeling dataset comprised data from 36 subjects, including data from 10 normal sleepers and 26 DSWPD patients. Additional data from three normal sleepers and 3 DSWPD patients were used for model application (Fig. 1).Fig. 1.Flowchart of model building process.Modeling data for 10 normal sleepers (mean age 33.60 ± 11.22 years) were obtained from a clinical study (registration number: NCT06467851). Subjects were admitted to the study ward before 16:00, and saliva samples were collected between 16:00 and 11:00 the following day. Samples were collected at 1-h intervals between 16:00 and 7:00 the following day, and at 2-h intervals between 7:00 and 11:00. Saliva samples were collected in the study ward in dim light. Furthermore, the use of light-producing electronic devices, such as mobile phones, was prohibited. Products that could affect salivary melatonin concentrations, such as alcohol, coffee, and strong tea, are prohibited throughout the study period. Melatonin concentrations were determined using the salivary melatonin ELISA (IBL International GmbH, Hamburg, Germany, RE54041). In daytime samples (<8 pg/mL), the mean coefficients of variation within and between groups were 17.0% and 20.5%, respectively. In nighttime samples (>10 pg/mL), the mean coefficients of variation within and between groups were 13.9% and 18.4%, respectively.Modeling data for 26 DSWPD patients were digitally extracted from individual patient profiles reported by Micic et al. [15]. The original study by Micic et al. [15] met the following screening criteria: (1) patients with DSWPD were included; (2) patients were not treated with any medication that affects endogenous melatonin levels or with phototherapy interventions; (3) salivary melatonin was measured as the study endpoint; (4) the characteristics of a complete 24-h melatonin secretion were assessed; (5) subjects were sampled in dim light; (6) strictly control the intake of products that affect salivary melatonin level, such as alcohol, coffee, and cigarettes.Data collection for both normal sleepers and DSWPD patients was confined to autumn and winter, limiting the assessment of seasonal effects on the model. Complementing our findings, a comprehensive study by Titone et al. [23] examined salivary melatonin circadian rhythms across diverse racial groups, including Caucasians, African Americans, Asians, Pacific Islanders, and other ethnicities (e.g., individuals of Mexican, Native American, and Hispanic descent). Their results revealed no statistically significant racial differences in melatonin rhythms. Consequently, racial variation was not incorporated into the final model in this study. The data used for model application, employing the Bayesian estimation method, included data from three normal sleepers sourced from Burgess et al. [24] and data from 3 DSWPD patients sourced from Burgess et al. [25].Model DevelopmentThis study describes the circadian rhythm of endogenous melatonin using an indirect effect model that includes surge functions. Based on the final model, Monte Carlo simulations were conducted to quantitatively analyze the circadian rhythm characteristics of normal sleepers and patients with DSWPD. Lastly, Bayesian feedback was used to predict DLMO and melatonin concentration at different time points based on sparse data. The process of model building and application is shown in Figure 1.Structural ModelIn normal sleepers, melatonin synthesis begins at night, between 20:00 p.m. and 22:00 p.m., and reaches a maximum at midnight or between 2:00 a.m. and 3:00 a.m. After that, melatonin production gradually decreases and remains very low during the day [26]. Melatonin production and elimination can be described by an indirect effect model, as shown in formula (1) [20, 22, 27]:(1)dRdt=kin−kout×Rwhere R is the concentration of melatonin, Kin is the zero-order release rate constant of melatonin, and Kout is the first-order elimination rate constant of melatonin. Melatonin is a hormone secreted by the pineal gland. In mammals, the pineal gland is inhibited by light during the day and activated at night through a polysynaptic pathway in the SCN, which promotes melatonin secretion [28]. Therefore, the surge in melatonin at night may be attributed to an increase in its production rate. Charles et al. [22] have successfully applied a surge function to characterize the melatonin surge at night, so we chose an indirect effect model that included a surge function to represent the changes in melatonin over time (shown in Fig. 2), as shown in Equation (2-4) [20]:(2)dRdt=kin×Surget−Kout×R(3)Surget=1+AMPt−T0WIDN+1(4)Kin=Kout×BaselineIn Equations (2-4), Surge(t) denotes the surge function; Baseline refers the baseline level of melatonin; AMP indicates the amplitude of the peak, which is a dimensionless parameter related to the amplitude of the peak; t represents time; T0 is the time to peak; WID denotes the width of the peak; and N signifies the exponential of the slope of the peak, where N is an even number (e.g., 2, 4, 6, etc.) to ensure that the melatonin concentration remains positive both before and after the peak occurs. If t ≫ T0 or t ≪ T0, melatonin concentration is close to baseline, and if t ⟶ T0, there is a sharp increase in melatonin concentration.Fig. 2.An indirect effect model with surge function. R represents the PD response variable. Kin(t) denotes periodic production rates. Kout denotes the elimination rate. Dashed red arrow indicates stimulation. The surge function illustrated in this figure is adapted from the study by Olta Tarko [29].Random Effects ModelThe nonlinear mixed-effects model was established using Phoenix NLME (version 8.3.5, Certara, Co., Princeton, NJ, USA). The model parameters were estimated using first-order conditional estimation-extended least squares. The inter-individual variation (IIV) was expressed using an exponential model (Equation 5):(5)Pij=Pj×eηijPj represents the typical value of the jth parameter in the population, while Pij denotes the true value of the jth parameter for the ith subject. ηij denotes the IIV in the jth parameter for the ith subject. The IIV (η) conforms to a normal distribution with mean 0 and variance ω2. The residual variances were expressed using additive, proportional, and combined additive and proportional models (Equations 6, 7, and 8):(6)Cij=IPREDij+εij(7)Cij=IPREDij×1+εij(8)Cij=IPREDij×1+εij,1+εij,2Cij represents the observed concentration of the ith subject at the jth sampling point, whereas IPREDij denotes the model-predicted value for the subject. εij denotes the residual variance at the jth sampling point for the ith subject. The residual variance (ε) conforms to a normal distribution with mean 0 and variance σ2.Covariate and Covariance ModelsThe model examined the effect of disease state on melatonin circadian rhythms, so disease state was included in the model as a categorical variable. Assess the correlation between the IIV distributions of T0, Kout, AMP, Baseline, and WID through graphical analysis, and examine this relationship by incorporating the covariance between IIV parameters within the model.Model EvaluationThe degree of fit of the model to the observed data and the accuracy of the model were evaluated by diagnostic graphs (e.g., observations vs. population predictions [PRED], observations vs. individual predictions [IPRED], conditional weighted residuals [CWRES] vs. PRED, and CWRES vs. time). The predictive power of the final model was evaluated using the visual predictive check. A total of 1,000 simulations were performed for the final model, calculating the 5%, 50%, and 95% percentiles of both the observed data set and the simulated data set at various time points to compare their distribution characteristics. The bootstrap resampling method was employed 1,000 times to calculate the median values and 95% confidence intervals (2.5–97.5%) of the parameters, which were then compared to the original parameter estimates from the model to assess the robustness of the final model.SimulationMonte Carlo simulation (n = 1,000) was performed based on the final model, comparing the circadian rhythm characteristics of melatonin in normal sleepers and DSWPD patients, and DLMO was calculated from the simulation results. The absolute threshold method was used to calculate DLMO. DLMO is considered as the time when the salivary melatonin concentration reaches a threshold of 7 pg/mL [30].The external data used for Bayesian estimation consisted of data from three normal sleepers and 3 DSWPD patients, which were obtained from the references Burgess et al. [24] and Burgess et al. [25], respectively. The sparse data of three normal sleepers and 3 DSWPD patients were imported into Phoenix, respectively. Based on the established final model, the maximum a posteriori Bayesian method was employed to estimate melatonin concentrations at different time points and DLMO.ResultsSubjectsThe study included 1,740 concentrations from 36 subjects. A total of 180 concentrations were obtained from 10 normal sleepers, while 1,560 concentrations were obtained from 26 patients with DSWPD. The demographic characteristics and sleep quality scores of the 36 subjects are summarized in Table 1. The 10 normal sleepers consisted of 4 males and 6 females with a mean age of 33.6 ± 11.22 years, mean weight of 66.43 ± 7.07 kg, and a mean BMI of 24.57 ± 3.16 kg/m2. The DSWPD patients included 18 males and 8 females with an average age of 21.73 ± 4.98 years [15].Table 1.Summary of demographic characteristics and sleep quality scoresNormal sleeper (N = 10)DSWPD patients (N = 26)aMale, n (%)4 (40%)18 (69%)Female, n (%)6 (60%)8 (31%)Age, years33.6±11.2221.73±4.98Weight, kg66.43±7.07NABMI, kg/m224.57±3.16NAMEQ score49.4±7.2933.09±5.76PSQI score3.4±1.906.15±3.32Values are expressed as means ± SD.BMI, body mass index; NA, not applicable.aThe demographic characteristics and sleep quality scores of DSWPD patients were obtained from Micic et al. [18].Melatonin Population ModelAn indirect effect model incorporating a surge function was employed as the structural model, which effectively characterized the circadian rhythm distribution property of endogenous melatonin. Compared to additive models and additive-proportional error models, the proportional error models provided a better description of residual variance. An exponential error model was used to describe IIV. The inter-individual variants were included in the T0, Kout, AMP, Baseline, and WID. The final model estimation results are shown in Table 2. The model parameters were estimated with high precision, and their RES% values were low. Typical values of T0, Kout, AMP, Baseline, and WID estimated by the final model were 23:59, 1.23 h−1, 7.95, 3.21 pg/mL, and 4.12 respectively. The proportional-type residual variance of the model estimate was 0.21.Table 2.Parameter estimates of the final model and bootstrap evaluationParameterFinal modelBootstrapestimateRES, %95% CImedian95% CItv AMP7.9515.475.54∼10.368.005.79∼10.86tv T0, hh:mm23:594.1322:37∼01:2200:0222:44∼01:23tv WID4.125.783.65∼4.584.163.63∼4.52tv Kout, h−11.2321.820.70∼1.751.210.80∼1.93tv Baseline, pg/mL3.2123.271.75∼4.683.221.77∼6.45Inter-individual variabilityω2 T00.0119.950.006∼0.010.010.00∼0.01ω2 Kout1.0133.020.36∼1.660.910.32∼1.63ω2 AMP0.8526.670.41∼1.300.820.38∼1.28ω2 Baseline0.2153.38−0.01∼0.430.18−0.02∼0.42ω2 WID0.0672.33−0.03∼0.150.06−0.02∼0.16Residual variability (σ)Stev0a0.216.720.18∼0.240.210.18∼0.24AMP, amplitude; T0, time of peak melatonin concentration; WID, a peak width parameter; Kout, first-order elimination rate constant of melatonin; Baseline, baseline melatonin concentration.aStev0 represents residual variability (or intra-individual variability).The disease status was included as a categorical variable affecting T0 and Baseline in the model, resulting in a reduction of −2LL by approximately 146 (3,769 vs. 3,623). The results of the graphical evaluation indicate a correlation among the IIVs of AMP, Baseline, and WID. Incorporating their covariance structure enhanced model fit, reducing −2LL by 17 (3,623 vs. 3,606). Therefore, the disease status was included as a categorical variable in the final model, and the covariances among AMP, Baseline, and WID were examined.Model EvaluationThe final model was diagnosed by goodness-of-fit plots, including observations versus PRED, observations versus IPRED, CWRES versus PRED, and CWRES versus time (shown in Fig. 3). Although the trend lines of the observations versus PRED scatterplot deviate slightly from the reference line, the trend line of the observations versus IPRED scatterplots is highly coincident with the reference line, indicating that the observations are well fitted to the model predictions. The CWRES was distributed symmetrically on both sides of the reference line (y = 0) and did not exhibit a significant trend change with group predicted values and time. This indicates that the model was unbiased and well fitted. The results of the visual predictive check indicate that observations at the 5th, 50th, and 95th percentiles exhibit similar distributional characteristics to those of the modeled data (shown in Fig. 4). The results of the nonparametric bootstrap are shown in Table 2. The parameter values estimated by the final model are close to the median of the bootstrap results and are within the bootstrap 95% confidence interval, indicating that the final model is stable.Fig. 3.Goodness-of-fit plots for the final model. a Observations (DV) versus PRED. b DV versus PRED logarithm. c DV versus IPRED. d DV versus IPRED logarithm. e CWRES versus time. f CWRES versus PRED.Fig. 4.VPC plot of the final model. The blue lines represent the 5th, 50th, and 95th percentiles of observed melatonin concentration over time. The black lines represent the 5th, 50th, and 95th percentiles of the predicted value over time. The shaded region describes the 95% prediction interval of the 5th, 50th, and 95th percentiles for predicting melatonin concentration (n = 1,000 simulations). VPC, visual predictive check.SimulationAfter adding disease state as a covariate to the model, the −2LL of the model decreased (3,769 vs. 3,623), suggesting that disease state affects the circadian rhythm of melatonin. Therefore, the resampling of normal sleepers and DSWPD patients produced two virtual groups. Figure 5 illustrates that endogenous melatonin secretion in both normal sleepers and DSWPD patients exhibits a distinctive circadian rhythm feature. The area under the curve (AUC) in normal sleepers was significantly greater than the AUC in patients with DSWPD, which was 3.1 times greater than the AUC in patients with DSWPD. The results of the final model simulation, as presented in Table 3, indicated that patients with DSWPD had a delayed T0 of approximately 4.06 h, their Baseline was approximately 2.44 pg/mL lower than that of normal sleepers, and the AMP and WID were similar to that of normal sleepers. Additionally, the DLMO for DSWPD patients was delayed by approximately 7 h compared to normal sleepers. The rate of melatonin elimination in DSWPD patients was similar to that of normal sleepers, but their production rate constant was significantly lower. One potential explanation for the delayed melatonin circadian rhythm observed in patients with DSWPD is a reduced rate of melatonin production.Fig. 5.Average melatonin concentration-time curve after 1,000 simulations of the final model. Red line represents normal sleeping subjects, and blue line represents DSWPD patients.Table 3.The parameter values outputted after 1,000 simulations of the final modelParametersNormal sleepersDSWPD patientsAMP12.13±14.1012.20±14.07T0, hh:mm00:05±1.6904:08±2.11WID4.26±1.104.25±1.09Baseline, pg/mL3.58±1.741.14±0.56Kout, h−12.04±2.582.03±2.73Kin, pg/h7.35±10.912.33±3.64t1/2, ha0.93±1.200.94±1.30DLMO, hh:mmb17:4200:54AUC, pg/mL*hb422.30134.50Values for characteristics are presented as mean ± SD unless otherwise indicated.AMP, amplitude; T0, time to peak; WID, the peak width; Baseline, baseline levels of melatonin; Kout, first-order elimination rate constant of melatonin; Kin, zero-order release rate constant of melatonin; t1/2, half-life; DLMO, dim light melatonin onset; AUC, area under the melatonin concentration-time curve.aThe half-life is calculated based on the Kout, using the formula: t1/2 = 0.693/Kout.bDLMO and AUC were calculated from the average melatonin concentration-time curves generated using simulated data for normal sleepers and DSWPD patients.Sparse data from 3 normal sleepers and 3 DSWPD patients were imported into Phoenix, respectively. Based on the established final model, the maximum posteriori Bayesian method was used to estimate melatonin concentration at different time points and DLMO. The final results indicated that the DLMO estimated by the model for DSWPD patients and normal sleepers was similar to the DLMO calculated from the original data, with an error margin within ±10% (Table 4). Additionally, the circadian rhythm characteristics of melatonin for both DSWPD patients and normal sleepers derived from the model show a high degree of similarity to the literature data, with a significant overlap in the melatonin concentration-time profiles (shown in Fig. 6).Table 4.The DLMO of the final model estimate and the DLMO of the original data calculationClassSubjectsDLMOapredicted data, hh:mmObserved data, hh:mmPEbDSWPD patients101:2601:28−0.11%201:2801:250.33%322:0221:570.54%Normal sleepers423:4723:390.78%523:4523:380.72%622:1321:265.47%aDLMO denotes the dim light melatonin onset.bPE denotes the prediction error, which is given by the formula: PE = [(prediction data−observed data)/(observed data)] × 100%.Fig. 6.Comparison of simulated and observed concentration-time profile of melatonin in saliva. The blue line represents the simulated profile, and the red line represents the observed data. a–c DSWPD patients. d–f Normal sleeping subjects.DiscussionThe population PK/PD model successfully characterized the PK and PD profile of drug, capturing the variability in drug exposure and response across the study population [17]. In this study, a comprehensive population PD model was developed to characterize the properties of endogenous melatonin using the population PK/PD approach, leveraging pooled data from 10 normal sleepers and 26 DSWPD patients. Unlike exogenous drugs, endogenous melatonin cannot be studied using traditional PK non-compartment analysis or standard compartment models, as these approaches rely on the administration of an exogenous drug and are not suitable for capturing the unique dynamics of endogenous substances. To address this, an indirect response PD model combined with a surge function was employed to accurately describe the production and elimination of melatonin. The surge function, Surge(t), effectively captures the time-dependent surge in melatonin production, reflecting the natural rise and fall of melatonin levels during the night. Charles et al. [22] were the first to describe the circadian rhythm characteristics of melatonin in the blood of healthy males using an indirect effect model that incorporated a surge function. It was reported that the concentration of melatonin detected in saliva is 31.6 ± 2.2% of that in blood, with a significant correlation observed between saliva and blood melatonin levels during both the production and elimination phases [31]. Therefore, the indirect effect model incorporating a surge function was equally applicable to the fitting of melatonin data from saliva. Compared with study by Charles et al. [22], our study had several advantages. We collected saliva samples from the subjects instead of blood samples, increasing patient compliance. We used an indirect effect model with a surge function not only to characterize the circadian profile of melatonin, but also to enable quantitative comparative analyses between the circadian profiles of normal sleepers and DSWPD patients, as well as advanced our understanding of DSWPD etiology and enabling prediction of melatonin concentrations at different time points based on sparsely sampled data.Previous research indicated that melatonin production began at 3 to 4 months of age and increased throughout childhood, typically peaking between 8 and 10 years of age, followed by a sudden decline during puberty [26]. Waldhauser et al. [32] and Arendt et al. [33] found no significant age dependency of melatonin concentration in adults aged 20 to 39 years. In the present study, although there was some difference in age between normal sleepers (mean age: 33.6 ± 11.22 years) and DSWPD patients (mean age: 21.73 ± 4.98 years), the effect of this difference on melatonin levels was limited. In the indirect effect model incorporating the surge function (Equation 3), N represented the exponent of the surge slope. The shape of the surge function curve could be modified by adjusting the value of N to better fit the observed data and improve the accuracy of the model. In our study, when N = 12, our model exhibited the optimal fit, with the CV% of the typical values of the model estimates all being less than 30%. Furthermore, the low values of −2LL, AIC, and BIC indicated a high degree of precision and reliability. In addition, we investigated the correlations between T0, Kout, AMP, Baseline, and WID using a covariance model. The results indicated that considering the correlations among AMP, Baseline, and WID led to a significant optimization of our model, with the −2LL value decreasing by 17 compared to the model that did not account for parameter IIV correlations. The inter-individual variability correlation between AMP, Baseline, and WID suggested the individual’s secretion and metabolism of melatonin during which these parameters did not exist independently but rather interacted with each other. For example, higher Baseline concentrations were associated with greater amplitudes, suggesting that an individual’s physiological state in circadian regulation may have influenced their melatonin secretion pattern. Conversely, alterations in WID may also have interacted with changes in amplitude and baseline concentration to affect overall melatonin levels.It had been shown that DLMO is delayed by approximately 2–6 h in DSWPD patients compared to normal sleepers [15, 16, 34‒36]. Our findings indicated that, compared to normal sleeper, individuals with DSWPD exhibited a significant delay in their melatonin circadian rhythm, with the DLMO delayed by approximately 7 h. It was currently postulated that a delay in the timing of the melatonin circadian rhythm was responsible for DSWPD. The International Classification of Sleep Disorders-Third Edition (ICSD-3) recommended the measurement of a patient’s DLMO for the diagnosis of circadian rhythm sleep-wake disorders [37]. In addition, it was reported that sleep latency in DSWPD patients might have been related to a reduced rate of melatonin secretion [15, 34]. However, no study was reported to compare the rates of melatonin production and elimination between normal sleepers and DSWPD patients. In contrast, we applied an indirect effect model including a surge function to more directly relate endogenous melatonin production and elimination to melatonin concentration. Our findings suggested that the rate of melatonin production was significantly lower in DSWPD patients compared to normal sleepers (2.33 ± 3.64 pg/h vs. 7.35 ± 10.91 pg/h). The lower rate of melatonin production in DSWPD patients might have been due to a variety of factors including dysregulation of the biological clock, the effects of light, genetic factors, and psychological stress. Therefore, it was possible that the cause of DSWPD, in addition to the delay in the timing of the melatonin circadian rhythm, might have been due to the reduced rate of melatonin production.The endogenous melatonin profile was a reliable marker of circadian rhythm phase [38]. DSWPD treatment was based on the principle of aligning the endogenous circadian rhythms with the timing of sleep. Exogenous melatonin affected circadian rhythms in a phase-dependent manner, making oral melatonin a common treatment for DSWPD. The administration of melatonin doses of approximately 1–5 mg had been demonstrated to be the most efficacious in the treatment of DSWPD [15]. Due to the chronopharmacological properties of melatonin, the timing of melatonin administration was critical in the treatment of DSWPD. The optimal timing of melatonin administration was contingent upon DLMO. Oral administration of melatonin between 1.5 and 6.5 h prior to DLMO had been demonstrated to advance circadian rhythm phases, with the magnitude of phase advancement exhibiting a linear correlation with the timing of melatonin administration [39]. It had been demonstrated that a discrepancy of merely a few hours in the time of initial administration could result in a threefold difference in the phase-advancing effect of melatonin on circadian rhythms. This variability represents the distinction between therapeutic success and failure [40]. Therefore, individualized treatment for DSWPD was essential. However, current studies on determining DLMO and melatonin secretion characteristics predominantly relied on collecting saliva or plasma samples from subjects every hour or 30 min. This method required repeated sampling (usually lasting at least 6 h) to assess the time at which an individual’s melatonin secretion exceeds a certain threshold [41]. The intensive sampling of subjects at night not only reduces patient compliance, but also requires stringent environmental and personnel requirements, and expensive testing, making clinical diagnosis as well as individualization of patients’ medication very difficult. Our approach incorporated sparse data into model, and the estimated DLMO of DSWPD patients calculated through the established final model showed less than a 10% error compared to original DLMO observed values (Table 4; Fig. 6).These results suggested that our final indirect effect model can effectively predict melatonin circadian rhythm characteristics and DLMO in patients with DSWPD using sparse sample data. This approach offers a valuable tool for personalizing medication strategies for DSWPD patients, with the potential to enhance treatment outcomes. The model was developed using intensive data from 36 subjects, comprising a total of 1,740 sample points, which allowed for robust population modeling and meaningful results. However, the relatively small number of subjects may limit the generalizability of our findings. Furthermore, the lack of individual covariate information, such as age, weight, ethnicity, and living conditions which were not available, might constrain the broader applicability of the results. Future studies should aim to include larger and more diverse populations, incorporate comprehensive demographic data to validate and extend the applicability of our model, and strictly control dietary intake, as melatonin is present in natural products such as honey and plant-derived substances [42], in addition to alcohol, coffee, and cigarettes.ConclusionIn this study, a population PK/PD approach was employed to quantitatively compare normal sleepers’ melatonin circadian rhythm characteristics with those of DSWPD patients. In addition to delayed DLMO, patients with DSWPD had significantly lower rate of melatonin production compared to normal sleepers. Based on the final model, the Bayesian method can be used to estimate melatonin concentration at different time points and DLMO in DSWPD patients using sparsely sampled data, which helps in individualizing medication for DSWPD patients.AcknowledgmentsWe thank the doctors, nurses, and sample collection personnel involved in this study. Special thanks to the subjects who consented to take part in this study.Statement of EthicsThis study protocol was reviewed and approved by Ethics Committee of The First Affiliated Hospital of Nanjing Medical University, Approval No. 2023-SR-431. The age of the subjects ranges from 18 to 60 years. All subjects have signed an informed consent form prior to their enrollment in the study.Conflict of Interest StatementThe authors have no conflicts of interest to declare.Funding SourcesThis study was not supported by any sponsor or funder.Author ContributionsZhanhui Lv: writing – original draft, data curation, and methodology. Zimo Zhang: writing – original draft, data curation, and data analysis. Lu Wang: methodology. Sufeng Zhou: project administration. Jianguo Sun: writing – review. Feng Shao: writing – review and editing, supervision, methodology, resources, and project administration. Xiuqin Wang: writing – review, resources, and project administration.

Journal

Neuroendocrinology – Karger

Published: Aug 1, 2025

Keywords: Melatonin; Circadian rhythms; Delayed sleep-wake phase disorder; Quantitative comparison; Nonlinear mixed-effect modeling

There are no references for this article.