A stochastic differential equation model of diurnal cortisol patterns

Emery N. Brown1, Patricia M. Meehan1, and Arthur P. Dempster2

1 Neuroscience 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 02114; and 2 Department of Statistics, Harvard University, Cambridge, Massachusetts 02138


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

Circadian modulation of episodic bursts is recognized as the normal physiological pattern of diurnal variation in plasma cortisol levels. The primary physiological factors underlying these diurnal patterns are the ultradian timing of secretory events, circadian modulation of the amplitude of secretory events, infusion of the hormone from the adrenal gland into the plasma, and clearance of the hormone from the plasma by the liver. Each measured plasma cortisol level has an error arising from the cortisol immunoassay. We demonstrate that all of these three physiological principles can be succinctly summarized in a single stochastic differential equation plus measurement error model and show that physiologically consistent ranges of the model parameters can be determined from published reports. We summarize the model parameters in terms of the multivariate Gaussian probability density and establish the plausibility of the model with a series of simulation studies. Our framework makes possible a sensitivity analysis in which all model parameters are allowed to vary simultaneously. The model offers an approach for simultaneously representing cortisol's ultradian, circadian, and kinetic properties. Our modeling paradigm provides a framework for simulation studies and data analysis that should be readily adaptable to the analysis of other endocrine hormone systems.

circadian rhythms; deconvolution; kinetics; sensitivity analysis; ultradian rhythms


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

CORTISOL IS A STEROID HORMONE that in humans is primarily responsible for regulating metabolism and the body's response to inflammation and stress. The 24-h plasma profile of cortisol is comprised of episodic release of 15-21 secretory events whose magnitudes vary in a regular diurnal pattern (30, 52, 53, 56, 57). Plasma levels of the hormone are lowest from 2000 to 0200, climb rapidly through the late night and predawn hours, reach a maximum between 0800 and 1000, and decline throughout the course of the day into the evening (56, 57). Findings in human studies suggest that the circadian pacemaker governs a significant component of the diurnal variation in plasma cortisol levels (13, 53). For this reason, cortisol is often used along with core temperature and plasma melatonin levels to study the properties of the human circadian system (5, 6, 13, 14).

An accurate description of cortisol's diurnal patterns requires a representation that models the relation between the circadian and noncircadian components underlying the diurnal variation. Biochemical and physiological evidence from human investigations suggests that diurnal variation in plasma cortisol levels is governed primarily by the following three factors: 1) ultradian timing of the cortisol secretory episodes (24, 30, 52, 53, 56, 57); 2) circadian control of cortisol secretory amplitudes (24, 52, 53); and 3) the kinetics of cortisol synthesis in the adrenal glands and infusion into (20, 22, 24, 53) and clearance from the plasma by the liver. There is also measurement error variation as a result of the cortisol immunoassay (30, 53).

Cortisol secretion is under the control of the hypothalamic-pituitary-adrenal axis (Fig. 1). In the hypothalamus, communication between the suprachiasmatic nuclei, the site of the circadian pacemaker, and the paraventricular nuclei, the site of corticotropin-releasing hormone (CRH), occurs most probably by direct neural connections (40). CRH, secreted from the paraventricular nuclei into the hypophyseal portal blood vessels, is the principal hypothalamic factor responsible for inducing release of ACTH from the anterior pituitary (41, 49). The frequency of cortisol secretory events by the adrenal gland is tightly coupled to the episodic release of ACTH from the anterior pituitary (30). The amount of cortisol released in each secretory episode appears to be regulated primarily by changes in the amplitude rather than the frequency of the secretory episodes (53). This amplitude modulation is believed to be controlled by the circadian pacemaker through modulation of ACTH release (17, 52).


View larger version (28K):
[in this window]
[in a new window]
 
Fig. 1.   Compartment representation of the 2-dimensional linear stochastic differential equation of plasma cortisol levels. Arrows, direction in which neural and humoral signals are propagated through the hypothalamic-pituitary-adrenal (HPA) axis; +, excitatory interactions; -, inhibitory interactions. Direct neural connections between the suprachiasmatic nucleus (SCN), site of the circadian pacemaker, and the paraventricular nucleus (PVN) are most likely responsible for circadian modulation of PVN release of corticotropin-releasing hormone (CRH). CRH stimulates anterior pituitary release of ACTH, which induces the adrenal gland to synthesize and release cortisol. The adrenal compartment (Eq. 7) includes all of the HPA axis elements within the dotted line. The plasma compartment (Eq. 8) is where the diurnal cortisol rhythm is observed. The rate of cortisol movement from the adrenal gland to the plasma compartment is governed by the infusion rate constant (beta I), whereas the clearance rate constant (beta C) governs cortisol removal from the plasma compartment by the liver. Cortisol exerts a negative feedback effect at the levels of the hypothalamus and the anterior pituitary. This negative feedback is not part of our current model. See Glossary for other definitions.

Because no cortisol is stored in the adrenal gland, initiation of cortisol secretory episodes by ACTH is due to induction of de novo cortisol synthesis from cholesterol by the G protein-coupled receptor-mediated increase in cholesterol desmolase activity and transcription of genes encoding the enzymes required to synthesize cortisol (57). After synthesis, cortisol diffuses into the circulation where ~3-10% of the hormone is free and the remainder is transported bound to either albumin or cortisol-binding globulin (57). Cortisol is absorbed from the plasma in various tissues where it executes its regulatory functions as a steroid hormone. In the hypothalamus and the anterior pituitary, cortisol inhibits, respectively, release of CRH and ACTH by negative feedback (Fig. 1) (57). The hormone is cleared by the liver through reduction of the A ring in the steroid backbone followed by conjugation with glucuronic acid to form several water-soluble compounds that are excreted in the urine (20). The kinetics of cortisol infusion into and clearance from the plasma have been empirically classified as first order (22, 53). The biochemistry of cortisol suggests that a minimum of two compartments must be considered to represent its diurnal pattern (Fig. 1).

Reported concentrations of plasma cortisol depend critically on the reliability of the immunoassay used to measure it. The minimal detectable concentration of the cortisol immunoassays has been reported as 0.5 µg/dl with both intra- and interassay percent coefficients of variation ranging from 5 to 10% (8, 30, 53). The coefficients of variation change appreciably with the number of replicates assayed per sample. Therefore, the immunoassay error represents an important source of measurement variation and must be included explicitly in a model of the hormone's diurnal variation (5, 7).

A complete mathematical description of diurnal cortisol variation should include all of the neural and humoral feedforward and feedback linkages between the hypothalamus, anterior pituitary, and the adrenal gland (40), as well as the effects of exogenous factors such as sleep state, stress, and meals (28, 46). Simultaneous measurement of the variables necessary to specify correctly all components of such a model is not possible in humans and is only possible to a limited extent in other species. Instead of attempting to specify all of these relations, an alternative approach is to use known physiology to define a minimal model of the features necessary to describe the observed diurnal patterns. If the model is sufficiently parsimonious, then the parameters can be estimated from the experimental data. If the minimal model can be estimated for individual subjects, then this information can be used to define individual and population differences in normal and diseased conditions. The minimal model may also help define and set limits on some of the structures in the more comprehensive model. Specifying a minimal model of cortisol's diurnal variation is the approach we take in this report.

Current data analysis methods and mathematical models of cortisol have described some subsets of the three physiological factors underlying cortisol's diurnal variation. None characterizes the relation among the essential components in a single equation system. These models are primarily deterministic in structure and do not take account of the stochastic features of the hormone's diurnal pattern, which cannot be attributed to measurement error from the immunoassay. We formulate a minimal stochastic differential equation model based on these accepted physiological properties of cortisol and use it to describe diurnal variation in the hormone's plasma levels. We determine the model parameters from previous studies of the hormone and show that the joint physiological range assumed by these parameters may be succinctly summarized by a multivariate Gaussian probability density. We establish the plausibility of the model in a series of simulations and describe the relation between our model and current quantitative methods used to study diurnal cortisol patterns.

Glossary


µ(t)   Average circadian amplitude function at time t
Aj   Total amount of cortisol in a secretory event initiated at time uj (µg/dl)
A   Vector of secretory event amplitudes
[A1, A2,..., AN(T)]T
cr   Coefficient of the rth cosine harmonic of the mean circadian amplitude function
dr   Coefficient of the rth sine harmonic of the mean circadian amplitude function
H1(t)   Cortisol concentration in the adrenal glands at time t
H2(t)   Cortisol concentration in the plasma space at time t
N(t)   Number of secretory events occurring in an interval (0,t]
uj   Time of the jth secretory event
u   Vector of secretory event times [u1,u2,...,uN(t)]
wk   Waiting until the kth secretory event given the time of the k - 1st secretory event
yt   Cortisol concentration measured in the plasma at time t (µg/dl)
 beta C   Rate constant for clearance of cortisol from the plasma (min-1)
 beta I   Rate constant for infusion of cortisol from the adrenal glands into the plasma (min-1)
 gamma A   Circadian amplitude function coefficient of variation
 lambda 1   Location parameter of the gamma waiting time probability density
 lambda 2   Scale parameter of the gamma waiting time probability density
µw   Mean waiting time between secretory events
 sigma 2µ(t)   Variance of the circadian amplitude function at time t
 sigma 2w   Variance of the intersecretory event times
 Sigma beta    Covariance matrix of the cortisol infusion and clearance rate constants
 Sigma xi    Covariance matrix of the mean circadian function parameters
 Sigma lambda    Covariance matrix of intersecretory event waiting time parameters
 theta    Vector of model parameters:
theta  = [beta I,beta c,lambda 1,lambda 2,c0,c1,d1,c2,d2]
 xi    Vector of the mean circadian function parameters: xi  = [c0,c1,d1,c2,d2]


    METHODS
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

Model formulation. Given an observation interval (0,T], we define at time t, H1(t), the concentration of cortisol in the adrenal gland, and H2(t) the concentration of cortisol in the plasma. We let N(t) be the counting process that defines the number of secretory events occurring in an interval (0,t) for t is in  (0,T]. We let uj, j = 1, ... , N(T) be the times in (0,T] at which secretory episodes are initiated, and we let wk, k = 1, ... , N(T) be the waiting times between secretory events. It follows from the definitions of uj and N(t) that
u<SUB>j</SUB>=<LIM><OP>∑</OP><LL>k=1</LL><UL>j</UL></LIM> w<SUB>k</SUB> (1)

N(t)=<LIM><OP>∫</OP><LL>0</LL><UL>t</UL></LIM> d<IT>N</IT>(<IT>u</IT>) (2)
where dN(u) indicates a change in the counting process at u. Heuristically speaking, dN(u) is 1 if there is a secretory event at u, and 0 otherwise. We assume the value for wk is a renewal process defined as independent, identically distributed gamma random variables whose probability density has location and shape parameters lambda 1 and lambda 2, respectively (23). We use the gamma probability density because it is a two-parameter model, which offers a flexible description of both the mean and variance of the intersecretory event interval. The expected value and variance of wk are respectively
&mgr;<SUB>w</SUB>=&lgr;<SUB>1</SUB>/&lgr;<SUB>2</SUB> (3)

&sfgr;<SUP>2</SUP><SUB>w</SUB>=&lgr;<SUB>1</SUB>/&lgr;<SUP>2</SUP><SUB>2</SUB> (4)
Let Aj, j = 1, ..., N(T) be the amplitude or total amount of hormone contained in the secretory event initiated at time uj. We assume that the magnitude of the Aj values is a random variable modulated in a time-dependent manner by the circadian system. To represent the stochastic nature of this modulation, we define the two-harmonic mean circadian amplitude function µ(t) and variance function sigma <SUB>&mgr;(<IT>t</IT>)</SUB><SUP>2</SUP> as
&mgr;(t)=c<SUB>0</SUB>+<LIM><OP>∑</OP><LL>r=1</LL><UL>2</UL></LIM> <FENCE>c<SUB>r</SUB> cos <FENCE><FR><NU><IT>2&pgr;</IT>rt</NU><DE><IT>24</IT></DE></FR></FENCE><IT>+d<SUB>r</SUB> </IT>sin <FENCE><FR><NU><IT>2&pgr;</IT>rt</NU><DE><IT>24</IT></DE></FR></FENCE></FENCE> (5)

&sfgr;<SUP>2</SUP><SUB>&mgr;(t)</SUB>=&ggr;<SUP><IT>2</IT></SUP><SUB>A</SUB><IT>&mgr;</IT>(<IT>t</IT>) (6)
where gamma A > 0 is a coefficient of variation. We assume that each At is a Gaussian random variable with mean µ(t) and variance sigma <SUB>&mgr;(<IT>t</IT>)</SUB><SUP>2</SUP>.

Let beta I be the infusion constant governing the rate at which cortisol enters the plasma from the adrenal gland and beta C be the clearance parameter describing the rate at which cortisol is cleared from the plasma by the liver. We assume the kinetics of the infusion and clearance processes to be first order (Fig. 1).

The rate of change in the concentration of adrenal cortisol is equal to the rate of synthesis minus the rate of infusion from the adrenal gland into the plasma. Similarly, the rate of change in the concentration of plasma cortisol is equal to the rate of infusion of cortisol into the plasma from the adrenal minus the rate of its clearance from the plasma by the liver. Under the assumption of first-order kinetics for cortisol infusion and clearance, the rates of change in cortisol concentration in the adrenal gland and plasma may be written as
<FR><NU>d<IT>H<SUB>1</SUB></IT>(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR><IT>=</IT>−<IT>&bgr;</IT><SUB>I</SUB><IT>H<SUB>1</SUB></IT>(<IT>t</IT>)<IT>+A<SUB>t</SUB></IT>d<IT>N</IT>(<IT>t</IT>) (adrenal gland) (7)

<FR><NU>d<IT>H<SUB>2</SUB></IT>(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR><IT>=</IT><IT>&bgr;</IT><SUB>I</SUB><IT>H<SUB>1</SUB></IT>(<IT>t</IT>)<IT>−&bgr;</IT><SUB>C</SUB><IT>H<SUB>2</SUB></IT>(<IT>t</IT>) (plasma) (8)
The solution to this equation system on the interval (t0,t] is
<IT>H</IT>(<IT>t</IT>)<IT>=T</IT>(<IT>t−t<SUB>0</SUB></IT>)<IT>H</IT>(<IT>t<SUB>0</SUB></IT>)<IT>+</IT><LIM><OP>∫</OP><LL><IT>t<SUB>0</SUB></IT></LL><UL><IT>t</IT></UL></LIM><IT> T</IT>(<IT>t−u</IT>)<IT>G</IT>(<IT>u</IT>)d<IT>N</IT>(<IT>u</IT>) (9)
where H(t) = [H1(t),H2(t)]', G(t) = [At,0]' and
<IT>T</IT>(<IT>t−u</IT>)<IT>=</IT><FENCE><AR><R><C>e<SUP>−&bgr;<SUB>I</SUB>(<IT>t−u</IT>)</SUP></C><C>0</C></R><R><C>&kgr;[e<SUP>−&bgr;<SUB>C</SUB>(<IT>t−u</IT>)</SUP><IT>−e</IT><SUP><IT>−&bgr;</IT><SUB>I</SUB>(<IT>t−u</IT>)</SUP>]</C><C>e<SUP>−&bgr;<SUB>C</SUB>(<IT>t−u</IT>)</SUP></C></R></AR></FENCE> (10)
where kappa  = beta I/(beta I - beta C).

To complete the model description, we assume that during a time interval of length T we collect n blood samples, which are assayed for cortisol. We let yt1, ... , ytn denote the concentration of cortisol at times 0 < t1 < t2 <...< tn - 1 tn<= T and assume the yti satisfy the equation
y<SUB>t<SUB>i</SUB></SUB>=H<SUB><IT>2</IT></SUB>(<IT>t<SUB>i</SUB></IT>)<IT>+&egr;</IT><SUB><IT>t<SUB>i</SUB></IT></SUB> (11)
where epsilon ti are independent, zero mean Gaussian random variables with variance sigma 2epsilon ti defined primarily by the immunoassay measurement error.

Taken together, Eqs. 9 and 11 describe a two-dimensional state space model; the former is the state equation and the latter is the observation equation. The state vector of this model is H(t), and its components are the concentrations of cortisol in the adrenal gland and in the plasma. Equations 7 and 8 define how these variables change over time. It is a stochastic differential equation system because its forcing term AtdN(t) is a random amount of cortisol input at a random time t. The solution to Eqs. 7 and 8 is the stochastic convolution integral in Eq. 9. This model also belongs to a class of stochastic processes known as filtered point processes (45) in which the counting process is defined implicitly by the gamma probability model we assumed for the intersecretory event times. The mark process, At, is a Gaussian random variable whose mean (Eq. 5) and variance (Eq. 6) are periodic functions of time.

Characterizing the probability density of the model parameters. We develop a two-step approach to simulating the model in Eqs. 1-11. We represent the joint probability density of the observed plasma cortisol levels, the secretory amplitudes, the secretory event times as
f(y,A,u<IT>‖&thgr;</IT>)<IT>=f</IT>(<IT>y‖A,u,&thgr;</IT>)<IT>f</IT>(<IT>A‖u,&thgr;</IT>)<IT>f</IT>(<IT>u‖&thgr;</IT>) (12)
where y = [yt1, ... , ytn], A = [A1, ... , AN(T)], u = [u1, ... , uN(T)], and theta  = [beta I,beta C,lambda 1,lambda 2,c0,c1,d1,c2,d2]; f(y|A,u,theta ) is the multivariate Gaussian probability density defined in Eq. 11; f(A|u,theta ) is the joint probability density of the secretory event amplitudes whose mean and variance follow implicitly from Eqs. 5 and 6; and f(u|theta ) is the joint probability density of the secretory event times defined implicitly by the gamma probability density for the intersecretory event times given in Eqs. 1-4.

If we assume that different individuals have different values of theta , then simulation of plasma cortisol data from Eq. 12 with theta  fixed is equivalent to studying variation in the hormone's plasma levels for an individual. Therefore, we consider the expanded model
f(y,A,u,&thgr;)=f(y,A,u‖&thgr;)f(&thgr;) (13)
where f(theta ) denotes the probability density that summarizes the between-individual variation in theta . We simulate our cortisol model by first drawing theta  from f(theta ) and then simulating y, A, and u given theta  from f(y,A,u|theta ). This two-step approach is equivalent to simulating between- and within-individual variation in plasma cortisol patterns. The between-individual variation is due to interindividual differences in ultradian, circadian, and kinetic properties of cortisol and is summarized by f(theta ). The within-individual variation is summarized by f(y,A,u|theta ). A standard approach for assessing the sensitivity of a model to the choice of parameters is to compare the results of simulations from the model as the parameters are changed one at a time. Along with characterizing between-subject variation in plasma cortisol levels, this two-step simulation approach suggested by Eq. 13 allows all model parameters to be changed simultaneously.

We assume that f(theta ) is a multivariate Gaussian probability density, which has the form
f(&thgr;)=f(&lgr;)f(&xgr;)f(&bgr;) (14)
where lambda  = [lambda 1,lambda 2], xi  = [c0,c1,d1,c2,d2], and beta  = [beta I,beta C]. The probability densities f(lambda ), f(xi ), and f(beta ) summarize, respectively, the uncertainty in the ultradian, the circadian, and the kinetic parameters. We make the Gaussian assumption because a Gaussian probability density is completely defined by its mean vector and covariance matrix. We will derive its mean vector and covariance matrix from published reports that have quantitatively modeled diurnal cortisol patterns. We limited our analysis to those investigations in which subjects were monitored for at least 24 h and intersample intervals were 20 min or less. We excluded the study of Mortola and colleagues (35), which analyzed diurnal cortisol patterns in healthy female subjects, because they did not report the phase of the menstrual cycle during which the subjects were studied. It is now appreciated that the character of women's circadian rhythms is different at different phases of the menstrual cycle (42). We summarize below the analysis used to develop each of these three components of f(theta ).

Ultradian parameter probability density: f(lambda ). Linkowski and colleagues (30) reported the number of cortisol secretory events in 24 h for each of seven subjects. We computed for each subject the average intersecretory event interval by dividing 24 h by the number of secretory events. Veldhuis and colleagues (53) reported the individual average and SD of the intersecretory event intervals for six subjects. Assuming equal weighting of each of these 13 subjects, we computed µw and sigma <SUB><IT>w</IT></SUB><SUP>2</SUP>, the estimated between-subject mean and variance, respectively. On the basis of the assumption that the intersecretory event interval process is described by a gamma probability density, the method of moments estimate of the mean values of lambda 1 and lambda 2 can be computed as (23)
&lgr;<SUB>1</SUB>=&mgr;<SUP>2</SUP><SUB>w</SUB>/&sfgr;<SUP>2</SUP><SUB>w</SUB> (15)

&lgr;<SUB>2</SUB>=&mgr;<SUB>w</SUB>/&sfgr;<SUP>2</SUP><SUB>w</SUB>
We derive the covariance matrix for lambda 1 and lambda 2 from the data presented by Veldhuis and colleagues, since they provided estimates of µw and sigma <SUB><IT>w</IT></SUB><SUP>2</SUP> for six individual subjects (53). We used Eq. 15 to convert each subject's estimate of µw and sigma <SUB><IT>w</IT></SUB><SUP>2</SUP> into estimates of lambda 1 and lambda 2. Using the six pairs of lambda 1 and lambda 2 along with estimates of the means of lambda 1 and lambda 2 obtained by combining the data from Refs. 30 and 53, we computed the sample covariance matrix using the standard formula (1). The estimates of the mean vector and covariance matrix are
&mgr;<SUB>&lgr;</SUB>=<FENCE><AR><R><C>54</C></R><R><C>39</C></R></AR></FENCE> &Sgr;<SUB>&lgr;</SUB>=<FENCE><AR><R><C>246</C><C>135</C></R><R><C>135</C><C>100</C></R></AR></FENCE> (16)
Hence, we take f(lambda ) to be the bivariate Gaussian probability density whose mean vector and covariance matrix are given in Eq. 16.

Kinetic parameter probability density: f(beta ). Hellman et al. (22), Weitzman et al. (56), Linkowski et al. (30), and Veldhuis et al. (53) reported individual clearance half-lives for two, seven, seven, and six subjects, respectively. Jusko et al. (24) reported the average half-life and its SD for seven subjects. We combined the 22 subjects from the studies in Refs. 22, 57, 30, and 53 to compute the mean and SD of the clearance half-life. We converted these mean and SD estimates into estimates of the mean and variance of beta C using the delta method formula (38)
E[<IT>&bgr;</IT><SUB>C</SUB>]<IT>=</IT>−<FR><NU>log<IT> 2</IT></NU><DE><A><AC><IT>t</IT></AC><AC>&cjs1171;</AC></A><SUB><IT>1/2</IT></SUB></DE></FR><IT>+</IT><FENCE><FR><NU>log 2s<SUP><IT>2</IT></SUP><SUB><IT>t<SUB>1/2</SUB></IT></SUB></NU><DE><A><AC><IT>t</IT></AC><AC>&cjs1171;</AC></A><SUP><IT>3</IT></SUP><SUB><IT>1/2</IT></SUB></DE></FR></FENCE> (17)

Var [<IT>&bgr;</IT><SUB>C</SUB>]<IT>=</IT><FENCE><FR><NU>log<IT> 2</IT></NU><DE><A><AC><IT>t</IT></AC><AC>&cjs1171;</AC></A><SUP><IT>2</IT></SUP><SUB><IT>1/2</IT></SUB></DE></FR></FENCE><SUP><IT>2</IT></SUP>s<SUP><IT>2</IT></SUP><SUB><IT>t<SUB>1/2</SUB></IT></SUB> (18)
where <A><AC><IT>t</IT></AC><AC>&cjs1171;</AC></A><SUB>1/2</SUB> and s<SUB><IT>t</IT><SUB>1/2</SUB></SUB><SUP>2</SUP> are the mean and variance of the half-life estimates, respectively. We also converted the estimates of the average and SD of the half-life from the seven subjects in Ref. 23 to mean and variance estimates of beta C using Eq. 16. We combined the estimates of the mean of beta C and the variance of beta C from the 22 subjects with the ones computed from the 7 subjects by taking a weighted average. This yielded mean and variance estimates for beta C of 0.645 min-1 and 0.015 min-2.

The analyses of Veldhuis et al. (53) provide information on beta I and its variability. These authors reported an average infusion half-life of 16 min with a SD of 0.61 min in six subjects. Using the delta method formula in Eqs. 17 and 18, we computed the mean and variance of beta I to be 2.71 min-1 and 0.074 min-2, respectively. Although the mean estimate seemed reasonable, we suspected that this variance estimate may understate the uncertainty in beta I because it was based on only six subjects. Therefore, because the sample variance estimate is approximately distributed as a chi 2 random variable with five degrees of freedom, and the 99th percentile of this probability distribution is 0.669, we used this value as the estimate of the variance of beta I. Because we identified no reports that jointly estimated beta C and beta I, we assumed the covariance between beta C and beta I was zero and represented the joint probability density of these two parameters as the bivariate Gaussian density whose mean vector and covariance matrix are, respectively
&mgr;<SUB>&bgr;</SUB>=<FENCE><AR><R><C>2.71</C></R><R><C>0.646</C></R></AR></FENCE> &Sgr;<SUB>&bgr;</SUB>=<FENCE><AR><R><C>0.669</C><C>0</C></R><R><C>0</C><C>0.015</C></R></AR></FENCE> (19)

Circadian parameter probability density: f(xi ). We estimated the probability density for the circadian parameters by the analysis of plasma cortisol data collected from four healthy male subjects studied on 24-40 h constant routines (13). We fit an approximate form of the model in Eqs. 9 and 11 to each subject's plasma cortisol data by Bayesian Markov chain Monte Carlo methods (8, 9, 32). Our implementation of this algorithm provided an approximate Gaussian probability density for the circadian parameters for each subject. The average of the mean vectors and covariance matrices across the four subjects was
&mgr;<SUB>&xgr;</SUB>=<FENCE><AR><R><C>6.10</C></R><R><C>−<IT>4.75</IT></C></R><R><C>3.93</C></R><R><C>−<IT>3.76</IT></C></R><R><C>−<IT>2.53</IT></C></R></AR></FENCE> (20)

<IT>&Sgr;<SUB>&xgr;</SUB>=</IT><FENCE><AR><R><C>1.937</C><C>0.423</C><C>−<IT>0.095</IT></C><C>0.573</C><C>−<IT>0.214</IT></C></R><R><C></C><C>1.253</C><C>−<IT>0.197</IT></C><C>0.514</C><C>0.296</C></R><R><C></C><C></C><C>0.631</C><C>−<IT>0.254</IT></C><C>−<IT>0.125</IT></C></R><R><C></C><C></C><C></C><C>1.276</C><C>−<IT>0.123</IT></C></R><R><C></C><C></C><C></C><C></C><C>0.509</C></R></AR></FENCE>
where the lower triangle is omitted because of symmetry. Therefore, we took f(xi ) to be the Gaussian probability density whose mean vector and covariance matrix are defined in Eq. 20. The value of gamma A, the coefficient of variation, could not be determined from previous studies. Therefore, we simulate our model for the following three values of gamma A: 0.1, 0.3, and 0.5.

Joint probability density of the model parameters: f(theta ). Because no single study provided information about all of the parameters in our model, the covariances between the kinetic, ultradian, and circadian parameters could not be estimated. Therefore, we assumed these covariances to be zero. Combining Eqs. 16, 19, and 20, it follows that the joint probability density of the parameters is the multivariate Gaussian density whose mean and covariance matrix are
&mgr;<SUB>&thgr;</SUB>=<FENCE><AR><R><C>&mgr;<SUB>&lgr;</SUB></C></R><R><C>&mgr;<SUB>&bgr;</SUB></C></R><R><C>&mgr;<SUB>&xgr;</SUB></C></R></AR></FENCE> &Sgr;<SUB>&thgr;</SUB>=<FENCE><AR><R><C>&Sgr;<SUB>&lgr;</SUB></C><C>0</C><C>0</C></R><R><C>0</C><C>&Sgr;<SUB>&bgr;</SUB></C><C>0</C></R><R><C>0</C><C>0</C><C>&Sgr;<SUB>&xgr;</SUB></C></R></AR></FENCE> (21)

Immunoassay error. As stated in the INTRODUCTION the intra-assay coefficients of variation for the cortisol immunoassay measurement error range between 0.05 and 0.10 when multiple replicates of the blood sample are assayed and may be as high as 0.15 when a single replicate is assayed (7, 8). In either case, the probability density of the immunoassay error can be reasonably well approximated by a Gaussian density (5, 7). Under the assumption that each blood specimen is assayed in duplicate, we take the coefficient of variation of the errors in Eq. 11 to be 0.07.

Simulation algorithm of diurnal cortisol patterns. Given measurement times 0 < t1 < t2 <...< tn-1 < tn <=   T, the simulation algorithm for the cortisol model defined in Eqs. 1-9 proceeds through the following seven steps
(22)
All simulations were carried out using the Splus programming language.


    RESULTS
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

Figure 2 gives a graphic illustration of the steps in the simulation algorithm given in Eq. 22. In this simulation and all subsequent ones, we assume that the initial cortisol measurement is made at midnight and therefore that the initial concentrations of cortisol in the adrenal gland and the plasma compartments are zero. An illustration of the secretory event times and amounts of cortisol produced in the adrenal compartment as generated by steps 1-4 of the algorithm are shown in Fig. 2A. The circadian modulation of the secretory process is shown in the variation of the amplitudes of the secretory events. Ultradian variation in secretory event timing is less pronounced because the values of the mean and variance of the intersecretory event times in this example are 72 and 9 min, respectively.


View larger version (14K):
[in this window]
[in a new window]
 
Fig. 2.   Graphic summary of a 24-h realization of the cortisol simulation algorithm in Eq. 22 for a given value of theta . A: times and amounts of cortisol produced in the adrenal compartment in the 24-h period (steps 2-4 of the algorithm). B: time course of the impulse response function of a secretory event initiated at time 0 as defined by Eq. 9. C: evaluation of Eq. 9 to obtain the true plasma cortisol levels for the secretory episodes in A (step 5). D: observed plasma cortisol level is composed of the true plasma cortisol in C plus the error due to the cortisol immunoassay (steps 6 and 7).

Step 5 of the algorithm integrates the stochastic differential model in Eqs. 7 and 8 to obtain the solution in Eq. 9 (Fig. 2C). The impulse response function for this system, as shown in Fig. 2B, is derived directly from Eq. 9. That is, for a given secretory event, i.e., a uk and an Ak, the infusion and clearance parameters govern the fraction of Ak that appears in the plasma at a given time after initiation of the secretory event. Within ~30 to 40 min is when the largest fraction of the Ak enters the plasma. By 6 h, all of the cortisol in this secretory event has reached the plasma. In step 6 of the algorithm, Gaussian noise based on the properties of the cortisol immunoassay is added to the true plasma cortisol levels. The coefficient of variation is set at 0.07, making the variance proportional to the square of the true plasma cortisol level. Finally, in step 7, the immunoassay noise is added to the simulated true plasma level (Fig. 2D). The difference between C and D in Fig. 2, shows the effect of the immunoassay error on the measured plasma cortisol levels.

To simulate the cortisol model, we chose the observation interval T to be 48 h, the time between measurements to be 10 min, and gamma A to be 0.1 (Fig. 3), 0.3 (Fig. 4), and 0.5 (Fig. 5). As stated above, all simulations start at midnight when both the initial adrenal and plasma cortisol concentrations are assumed to be zero. [Initial condition estimates for times other than midnight can be obtained by beginning at midnight (zero initial condition) and simulating the model over several days. The value of the cortisol series taken at the time of interest several days into the simulation should be a reasonable initial guess.] We drew different valves of theta  from f(theta ), and for each value of gamma A, we simulated six 48-h cortisol series. Figures 3, 4, and 5 show the same value of theta . That is, for example, Figs. 3A, 4A, and 5A were all simulated with the same random draw of theta . Figures 3-5 differ only in their choice of gamma A. Therefore, Figs. 3-5 each represent 48 h of plasma cortisol levels for the same individual with three different amplitude coefficients of variation. Different panels in Figs. 3-5 represent different individuals. Together, the 18 panels in Figs. 3-5 represent both between- and within-subject variation in plasma cortisol diurnal patterns.


View larger version (18K):
[in this window]
[in a new window]
 
Fig. 3.   Forty-eight-hour simulations of the population variation in cortisol plasma levels with theta  chosen randomly from f(theta ) in A-F and gamma A = 0.1. Each panel (A-F) is simulated with a different value of theta  chosen randomly from f(theta ). Corresponding panels in Figs. 4, 5, and 6 were simulated with the same values of theta  but a different value of gamma A.



View larger version (18K):
[in this window]
[in a new window]
 
Fig. 4.   Forty-eight-hour simulations of the population variation in cortisol plasma levels with theta  chosen randomly in each panel from f(theta ) and gamma A = 0.3. Each panel (A-F) is simulated with a different value of theta  chosen randomly from f(theta ). Corresponding panels in Figs. 3 and 5 were simulated with the same values of theta  but a different value of gamma A.



View larger version (19K):
[in this window]
[in a new window]
 
Fig. 5.   Forty-eight-hour simulations of the population variation in plasma cortisol levels with f(theta ) chosen randomly in each panel from f(theta ) and gamma A = 0.5. Each panel (A-F) is simulated with a different value of theta  chosen randomly from f(theta ). Corresponding panels in Figs. 3 and 4 were simulated with the same values of theta  but a different value of gamma A.

As expected, the variability in the plasma cortisol levels increases as gamma A increases from 0.1 (Fig. 3) to 0.3 (Fig. 4) to 0.5 (Fig. 5). Between-day variation within subjects can be appreciated by comparing within a panel the first 24 h with the second 24 h. For example, in Fig. 4, A and F, the plasma levels in the first 24 h are higher than those in the second, whereas in Fig. 3A the plasma levels in the second 24 h are higher. For gamma A = 0.1 (Fig. 3), there is less between-day variation in the plasma cortisol levels. The simulated data from each value of gamma A show a good resemblance to the diurnal patterns seen in actual cortisol series (see Fig. 8).

The simulated circadian amplitude function for the panels in Figs. 3 to 5 is shown in Fig. 6, along with the locations and amplitudes of the secretory events. The amplitude functions have the same asymmetric bimodal shape due to the two-harmonic model for this process defined in Eqs. 5 and 6. The largest mode of the circadian function occurs between 0600 and 0900, and the second smaller mode occurs between 2100 and 2200. The amplitude is near zero around midnight for nearly all of the series. A close comparison of the largest mode of the circadian amplitude function with the time of maximum plasma level in the early morning (Figs. 3-5) shows that the latter tends to lag behind the former by ~45-60 min. This lag is due to the infusion and clearance processes, i.e., the impulse response function (Fig. 2B). Because there is infusion of cortisol from the adrenal gland into the plasma, all of the hormone cannot enter the plasma instantaneously. Moreover, concomitant with its entry into the plasma, the cortisol is being metabolized in the liver. The 45- to 60-min lag is consistent with the time for a secretory event to reach its maximum effect on the plasma as shown in Fig. 2B. Although each simulated circadian function has the same general shape and time course, there is significant between-subject variation in its structure. Between days, for a given subject, the average amplitude function is the same because it is periodic.


View larger version (25K):
[in this window]
[in a new window]
 
Fig. 6.   Simulated circadian amplitude functions and pulse locations of the plasma cortisol series in Figs. 3-5. Each panel (A-F) is simulated with a different value of theta  chosen randomly from f(theta ). A-F of the simulated secretory event amplitudes and pulse locations correspond to A-F for the cortisol series in Figs. 3-5. Plasma cortisol concentrations were converted to amounts of cortisol using the volumes of distribution estimates reported in Ref. 53 and normal total daily cortisol production reported in Ref. 19.

To appreciate better within-subject variation in plasma cortisol levels, we simulated an additional set of six cortisol series with theta  fixed at theta 0 and gamma A = 0.3 (Fig. 7). That is, for all six draws, theta  is set equal to theta 0 in step 1 of the simulation algorithm (Eq. 22). These simulated series show that a given subject can have appreciable day-to-day variation in plasma cortisol levels.


View larger version (18K):
[in this window]
[in a new window]
 
Fig. 7.   Forty-eight-hour simulations of the within-individual variation in plasma cortisol levels obtained by simulating the model with theta  fixed at theta 0, the mean of f(theta ), and gamma A = 0.3. Each panel (A-F) represents a different simulation of the cortisol model using theta  = theta 0.

As a qualitative assessment of the ability of our model to reproduce the diurnal patterns in actual cortisol series, we plot in Fig. 8 cortisol data from eight healthy male subjects measured on a 24-h constant-routine protocol (6, 13). The data were collected as part of the study conducted by Dr. Charles A. Czeisler to measure the ability of bright light to shift the phase of the human circadian pacemaker at a fixed phase of the circadian cycle (3). In this study, five groups of either seven or eight subjects received a 5-h light treatment of either 0.03, 10, 180, 1,260, or 9,500 lux. The subjects whose cortisol data are shown in Fig. 8 were in the 9,500-lux group. These cortisol measurements were made at 20-min intervals during the baseline constant (routines the subjects underwent before receiving light therapy). These cortisol series are representative of the data from the 39 subjects studied under this protocol.


View larger version (32K):
[in this window]
[in a new window]
 
Fig. 8.   Actual cortisol series collected from 8 healthy male subjects on a 24-h constant routine. The cortisol samples were taken at 20-min intervals. Each panel (A-H) indicates a different subject.

These experimental cortisol plasma levels are consistent with the asymmetric circadian amplitude modulation defined in Eqs. 5 and 6 and shown in Fig. 6. The rising phase of the cortisol cycles requires 6 h or more, and the decline phase is ~10-12 h. Three different patterns are evident both between and within subjects. These are fine pulsatile activity (Fig. 8, C and F), extended periods of short or low pulsatile activity (Fig. 8, B, D, F, and H), and periods of apparently large pulses (Fig. 8A). Several of the series appear slightly smoother than the simulated series in Figs. 3-5. This may in part be due to the difference in sampling intervals, which is every 20 min for the actual cortisol series and every 10 min in our simulated series. Overall, the qualitative agreement between the model and actual cortisol series is good.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

Development of methods for the analysis of endocrine hormone data is an active area of research. These procedures fall typically into one of four categories: pulse-finding algorithms, deconvolution procedures, Fourier and harmonic regression methods, and approximate entropy (ApEn). Pulse-finding algorithms use local statistical criteria such as Bonferroni bounds, t-tests, or a specific multiple of the immunoassay coefficient of variation to identify in hormone data series the times of secretory events (11, 33, 34, 36, 43, 51, 54). Deconvolution procedures use deterministic, usually linear time-invariant convolution integrals to represent the relation between inputs, secretory event times and amplitudes, and observed output hormone levels (55). The deconvolution analyses use convolution integral models to determine the secretory event times and secretory amplitudes from the observed hormone levels. If the hormone's diurnal pattern has an important circadian component, this is generally extracted in a separate analysis using harmonic regression methods (44, 46, 50, 53). ApEn is a technique that may be used to assess regularity and complexity of any biological series (36). Differences in regularity and complexity detected by ApEn can help characterize differences in physiological and pathological states. With current methods, separate analyses are required to completely analyze circadian, ultradian, and kinetic components of a hormone data series.

The first step toward obviating multiple analyses is to devise a single equation system that captures the principal physiological components in the experimental data. We studied cortisol because much is known about the physiology of its diurnal patterns. Our work builds on the two-compartment deterministic differential equation model of cortisol's diurnal patterns proposed by Jusko et al. (24). The model by Jusko et al. represents circadian input as a simple sine wave and the kinetic process as first order. Secretory events are defined as local percentage changes in plasma hormone levels, and the model is fit to experimental data by using nonlinear least squares analysis. The model of Jusko et al. is, to our knowledge, the first attempt at a single equation system description of the ultradian, circadian, and kinetic properties of cortisol's diurnal patterns. This model does not consider immunoassay uncertainty.

To improve upon this model, we represented the physiological properties summarized by Jusko et al. as a stochastic rather than a deterministic model. We modeled the ultradian process after a gamma renewal process because this probability model is widely used to describe waiting time phenomena (16, 23). From Eqs. 3, 4, 15, and 16, the average intersecretory event time is ~83 min with a SD of 11 min. Hence, the probability of a first secretory event followed within 60 min or less by a second one is small to negligible. The low probability of rapid-succession secretory events may reflect a refractory period of cortisol synthesis induced by negative feedback at the hypothalamus, anterior pituitary, and possibly the adrenal gland. Moreover, the ultradian and circadian processes may also interact in that the degree of feedback inhibition may depend on the amount of hormone released in a secretory episode. For these reasons, the assumptions of independence of the intersecretory event times (renewal process) and of independence between the circadian and ultradian processes are only plausible first approximations to be modified in future versions of the model.

Formulation of the ultradian input to our model is based on the reported high concomitance between ACTH and cortisol secretory episodes and the fact that the adrenal gland maintains no stores of cortisol (29-31, 57). Because these observations suggest that the ultradian process must have a substantial effect on the initiation of cortisol synthesis, the ultradian process defines in our model the times at which the secretory events begin. This representation differs from the model of Jusko et al. in which the ultradian process modulates hormone release. It differs also from the convolution integral model of Veldhuis and Johnson (55) in which the ultradian process identifies the time at which one-half the amount of hormone in a given secretory episode has been released. Whereas the ultradian process may also affect adrenal cortisol release, we have modeled infusion as a first-order kinetic process independent of the ultradian modulation.

The representation of our model's circadian component is based on the finding that the circadian system modulates primarily the amplitude rather than the frequency of the secretory events (53). The shape of the mean circadian waveform deduced by applying a preliminary version of the current model presented (8, 9) to the seven cortisol series in Ref. 56 closely resembles those previously derived using Fourier methods (30, 50) and those shown in Fig. 6. A second harmonic added to the circadian amplitude function accounts simply for the asymmetry in the underlying circadian process. Our decision to model the amplitude process at a given time as a Gaussian random variable is arbitrary. We also studied a frequency-modulation model of the circadian input. None of its simulated outputs resembled actual cortisol data.

Because neither the form nor the magnitude of the amplitude variance could be deduced from previous reports, we represented this parameter as a function of the circadian amplitude. We obtained an apparently realistic data series with coefficients of variation in the range from 0.1 to 0.5. This variability most probably reflects variability in both the ACTH stimuli delivered to the adrenal gland and in the cortisol synthetic response to the ACTH stimuli. Although there is evidence to suggest that cortisol clearance from the plasma may have some degree of diurnal modulation (15), we have found that the hormone's kinetic parameters are well described as a time-invariant first-order kinetic process.

The second step in defining a single equation system for the analysis of cortisol's diurnal patterns is to show that model simulations reliably reproduce experimental data. Our approach builds on the endocrine hormone simulation model of Guardabasso et al. (21). These authors report an autoregressive moving average process in which coefficients are nonlinear functions of the decay parameters. Secretory event times are generated from a discrete Poisson approximation, the amplitudes are sampled from either a uniform or truncated Gaussian probability density, and the observational errors are Gaussian and depend on the immunoassay uncertainty. Our model makes use of two extensions these authors suggested in their discussion: use of a more general renewal process model to describe intersecretory event intervals and a time-dependent circadian modulation. The simulation algorithm proposed by these authors can be greatly simplified and implemented in continuous rather than discrete time by reexpressing their model as in Eqs. 1-11.

Three other stochastic simulation models relate to our work. Straume et al. (47) reported a model for simulating hormone patterns in which intersecretory event times are based on a logistic model, secretory event amplitudes depend logarithmically on secretory event times, and secretion is modeled as a Gaussian-shaped rate process. Diurnal timing is either a square-wave (day-night) or sinusoidal (circadian) process, hormone elimination is either a mono- or biexponential process, and measurement uncertainty is due to assay variability. The authors do not relate the several model components in a single equation system; however, they do show that the model successfully simulates plasma growth hormone levels. Keenan (25) reports a model in which the true hormone level is a one-dimensional stochastic convolution of secretory events from a periodic, inhomogeneous Poisson process and a time-invariant gamma distribution. Integrated Brownian motion noise added to the convolution integral output defines the measured hormone level. Model fitting to simulated data is by nonlinear least squares. This model differs from ours in that the Poisson process governs secretory event times (frequency modulation), whereas the gamma distribution subsumes both the secretory event amplitudes and the kinetic processes in our model. Keenan and Veldhuis (26) report a stochastic differential equation model to simulate the hypothalamic-pituitary-gonadal axis. Unlike our minimal feedforward cortisol model, this model specifies all of the feedforward and feedback neural and humoral linkages in the hypothalamic-pituitary-gonadal axis, including ultradian, circadian, kinetic, and time-delay parameters. This model has the potential to give a nearly complete description of the dynamics of the hypothalamic-pituitary-gonadal system under normal and pathological conditions. However, it requires specification of relations among variables and postulating values of several parameters that are not known or directly measurable. This equation system would require simplification for use in data analysis. None of the other stochastic models of endocrine hormones exploits physiological properties in its design (4, 18, 37).

By summarizing in a probability density the values of the model parameters from previous reports, we characterize their joint physiological range. The model may be used to simulate formally within- and between-subject variation and to perform a sensitivity analysis in which all parameters vary simultaneously. We propose this approach to assess the model's sensitivity to its parameter values as a more appealing alternative to that of varying individual parameters one at a time. Simulated series generated with our model may be used to test current pulse-finding, deconvolution, Fourier, and ApEn methods.

The final step in developing a single system for the analysis of cortisol's diurnal patterns is to establish that the model can be fit to an experimental cortisol time series. In a report currently under preparation, we present a set of methods to fit our cortisol model to experimental data using Bayesian Markov Chain Monte Carlo methods. The joint parameter probability density in Eq. 21 can serve as a prior probability density for the Bayesian analysis. Preliminary work applying these methods has been reported previously (8, 9, 32). In addition to including feedback interactions, other extensions of our cortisol model that we will consider are the effect of factors exogenous to the hypothalamus-pituitary-adrenal axis, such sleep states and meals (28, 46), and simultaneous modeling of cortisol and ACTH.

In summary, we have developed a model of cortisol's diurnal patterns that unifies in a minimal physiologically based stochastic differential equation the principal elements underlying current pulse-finding, deconvolution, and Fourier data analysis methods. Our modeling paradigm provides a framework for both simulation studies and data analysis that should be readily adaptable to the analysis of other endocrine hormone systems. Analysis of physiological systems with stochastic models may help alleviate some presently unexplained discrepancies between experimental data and model descriptions derived from deterministic models.


    ACKNOWLEDGEMENTS

We thank three anonymous referees whose suggestions helped improve the content and presentation in this work; Dr. C. A. Czeisler for use of the cortisol data in Fig. 8; R. Barbieri, Y. Choe, and L. M. Frank for assistance with the data analysis and graphics; and B. Marshall for assistance with manuscript preparation.


    FOOTNOTES

This work was supported by Robert Wood Johnson Foundation Grants 19122 and 23397; National Aeronautics and Space Administration (NASA) Grant NAGW 4061 and NASA Cooperative Agreement NCC 9-58 with the National Space and Biomedical Research Institute; National Institutes of Health Grants 1-R01-GM-53559, 1-P01-AG-09975, 1-R01-AG-06072, 1-R01-MH-45130, and M01-RR-02635; an American Dissertation Fellowship from the American Association of University Women; and Army Research Office Grant DAAL03-91-6-0089.

Address for reprint requests and other correspondence: E. N. Brown, Dept. of Anesthesia and Critical Care, Massachusetts General Hospital, 55 Fruit St., Boston, MA 02114 (E-mail: brown{at}SRLB.mgh.harvard.edu).

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.

Received 21 December 1999; accepted in final form 5 October 2000.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

1.   Anderson, TW. Introduction to Multivariate Statistical Analysis. New York: Wiley, 1984.

2.   Batschelet, E. Circular Statistics in Biology. London: Academic, 1981.

3.   Boivin, DB, Duffy JF, Kronauer RE, and Czeisler CA. Dose-response relationships for resetting of human circadian clock by light. Nature 379: 540-542, 1996[ISI][Medline].

4.   Bolstad, WM. The multiprocess dynamic linear model with biased perturbations: a real time model for growth hormone level. Biometrika 75: 685-692, 1988[ISI].

5.   Brown, EN, Choe Y, Shanahan T, and Czeisler CA. A mathematical model of diurnal variation in human plasma melatonin levels. Am J Physiol Endocrinol Metab 272: E506-E516, 1997[Abstract/Free Full Text].

6.   Brown, EN, and Czeisler CA. The statistical analysis of circadian phase and amplitude in constant-routine core-temperature data. J Biol Rhythms 7: 177-202, 1992[ISI][Medline].

7.   Brown, EN, McDermott T, Bloch KJ, and McCollom AD. Defining the smallest analyte concentration an immunoassay can measure. Clin Chem 42: 893-903, 1996[Abstract/Free Full Text].

8.   Brown, EN, Meehan PM, and Dempster AP. A stochastic differential equation model for simulation and analysis of diurnal cortisol patterns (Abstract). In: Proc 3rd Meet Soc Res Biol Rhythms, 1992.

9.   Brown, EN, Meehan PM, Dempster AP, and Czeisler CA. A mathematical model for the analysis of diurnal cortisol patterns (Abstract). In: Proc 4th Meet Soc Res Biol Rhythms, Amelia Island, FL, 1994.

10.   Butler, JP, Spratt I, O'Dea STL, and Crowley WF. Interpulse interval sequence of LH in normal man essentially constitutes a renewal process. Am J Physiol Endocrinol Metab 250: E338-E340, 1986[Abstract/Free Full Text].

11.   Clifton, DK, and Steiner RA. Cycle detection: a technique for estimating the frequency and amplitude of episodic fluctuations in blood hormone and substrate concentrations. Endocrinology 112: 1057-1064, 1983[Abstract].

12.   Cone, RD, and Mountjoy KG. Molecular genetics of the ACTH and melanocyte-stimulation hormone receptors. Trends Endocrinol Metab 4: 242-247, 1993[ISI].

13.   Czeisler, CA, Allan JS, Strogatz SH, Ronda JM, Sánchez R, Ríos CD, Freitag WO, Richardson GS, and Kronauer RE. Bright light resets the human circadian pacemaker independent of the timing of the sleep-wake cycle. Science 233: 667-671, 1986[ISI][Medline].

14.   Czeisler, CA, Duffy JD, Shanahan TA, Brown EN, Mitchell JF, Dijk DJ, Rimmer DW, Ronda JM, Allan JS, Emens JS, and Kronauer RE. Stability, precision, and near 24 hour period of the human circadian pacemaker older subject. Science 284: 2177-2181, 1999[Abstract/Free Full Text].

15.   DeLacerda, L, Kowarski A, and Migeon CJ. Diurnal variation of the metabolic clearance rate of cortisol. Effect on measurement of cortisol production rate. J Clin Endocrinol Metab 36: 1043-1049, 1973[ISI][Medline].

16.   DeNicolao, G, Guardabasso V, and Rocchetti M. The relationship between rate of venous sampling and visible frequency of hormone pulse. Comput Methods Programs Biomed 33: 145-157, 1990[ISI][Medline].

17.   Désir, D, Van Cauter E, Golstein J, Fang VS, Leclercq R, Refetoff S, and Copinschi G. Circadian and ultradian variations of ACTH and cortisol secretion. Horm Res 23: 302-316, 1980.

18.   Diggle, PJ, and Zeger SL. A non-Gaussian model for time series with pulses. J Am Stat Assoc 84: 354-359, 1989[ISI].

19.   Esteban, NV, Loughlin T, Yergey AL, Zawadzki JK, Booth JD, Winterer JC, and Loriaux DL. Daily cortisol production rate in man determined by stable isotope dilution/mass spectrometry. J Clin Endocrinol Metab 71: 39-45, 1991.

20.   Goodman, LS, and Gilman AG. The Pharmacological Basis of Therapeutics (9th ed.). New York: McGraw-Hill, 1996.

21.   Guardabasso, V, DeNicolao G, Rocchetti M, and Rodbard D. Evaluation of pulse-detection algorithms by computer simulation of hormone secretion. Am J Physiol Endocrinol Metab 255: E775-E784, 1988[Abstract/Free Full Text].

22.   Hellman, L, Nakada F, Curti J, Weitzman ED, Kream J, Roffwarg H, Ellman S, Fukushima DK, and Gallagher TF. Cortisol is secreted episodically by normal man. J Clin Endocrinol Metab 30: 411-422, 1970[ISI][Medline].

23.   Johnson, NL, and Kotz S. Continuous univariate distributions-I. New York: Wiley, 1970.

24.   Jusko, WJ, Slaunwhite WR, Jr, and Aceto T, Jr. Partial pharmacodynamic model for the circadian-episodic secretion of cortisol in man. J Clin Endocrinol Metab 40: 278-289, 1975[Abstract].

25.   Keenan, DM. Implementation of a stochastic model of hormonal secretion. In: Quantitative Neuroendocriology, Methods in Neuroscience, edited by Johnson ML, and Veldhuis JD.. San Diego, CA: Academic, 1995, vol. 28, p. 311-323.

26.   Keenan, DM, and Veldhuis JD. A biomathematical model of time-delayed feedback in the human male hypothalamic-pituitary-Leydig cell axis. Am J Physiol Endocrinol Metab 275: E157-E176, 1998[Abstract/Free Full Text].

27.   Komaki, F. State-space modelling of time series sampled from continuous processes with pulses. Biometrika 80: 417-429, 1993[ISI].

28.   Leproult, R, Copinschi G, Buxtion O, and Van Cauter E. Sleep loss: sleep loss results in an elevation of cortisol levels the next evening. Sleep 20: 865-870, 1997[ISI][Medline].

29.   Linkowski, P, Mendlewicz J, Kerkhofs M, Leclercq R, Goldstein J, Brasseur M, Copinschi G, and Van Cauter E. 24-Hour profiles of adrenocorticotropin, cortisol, and growth hormone in major depressive illness: effect of antidepressant treatment. J Clin Endocrinol Metab 65: 141-152, 1987[Abstract].

30.   Linkowski, P, Mendlewicz J, Leclercq R, Brasseur M, Hubain P, Goldstein J, Copinschi G, and Van Cauter E. The 24-hour profile of adrenocorticotropin and cortisol in major depressive illness. J Clin Endocrinol Metab 61: 429-438, 1985[Abstract].

31.   Liu, JH, Kazer RR, and Rasmussen DD. Characterization of the twenty-four hour secretion patterns of adrenocorticotropin and cortisol in normal women and patients with Cushing's disease. J Clin Endocrinol Metab 64: 1027-1035, 1987[Abstract].

32.   Meehan, PM. A Bayesian Analysis of Diurnal Cortisol Patterns (PhD thesis). Cambridge, MA: Harvard University, 1993.

33.   Merriam, GR, Ma N, Liu L, Wachter KW, and Libre E. Methods for assessing the linkage between pulsatile hormone profiles. Acta Paediatr Scand Suppl 349: 167-172, 1989[Medline].

34.   Merriam, GR, and Wachter KW. Algorithms for the study of episodic hormone secretion. Am J Physiol Endocrinol Metab 243: E310-E318, 1982[Abstract/Free Full Text].

35.   Mortola, JF, Liu JH, Gillin JC, Rasmussen DD, and Yen SSC Pulsatile rhythms of adrenocorticotropin (ACTH) and cortisol in women with endogenous depression: evidence for increased ACTH pulse frequency. J Clin Endocrinol Metab 65: 962-968, 1987[Abstract].

36.  Munson PJ and Rodbard D. Pulse detection in hormone data: simplified, efficient algorithm. Proc Stat Comp Sec Am Stat Assoc Washington DC 1989, p. 295-300.

37.   O'Sullivan, F, and O'Sullivan J. Deconvolution of episodic hormone data: an analysis of the role of season on the onset of puberty in cows. Biometrics 44: 339-353, 1988[ISI][Medline].

38.   Papoulis, A. Probability, Random Variables, and Stochastic Processes (3rd ed.). Boston, MA: McGraw-Hill, 1991.

39.   Pincus, SM, and Goldberger AL. Physiological time-series analysis: what does regularity quantify? Am J Physiol Heart Circ Physiol 266: H1643-H1656, 1994[Abstract/Free Full Text].

40.   Reichlin, S. Neuroendocrinology. In: Williams Textbook of Endocrinology (9th ed.), edited by Wilson JD, Foster DW, Kronenberg HM, and Larsen PR.. Philadelphia, PA: Saunders, 1998, chapt. 8, p. 165-248.

41.   Rivier, C, Rivier J, and Vale W. Inhibition of adrenocorticotrophic hormone secretion in the rat by immunoneutralization of corticotropin-releasing factor (Abstract). Science 218: 377, 1982[ISI][Medline].

42.   Rogacz, S, Duffy JF, Ronda JM, and Czeisler CA. Endogeneous circadian temperature amplitude is reduced and phase is delayed at the luteal phase of the menstrual cycle in ovulating women (Abstract). In: Proc 1st Meet Soc Res Biol Rhythms, 1988.

43.   Santen, RJ, and Bardin CW. Episodic luteinizing hormone secretion in man. Pulse analysis, clinical interpretation, physiologic mechanisms. J Clin Invest 52: 2617-2628, 1973[ISI][Medline].

44.   Sherman, BC, Wysham C, and Pfohl B. Age-related changes in the circadian rhythm of plasma cortisol in man. J Clin Endocrinol Metab 61: 439-443, 1985[Abstract].

45.   Snyder, DL, and Miller MI. Random Point Processes in Time and Space. New York: Springer-Verlag, 1991.

46.   Spiegel, K, Leproult R, and Van Cauter E. Impact of sleep debt on metabolic and endocrine function. Lancet 354: 1435-1439, 1999[ISI][Medline].

47.   Straume, M, Johnson ML, and Veldhuis JD. Realistic emulation of highly irregular temporal patterns of hormone release: a computer-based pulse simulator. In: Quantitative Neuroendocrinology, Methods in Neurosciences, edited by Johnson ML, and Veldhuis JD.. San Diego, CA: Academic, 1995, vol. 28, p. 220-243.

48.   Vagunucci, AH. Analysis of circadian periodicity of plasma cortisol in normal man in Cushing's syndrome. Am J Physiol Regulatory Integrative Comp Physiol 236: R268-R281, 1979[ISI][Medline].

49.   Vale, W, Spiess J, Rivier C, and Rivier J. Characterization of a 41-residue ovine hypothalamic peptide that stimulates secretion of corticotropin and beta -endorphin (Abstract). Science 213: 1394, 1981[ISI][Medline].

50.   Van Cauter, E. Method for characterization of 24-h temporal variation of blood components. Am J Physiol Endocrinol Metab Gastrointest Physiol 237: E255-E264, 1979[Free Full Text].

51.   Van Cauter, E. Estimating false-positive and false-negative errors in analysis of hormonal pulsatility. Am J Physiol Endocrinol Metab 254: E1-E8, 1988[Abstract/Free Full Text].

52.   Veldhuis, JD, Iranmanesh A, Johnson ML, and Lizarralde G. Amplitude, but not frequency, modulation of adrenocorticotropin secretory bursts given rise to the nyctohemeral rhythm of the corticotrophic axis in man. J Clin Endocrinol Metab 71: 452-463, 1990[Abstract].

53.   Veldhuis, JD, Iranmanesh A, Lizarralde G, and Johnson ML. Amplitude modulation of a burstlike mode of cortisol secretion subserves the circadian glucocorticoid rhythm. Am J Physiol Endocrinol Metab 257: E6-E14, 1989[Abstract/Free Full Text].

54.   Veldhuis, JD, and Johnson ML. Cluster analysis: a simple, versatile, and robust algorithm for endocrine pulse detection. Am J Physiol Endocrinol Metab 250: E486-E493, 1986[Abstract/Free Full Text].

55.   Veldhuis, JD, and Johnson ML. A novel general biophysical model for simulating episodic endocrine gland signaling. Am J Physiol Endocrinol Metab 255: E749-E759, 1988[Abstract/Free Full Text].

56.   Weitzman, ED, Fukushima D, Nogeire C, Roffwarg H, Gallagher TF, and Hellman L. Twenty-four hour pattern of the episodic secretion of cortisol in normal subjects. J Clin Endocrinol Metab 33: 14-22, 1971[ISI][Medline].

57.   Williams, GH, and Dluhy RG. Diseases of the adrenal cortex. In: Harrison's Principles of Internal Medicine (14 ed.), edited by Fauci AS, Braunwald E, Isselbacher KJ, Wilson JD, Martin JB, Kasper DL, Hauser SL, and Longo DL.. New York: McGraw-Hill, 1998, chapt. 332, p. 2035-2057.


Am J Physiol Endocrinol Metab 280(3):E450-E461
0193-1849/01 $5.00 Copyright © 2001 the American Physiological Society