1Division of Sleep Medicine and 2Division of Endocrinology, Diabetes, and Hypertension, Department of Medicine, Brigham and Women's Hospital, Harvard Medical School, Boston 02115; and 3Neuroscience Statistics Research Laboratory, Department of Anesthesia and Critical Care, Massachusetts General Hospital, Division of Health Sciences and Technology, Harvard Medical School/Massachusetts Institute of Technology, Boston, Massachusetts 02114
Submitted 23 December 2002 ; accepted in final form 8 July 2003
![]() |
ABSTRACT |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
compartment model; constant routine
Under the control of the hypothalamic-pituitary axis (Fig. 1), GH is continuously synthesized and packaged into secretory vesicles in the anterior pituitary (2). The neuropeptides GH-releasing hormone (GHRH) and somatotropin releasing-inhibiting factor (SRIF) are the principal hypothalamic regulators of GH secretion. GHRH, synthesized primarily in the arcuate nucleus of hypothalamus, and SRIF, produced in the arcuate and periventricular nuclei of the hypothalamus, are released into the hypothalamic-hypophysial portal system from neuronal processes terminating at the median eminence. Ghrelin is a recently identified GH secretagogue produced in the gut that interacts with G protein-coupled GH secretagogue receptors in the arcuate nucleus and ventromedial nucleus of the hypothalamus and in the pituitary to stimulate GH secretion (17); its potency is similar to that of GHRH (23). GHRH and SRIF are transported through the portal system to the anterior pituitary, where they activate their respective G protein-coupled receptors located on anterior pituitary somatotrophs. GHRH stimulates an increase in cytosolic calcium and cAMP, an increase in GH gene transcription and GH biosynthesis, and a release of GH into the circulation. Individuals with an inactivating mutation of the GHRH receptor show a markedly reduced 24-h GH serum level with a continued episodic pattern (30). SRIF, acting as a noncompetitive antagonist of GHRH, reduces cAMP and cytosolic calcium and inhibits GH secretion (2, 9, 31). SRIF appears to affect the timing and amount of GH released in secretory episodes and basal GH secretion but not GH synthesis (10).
|
In response to GHRH, there is activation of the cAMP/protein kinase A pathway leading to an influx of calcium ions into somatotropic cells, fusion of the GH secretory vesicles with the cell membrane, and release of GH into the serum (6). GH dissolves readily in the serum and is transported to its target sites. In the liver and other target sites, GH stimulates production of insulin-like growth factors IGF-I and IGF-II, which exert growth-promoting effects in peripheral tissues. These growth factors exert negative feedback on GH serum levels by suppressing GHRH and stimulating SRIF release from the hypothalamus, as well as by directly inhibiting GH secretion from the anterior pituitary (10, 13). GH may also exert a direct negative feedback on its secretion by regulating secretion of GHRH and SRIF (1, 32). GH is cleared by the kidney and by receptor-mediated uptake by the liver (11, 18). The clearance of GH has been successfully modeled as first-order kinetics (15, 16, 26, 27).
A complete mathematical description of diurnal GH variation would include all of the neural and humoral feedforward and feedback linkages between the hypothalamus, the anterior pituitary, and the peripheral GH regulators, as well the effects of exogenous factors such as sleep state, stress, and meals. At present, simultaneous measurement of the variables necessary to define such a model is not possible in humans, and only possible to a limited degree in other species. An alternative approach that we have used successfully in the study of melatonin and cortisol (4, 5) is to construct a minimal model to describe the observed diurnal patterns based on known physiology. If this model is sufficiently parsimonious, then the parameters can be estimated from the experimental data of individual subjects and the information pooled to define population differences in normal and diseased conditions. Specifying a minimal model of the diurnal variation in GH serum levels is the approach we take here.
The essential features of the GH diurnal pattern to consider in formulating a model are the kinetic parameters of infusion into and clearance from the serum, the feedback of GH on its secretion, and the episodic structure of the secretory events. Although there is an observed circadian rhythm of GH secretion independent of sleep, the magnitude of the circadian effect is relatively small compared with the circadian effect on other hormones, such melatonin and cortisol (7). In this report, we present a two-dimensional linear differential equation model based on these physiological properties of GH and use it to analyze diurnal variation in GH serum levels in 12 normal, healthy women during sleep and during constant-routine conditions.
![]() |
MATERIALS AND METHODS |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() | (1) |
![]() | (2) |
To complete the model description, we assume that, during a time interval (0,T], we collect n blood samples, which are assayed for GH. We let yt1,...,ytn denote the concentration of GH at times 0 t1 < t2 <... < tn1 < tn
T and assume the yti values satisfy the equation
![]() | (3) |
Protocol and subjects. Healthy, adult, premenopausal women [age 2150 yr and body mass index (BMI) 19.030.0 kg/m2] were recruited using advertisements in local newspapers. All subjects underwent a detailed history and physical examination, including measurement of blood and urine chemistries and thyroid function tests. Subjects underwent a Structured Clinical Interview from the Diagnostic and Statistical Manual of Mental Disorders, 4th edition (DSM-IV) (33), and individuals with past or current medical or psychological problems were excluded from the study procedure. Any subject who had taken glucocorticoids within the past year or estrogen/progesterone in the past 4 mo was excluded. No subject was on medication. Any subject who had traveled across more than one time zone in the past 3 mo was excluded. All studies were reviewed and approved by the Human Research Committee of the Partners HealthCare System. Informed written consent was obtained from each subject. All studies took place at the General Clinical Research Center (GCRC) of the Brigham and Women's Hospital. Data from 12 healthy, adult women are included in this report. The mean (±SD) age of the subjects was 28.1 ± 7.5 yr (range 2345 yr), and the mean BMI was 24.0 ± 3.3 kg/m2 (range 19.529.6 kg/m2).
For 3 wk before a 5-day inpatient protocol, all subjects wore a wrist actigraph, kept a sleep/wake diary, and were instructed to adhere to an 8-h sleep/wake schedule based on their habitual schedule. Caffeine, vitamins, and herbal supplements were discontinued 2 wk before hospitalization. All events in the 5-day inpatient study were scheduled using the subjects' habitual sleep/wake time, which was calculated from the 3-wk sleep diary. Subjects were scheduled for admission to the GCRC during the follicular phase of their menstrual cycle. Within the GCRC, for the first three consecutive nights, subjects were scheduled to sleep for 8 h in the dark according to their habitual schedules. Electrocardiograms and polysomnographic recordings of sleep were obtained. Blood was sampled every 10 min through an indwelling intravenous catheter inserted 2 h before blood sampling. For GH, sampling began 7 h before the night 3 sleep episode and continued for the next 24 h through the first 16 h of the constant routine (CR, described below).
Beginning with the sleep episode on night 2, subjects were in an environment free of time cues (no clock, watch, radio, TV, or telephone) and in darkness (0.3 lux) during sleep episodes or dim light (<15 lux) when awake. Beginning 7 h before the night 3 sleep episode and continuing until the end of the protocol, subjects remained in bed, either supine (during the night 3 sleep episode) or semirecumbent (30° head-up tilt) during scheduled wake episodes. Beginning after awakening from the night 3 sleep episode at their habitual wake time, and continuing for 40 h, subjects were maintained in CR conditions of constant semirecumbent position, low light (<15 lux), continuous wakefulness, and frequent small meals.
Subjects consumed a controlled nutrient, isocaloric diet (125 meq sodium, 100 meq potassium, 1,000 mg calcium, and 2,500 ml of fluid) beginning 2 days before admission and continuing throughout the protocol. During inpatient days 1, 2, and 3, subjects were given three meals and two snacks. During the CR on days 4 and 5, food was evenly distributed.
GH assay. GH was assayed in duplicate using the Nichols Advantage Human Growth Hormone chemiluminescence immunoassay from Nichols Institute Diagnostics (San Juan Capistrano, CA). The lower limit of detection is 0.01 ng/ml with an intra-assay coefficient of variance (CV) of 4.28.0% and an interassay CV of 4.112.1%.
Model fitting. Let y = (yt1,...,ytn) be the vector of n GH observations collected on a subject. The number n may be different for each subject, because there may be isolated missing samples from some subjects. To keep our notation simple, we have not included these differences in our notation. Our objective was to estimate = (µ, A, ton,
), where A = (A1,...,AJ), ton = (ton,1,...,ton,J), and
= (
I,
C,
f), and J is the number of secretory events from each subject's GH data series. The sleep episode data analyses were started at the subject's habitual bedtime; the circadian data analyses were started at the subject's habitual wake time, when the CR protocol began. It follows from Eq. 3 that the joint probability density of y is the multivariate Gaussian density defined by
![]() | (4) |
![]() | (5) |
To decide on the best model fit, including the best estimates of the kinetic parameters and the best estimates of the number, amplitude, and onset times of the secretory events, we used the following model selection strategy. First, we made initial guesses at the secretory event locations by inspecting the plot of the GH data series. Second, for a given set of initial guesses of these locations, we fit the model by maximum likelihood to determine the best estimates of the model parameters, including the secretory event locations and their amplitudes. Third, the AIC for this model, the model parameter standard errors, and the change between the initial guess at each secretory onset time and its final estimate were computed. The AIC assesses goodness-of-fit by measuring the trade-off between minimizing the 2log likelihood function (closeness to the data) and the number of parameters required to achieve that closeness (cost). By repeating steps 13 with different locations and numbers of secretory events and choosing the model for which the AIC is a minimum, we have an objective guide to determining the location and number of the secretory events, along with the other model parameters. We also refit the model if 1) any estimated secretory event onset time differed by >30 min from the initial estimate; 2) any component of A was 0.0; or 3) the standard error of any parameter was greater than the value of the parameter itself. The refit procedure involved adding, removing, or reestimating one or more secretory event onset times. For each individual data set, the model was fit with and without feedback terms.
![]() |
RESULTS |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
|
Model parameter analysis. The parameter estimates for I ranged from 6.0 to 27.0 min, with a median of 13.8 min for data during the sleep episode (Table 1), whereas they ranged from 5.4 to 55.8 min, with a median of 12.6 min, for data during the CR (Table 2). With the exception of subjects 1, 2, 10, and 12, the estimates of
I were not significantly different between the sleep episode and the CR. Across all subjects, the estimates of
I were not significantly different between the sleep episode and the CR (paired t-test P = 0.53). The parameter estimates of
C ranged from 6.6 to 63.0 min, with a median of 13.8 min for the sleep episode, whereas they ranged from 10.2 to 19.2 min, with a median of 11.7 min, for the CR. The range of
C estimates was much wider for data from the sleep episode than for those from the CR. For the estimate of
C, six of the subjects (subjects 2, 3, 4, 7, 8, and 10) were significantly different for the sleep episode compared with the CR. Across all subjects, the estimates of
C were not significantly different between the sleep episode and the CR (paired t-test P = 0.31). For all of the subjects, the estimated level of basal secretion was small. The estimates of basal secretion, µ, ranged from <0.01 to 0.23 ng/l for the sleep episode, with a median of 0.11 ng/l, compared with a range of <0.01 to 0.15 ng/l and a median of 0.08 ng/l for the CR. Across all subjects, the level of basal secretion was not different between the two conditions (paired t-test P = 0.18).
|
|
The number of secretory events differed significantly between the 8-h sleep episode and the 16-h CR (paired t-test P = 0.0004). However, the number of secretory events per hour did not differ between the sleep episode and the CR (paired t-test P = 0.11). The number of secretory events was either two or three during the sleep episode, whereas it ranged from one to nine, with a median of five, for the CR. For each subject, we computed for both the sleep and CR episodes the mean GH secretory influx, which for a given subject is the average of the At in each of the two conditions. The mean GH secretion per secretory event between the two conditions was significantly different only for subject 8, for whom it was 7.6 ng · l1 · h1 during the sleep episode compared with 34.2 ng · l1 · h1 during the CR.
We used the AIC goodness-of-fit criterion to evaluate the utility of including or not including the feedback term, f, in the analysis. The AIC values were lower (implying a more parsimonious fit) when the the model was used without the feedback term for 10 of the 12 subjects during the sleep episode compared with 4 of the 12 subjects during the CR. For 20 of the 24 model fits, the AIC values for the models with and without feedback were within 10% of each other. The model fits in Fig. 2 are those estimated without the feedback term. In all cases, as Fig. 3 illustrates, the model fits with and without feedback were graphically indistinguishable. For several of the data series, primarily those from the sleep episodes, the model parameter estimates from the fits with the feedback terms had larger standard errors.
|
![]() |
DISCUSSION |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
As we stated in the introductory comments, a complete mathematical description of diurnal variation in serum GH levels, like that of any pituitary hormone, should include all neural and humoral feedforward and feedback linkages between the hypothalamus, anterior pituitary, and other tissue sites, as well as the effects of exogenous factors. A mathematical model for GH in rats in a simulation study has been reported (40); however, that model was not fit to experimental data. A complete model of the hypothalamic-pituitary Leydig cell system has been reported by Keenan and Veldhuis (22). These authors used likelihood methods to fit the model to experimental data. In general, such complex models are difficult to specify in humans, because most of the components needed to define them cannot be known. Despite the complex physiology that underlies the GH diurnal patterns, the fact that our minimal linear differential equation model fits all the data series well suggests that the essential features of this system are a two-compartment model with first-order infusion and clearance, along with discrete secretory events of differing amplitudes. This minimal model, based on known physiology, can be used to help constrain the more complete model.
Like many other endocrine hormones, GH has been analyzed using deconvolution methods, pulse-finding algorithms, and approximate entropy (ApEn) procedures (14, 19, 20, 21, 36, 37, 39). Pulse-finding algorithms use local statistical criteria, such as Bonferroni bounds, t-tests, or the immunoassay CV, to identify the times of secretory events in hormone data series. Deconvolution models use deterministic linear time-invariant convolution integrals to represent the relation between inputs, secretory event times and amplitudes, and observed output hormone levels. In deconvolution analyses, the convolution integral is used to estimate the secretory event times and secretory amplitudes from the time series of serum hormone levels. In this sense, these algorithms perform simultaneous deconvolution and pulse finding. ApEn is a technique used to assess regularity and complexity of a biological series. Although the ApEn statistic does not have any direct physiological interpretation, it can characterize differences in physiological and pathophysiological states if these differences are reflected in differences in the degree of regularity of the hormone time series between the two states (39).
Our approach includes five extensions of the models and algorithms currently used to perform simultaneous deconvolution and pulse finding. First, the algorithms DECONV (37), PULSE 1 (20, 21), and PULSE 2 (20, 21) apply to neuroendocrine hormone data series in general. In our framework, we specify a model for each hormone based on its known physiology instead of assuming that the same deconvolution model and algorithm apply to every system. Using the specifics of the GH physiology in the model formulation makes apparent the modeling assumptions and the extent to which the estimated model components can be interpreted. This is why our GH model differs from our previous models for melatonin and cortisol (4, 5).
Second, current deconvolution algorithms fit one-dimensional models in which either Gaussian-shaped secretory event functions (DECONV) or impulse response functions (PULSE 1 and PULSE 2) are convolved with a one-compartment or a two-compartment clearance function. Our algorithm performs simultaneous deconvolution and pulse finding with a two-dimensional equation system. This is because the solution to the differential equation system in Eqs. 1 and 2 is a bivariate convolution integral (4, 5). The GH physiology we reviewed in the introductory paragraphs suggests that a minimal model for this system is two dimensional.
Third, in the DECONV algorithm, each pulse is represented as a Gaussian function, characterized by an amplitude and a width at half-maximum. The width at half-maximum is identical for all of the pulses, and each pulse has an infinite extent in both the future and the past. The waveform-independent deconvolution algorithms PULSE 1 and PULSE 2 use impulse response functions as an alternative to avoid specifying explicit pulse model. However, to do so, the original form of this algorithm, PULSE 1, assumed a secretory event at each data measurement point. Our approach suggests a straightforward compromise between use of either a Gaussian pulse model or a large number of waveform-independent pulses. This is because in our model (Eqs. 1 and 2), the serum level at any time, and hence a secretory event, is always defined causally by the dynamic interaction among four physiological quantities. These are the times of secretory event initiation, the amount of GH that is expelled from the releasable pituitary pool with each secretory event, the time constants for infusion, clearance, and feedback, and the basal level of secretion. Thus our model is parameterized like the Gaussian pulses yet with the flexibility of the waveform-independent functions.
Fourth, like PULSE 1 and PULSE 2, our algorithm directly estimates secretory event onset times. In contrast, the DECONV algorithm estimates the time at which a pulse is maximal. The onset time is more physiologically interpretable, because it defines the point at which a discrete biological signal is initiated in the GH axis. The effect of this signal is transmitted only into the future. Similarly, whereas we define the rate constant for GH infusion from the releasable pool, the comparable parameter in the DECONV algorithm is the width at half-maximum of the Gaussian pulse. This parameter, like the time of the pulse maximum, is more challenging to interpret physiologically. The PULSE 1 and PULSE 2 algorithms have no parameter comparable to our infusion parameter.
Finally, DECONV, PULSE 1, and PULSE 2 use nonlinear least squares to estimate the pulse locations, their amplitudes, clearance, basal secretion, and the width at half-maximum. For the DECONV algorithm, the nonlinear least squares is maximum likelihood. The PULSE 1 and PULSE 2 algorithms identify the pulse onset times, their amplitudes, and basal secretion with the assumption that the clearance parameter is known. PULSE 2 is identical to PULSE 1 except that the former uses an approximate F-test in a forward selection procedure to decide on the number of pulses. PULSE 2 and DECONV are used together to estimate all of the model parameters. In contrast, our approach uses one model family and maximum likelihood estimation, with AIC as a model selection criterion, to estimate the number of secretory events, their locations, amplitudes, and values of the kinetic parameters, and in particular, the presence or absence of feedback.
Our estimates of the clearance rate constants and the average secretory influx are consistent with previous reports (8, 15, 16, 26, 27, 34, 35). Because the GH secretory pattern is highly state dependent, we expected to find different parameter estimates for the sleep episode and for the CR. Although the infusion time constant was the same for the two conditions, the clearance time constants were different for 6 of the 12 subjects. In addition to reflecting real differences between the sleep episode and the CR, these differences in parameter estimates may also be due to the fact that there are more data from which to carry out the estimation from the 16-h CR than from the 8 h of data from the sleep episode. Another potential confounding effect is that the data presented in this report were collected continuously during 8 h of sleep followed by the first 16 h of a 40-h protocol of wakefulness. Therefore, we cannot distinguish between effects of sleep and circadian phase, since we did not have GH data collected at the same clock hour in each individual, both when she was awake and when she was asleep.
The data from our sample of female subjects are different from those reported for young men, in whom increased GH secretory event activity is observed with sleep onset and during sleep. However, data in those reports of men were not analyzed with this model and were done under different experimental conditions. Further experiments are required to dissect circadian and sleep effects on GH in men and women when this model is used.
On the basis of the current model and data set, we are unable to reliably estimate the model's feedback component for all of the data series. Although the presence of feedback upon GH release is well established, our analyses with and without feedback are indistinguishable. Estimation of the feedback term did not affect the times of secretory event onset or the GH influx rates per secretory event. There are several possible reasons we were unable to identify feedback. One is that feedback in this system is modulated differently from the way it is formulated in our current model, which assumes it to be related to the current serum level of GH. Second, data from other experimental conditions may be required to estimate the feedback term for the current model formulation. Third, for a given set of experimental conditions, more data may be needed to accurately and stably estimate this model component. The longer data recording time must be balanced against the requirement that the experimental conditions remain constant. Finally, a different estimation approach, such as Bayesian Monte Carlo methods (29), may allow us to determine the feedback reliably with the current model.
Several improvements to the current model will be considered in our future work. First, we did not use information on the immunoassay error as part of our fitting procedure as we did in our analysis of melatonin (4). Second, we used a simple feedback model in which the serum GH level decreased the rate of secretion from the releasable pool. Other formulations of feedback we will consider are the rate of change of serum GH modulating secretory release and the level of GH or its rate of change modulating GH production. Third, we used our plots of each hormone data series to identify the putative locations and hence the putative number of secretory events. The AICs were compared for each fixed number of secretory events, and the models with the smallest AIC values guided us in deciding on the most likely number of secretory events and their amplitudes and locations. As an alternative approach, we are developing variable-dimension Markov Chain Monte Carlo (MCMC) methods (29) to conduct the parameter and model order estimation simultaneously. This requires a stochastic formulation of the model in Eqs. 1 and 2. We have used MCMC methods previously to estimate a stochastic model of diurnal variation in cortisol plasma levels for a fixed number of secretory events (24). The variable dimension MCMC will allow us to estimate the number of secretory events at the same time that we estimate the other model components.
Our paradigm suggests a way to characterize GH physiology quantitatively under normal conditions, in pathological states, and after pharmacological or genetic manipulations in terms of the model component estimates. For example, low GH levels have been suggested as contributing to some of the adverse physiological changes associated with aging. Administration of drugs such as pyridostigmine, which decreases somatostatin tone, most certainly has state-dependent effects on GH physiology. A quantitative study of GH in animals that have been genetically altered is a promising approach to studying GH regulation. Model-based analyses of these questions may help us to better understand GH and the many complexities of its physiological roles.
In summary, we have developed a parsimonious physiologically based model in which we can estimate simultaneously the kinetic parameters, times of secretory events and their amplitudes, and basal secretion using maximum likelihood for model fitting and AIC to aid in model selection. This analysis paradigm should facilitate systematic study of normal and abnormal physiology of the GH axis.
![]() |
DISCLOSURES |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Note: The Windows-based program to fit our GH model to experimental data is available upon request from Dr. Klerman.
![]() |
ACKNOWLEDGMENTS |
---|
![]() |
FOOTNOTES |
---|
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
* E. B. Klerman and G. K. Adler are joint first authors of this work.
![]() |
REFERENCES |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|