## Abstract

The diurnal pattern of growth hormone (GH) serum levels depends on the frequency and amplitude of GH secretory events, the kinetics of GH infusion into and clearance from the circulation, and the feedback of GH on its secretion. We present a two-dimensional linear differential equation model based on these physiological principles to describe GH diurnal patterns. The model characterizes the onset times of the secretory events, the secretory event amplitudes, as well as the infusion, clearance, and feedback half-lives of GH. We illustrate the model by using maximum likelihood methods to fit it to GH measurements collected in 12 normal, healthy women during 8 h of scheduled sleep and a 16-h circadian constant-routine protocol. We assess the importance of the model components by using parameter standard error estimates and Akaike's Information Criterion. During sleep, both the median infusion and clearance half-life estimates were 13.8 min, and the median number of secretory events was 2. During the constant routine, the median infusion half-life estimate was 12.6 min, the median clearance half-life estimate was 11.7 min, and the median number of secretory events was 5. The infusion and clearance half-life estimates and the number of secretory events are consistent with current published reports. Our model gave an excellent fit to each GH data series. Our analysis paradigm suggests an approach to decomposing GH diurnal patterns that can be used to characterize the physiological properties of this hormone under normal and pathological conditions.

- compartment model
- constant routine

growth hormone (GH) is a peptide hormone critical for regulation of growth, development, and metabolism that is synthesized and released from somatotropic cells of the anterior pituitary (2, 10, 25). The 24-h profile of serum GH is characterized by multiple GH peaks due to the episodic release of GH from these cells. The timing and magnitude of the secretory episodes are affected by the time of day, arousal and sleep states, level of stress, nutritional state, body composition, stage of development, age, and gender of the individual (38). Reduced effective GH levels result in lower linear growth velocity in children and contribute to bone loss, decreased muscle mass and strength, increased abdominal fat, altered lipid profiles, and an impaired quality of life in adults (34). As a consequence, there is a compelling need to have models of the diurnal pattern of GH serum levels that can be used to study the physiology of GH secretion under different physiological and pathophysiological situations.

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

*Model formulation.* As described in our introductory comments, the amount of hormone in the releasable pool over time depends on the amount of GH synthesized and packaged into secretory vesicles and on the amount released into the circulation. These secretory vesicles constitute the releasable pool of GH, the GH available for immediate release upon stimulation with GH secretagogues. In addition, there may be a pool of GH that is secreted basally. We assume that release of GH into the circulation and clearance of GH from the circulation obey first-order kinetics. We also assume that the negative feedback of the circulating GH level on GH release is proportional to the circulating GH concentration. Let H_{1}(*t*) be the level of GH in the releasable pool, and let H_{2}(*t*) be the concentration of GH in the circulation due to discrete secretory events. The relation between GH in the releasable pool of the pituitary and the circulation can be described by the following first-order differential equations (1) (2) where β* _{I}* is the infusion rate constant; β

*is the clearance rate constant; β*

_{C}*is the feedback rate constant of GH; A*

_{f}*is the influx rate of GH into the releasable pool at time*

_{t}*t*(ng · l

^{–1}· h

^{–1}); and d

*N(t)*is an indicator function that is 1.0 if there is a secretory event initiated at time

*t*, and 0.0 if otherwise. The quantity A

*represents the hypothalamic and nonhypothalamic secretory inputs to the somatotrophs, including GHRH, SRIF, ghrelin, IGF-I, and IGF-II, that result in GH release during a secretory event. Although there is experimental evidence to suggest that clearance of GH from the circulation is first order (15, 16, 26, 27), our representation of GH infusion and feedback as also first order are assumptions made to construct our minimal model. The validity of these and other assumptions leading to the model components will be assessed in our goodness-of-fit and model selection analyses detailed below.*

_{t}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 y_{t}_{1},...,y* _{tn}* denote the concentration of GH at times 0 ≤

*t*

_{1}<

*t*

_{2}<... <

*t*

_{n}_{–1}<

*t*≤

_{n}*T*and assume the y

*values satisfy the equation (3) where μ is the serum basal level of GH, H*

_{ti}_{2}(

*t*) is the GH serum level at time

_{i}*t*due to discrete secretory events, and the ϵ

_{i}*are independent, zero mean, Gaussian random variables with variance σ*

_{ti}^{2}.

*Protocol and subjects.* Healthy, adult, premenopausal women [age 21–50 yr and body mass index (BMI) 19.0–30.0 kg/m^{2}] 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 23–45 yr), and the mean BMI was 24.0 ± 3.3 kg/m^{2} (range 19.5–29.6 kg/m^{2}).

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.2–8.0% and an interassay CV of 4.1–12.1%.

*Model fitting.* Let y = (y_{t}_{1},...,y* _{tn}*) 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,

*t*, β), where A = (A

_{on}_{1},...,A

*),*

_{J}*t*= (

_{on}*t*

_{on,}_{1},...,

*t*), and β = (β

_{on,J}*, β*

_{I}*, β*

_{C}*), and*

_{f}*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) The –2log likelihood of the model is (5) where

**H**(

*t*

_{1}) = [H

_{1}(

*t*

_{1}),H

_{2}(

*t*)]. Because we estimate H

_{1}_{1}(

*t*

_{1}) = H

_{2}(

*t*

_{1}) = y(

*t*

_{1}), there are

*n*– 1 terms in the –2log likelihood instead of

*n*. We have found that with these assumed values for H

_{1}(

*t*

_{1}) and H

_{2}(

*t*

_{1}), we could reliably estimate the other parameters for all of the data series for all of the subjects. The maximum likelihood estimate of Θ, , is computed by minimizing

*Eq. 5*by use of Newton's method algorithm that we devised previously (3). To assess model goodness-of-fit, we plot for each subject the model fit with the original data, and we report the Akaike's Information Criteria (AIC) (28), defined as –2log

*f*[y|Θ,

**H**(

*t*)] – 2(2

_{1}*J*+ 4). To assess the sensitivity of the parameter estimation to the choice of initial guesses, multiple combinations of starting values for β

*, β*

_{I}_{C}, and β

*were chosen: all combinations of starting values of 1, 3, 5, and 7 h*

_{f}^{–1}for each of the three parameters were used, for a total of 64 different starting conditions for the model fit. The β-values were bounded between 0.05 and 8 h

^{–1}, corresponding to half-lives between 2 and 360 min. The initial estimates of the times of onset of the secretory events and, hence, a provisional guess at the number of secretory events were determined graphically. The maximum likelihood estimates of A and μ were computed in closed form at each step of Newton's procedure given estimates of β and

*t*as in

_{on}*Eq. 3*.

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 1–3* 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 estimation and goodness-of-fit analysis.* For each subject, the GH data collected during the 8-h sleep episode and the 16-h CR episode were fit separately. The model described each subject's data well, as evidenced by the goodness-of-fit statistics and the graphs of the fits (Fig. 2). The values for the rate constants β* _{I}*, β

*, and β*

_{C}*are reported in terms of their approximate half-lives, τ*

_{f}*, τ*

_{I}*, and τ*

_{C}*, where, for example, τ*

_{f}*= log(2)/β*

_{I}*.*

_{I}

*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 τ

*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*

_{I}*t*-test

*P*= 0.53). The parameter estimates of τ

*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 (*

_{C}*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 A* _{t}* 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 · l

^{–1}· h

^{–1}during the sleep episode compared with 34.2 ng · l

^{–1}· h

^{–1}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

We have developed a statistical model of diurnal variation in GH serum levels based on the physiological properties of the hormone. We have shown that this model can be fit to experimental GH data and used to estimate simultaneously the kinetics parameters, the locations and amplitudes of the secretory events, as well as their associated standard errors (12, 15). Our analysis approach is an extension of our previous work, in which we developed linear differential equation models to describe diurnal variation in plasma melatonin and cortisol (4, 5).

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

This study was supported in part by National Institutes of Health (NIH) Grants R01-AR-43130 (to G. K. Adler) and M01-RR-20635 (to the BWH GCRC), K01-AG-00661 (to E. B. Klerman), and NASA Cooperative Agreement NCC9–58 with the National Space and Biomedical Research Institute, and NIH Grants R01-GM-53559 and K02-MH-61637 to E. N. Brown.

Note: The Windows-based program to fit our GH model to experimental data is available upon request from Dr. Klerman.

## Acknowledgments

We thank the patients, nurses, and technicians of the Brigham and Women's Hospital General Clinical Research Center (BWH GCRC) for their efforts.

## 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.

- Copyright © 2003 by American Physiological Society