## Abstract

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

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

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* - A
_{j} - Total amount of cortisol in a secretory event initiated at time
*u*(μg/dl)_{j} - A
- Vector of secretory event amplitudes [
*A*_{1},*A*_{2},…,*A*_{N}_{(T)}]^{T} *c*_{r}- Coefficient of the
*r*th cosine harmonic of the mean circadian amplitude function *d*_{r}- Coefficient of the
*r*th sine harmonic of the mean circadian amplitude function *H*_{1}(*t*)- Cortisol concentration in the adrenal glands at time
*t* *H*_{2}(*t*)- Cortisol concentration in the plasma space at time
*t* *N*(*t*)- Number of secretory events occurring in an interval (0,
*t*] *u*_{j}- Time of the
*j*th secretory event *u*- Vector of secretory event times [u
_{1},u_{2},…,u_{N}_{(t)}] *w*_{k}- Waiting until the
*k*th secretory event given the time of the*k*− 1st secretory event *y*_{t}- Cortisol concentration measured in the plasma at time
*t*(μg/dl) - β
_{C} - Rate constant for clearance of cortisol from the plasma (min
^{−1}) - β
_{I} - Rate constant for infusion of cortisol from the adrenal glands into the plasma (min
^{−1}) - γ
_{A} - Circadian amplitude function coefficient of variation
- λ
_{1} - Location parameter of the gamma waiting time probability density
- λ
_{2} - Scale parameter of the gamma waiting time probability density
- μ
_{w} - Mean waiting time between secretory events
- ς
^{2}_{μ(t)} - Variance of the circadian amplitude function at time
*t* - ς
^{2}_{w} - Variance of the intersecretory event times
- Σ
_{β} - Covariance matrix of the cortisol infusion and clearance rate constants
- Σ
_{ξ} - Covariance matrix of the mean circadian function parameters
- Σ
_{λ} - Covariance matrix of intersecretory event waiting time parameters
- θ
- Vector of model parameters: θ = [β
_{I},β_{c},λ_{1},λ_{2},*c*_{0},*c*_{1},*d*_{1},*c*_{2},*d*_{2}] - ξ
- Vector of the mean circadian function parameters: ξ = [
*c*_{0},*c*_{1},*d*_{1},*c*_{2},*d*_{2}]

## METHODS

#### Model formulation.

Given an observation interval (0,*T*], we define at time*t*, H_{1}(*t*), the concentration of cortisol in the adrenal gland, and H_{2}(*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* ∈ (0,*T*]. We let*u _{j}
*,

*j*= 1, … ,

*N*(

*T*) be the times in (0,

*T*] at which secretory episodes are initiated, and we let

*w*,

_{k}*k*= 1, … ,

*N*(

*T*) be the waiting times between secretory events. It follows from the definitions of

*u*and

_{j}*N*(

*t*) that Equation 1 Equation 2where d

*N*(

*u*) indicates a change in the counting process at

*u*. Heuristically speaking, d

*N*(

*u*) is 1 if there is a secretory event at

*u*, and 0 otherwise. We assume the value for

*w*is a renewal process defined as independent, identically distributed gamma random variables whose probability density has location and shape parameters λ

_{k}_{1}and λ

_{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

*w*are respectively Equation 3 Equation 4Let

_{k}*A*,

_{j}*j*= 1, …,

*N*(

*T*) be the amplitude or total amount of hormone contained in the secretory event initiated at time

*u*. We assume that the magnitude of the

_{j}*A*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 μ(

_{j}*t*) and variance function ς as Equation 5 Equation 6where γ

_{A}> 0 is a coefficient of variation. We assume that each

*A*is a Gaussian random variable with mean μ(

_{t}*t*) and variance ς .

Let β_{I} be the infusion constant governing the rate at which cortisol enters the plasma from the adrenal gland and β_{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
Equation 7
Equation 8The solution to this equation system on the interval (*t*
_{0},*t*] is
Equation 9where *H*(*t*) = [*H*
_{1}(*t*),*H*
_{2}(*t*)]′,*G*(*t*) = [*A _{t}
*,0]′ and
Equation 10where κ = β

_{I}/(β

_{I}− β

_{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*y*
_{t1}, … ,*y*
_{tn
} denote the concentration of cortisol at times 0 <*t*
_{1} < *t*
_{2} <...<*t _{n}
*

_{ − 1}<

*t*≤

_{n}*T*and assume the

*y*

_{ti }satisfy the equation Equation 11where ε

_{ti }are independent, zero mean Gaussian random variables with variance ς

^{2}

_{ε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*A _{t}
*d

*N*(

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

*A*, is a Gaussian random variable whose mean (

_{t}*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
Equation 12where *y*= [*y*
_{t1}, … , *y*
_{tn
}],*A* = [*A*
_{1}, … ,*A _{N}
*

_{(T)}],

*u*= [

*u*

_{1}, … ,

*u*

_{N}_{(T)}], and θ = [β

_{I},β

_{C},λ

_{1},λ

_{2},

*c*

_{0},

*c*

_{1},

*d*

_{1},

*c*

_{2},

*d*

_{2}];

*f*(

*y*‖

*A*,

*u*,θ) is the multivariate Gaussian probability density defined in

*Eq.11*;

*f*(

*A‖u*,θ) is the joint probability density of the secretory event amplitudes whose mean and variance follow implicitly from

*Eqs. 5 and 6*; and

*f*(

*u*‖θ) 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 θ, then simulation of plasma cortisol data from *Eq. 12
* with θ fixed is equivalent to studying variation in the hormone's plasma levels for an individual. Therefore, we consider the expanded model
Equation 13where *f*(θ) denotes the probability density that summarizes the between-individual variation in θ. We simulate our cortisol model by first drawing θ from *f*(θ) and then simulating *y*, *A*, and *u* given θ from*f*(*y*,*A*,*u*‖θ). 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*(θ). The within-individual variation is summarized by*f*(*y*,*A*,*u*‖θ). 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*(θ) is a multivariate Gaussian probability density, which has the form
Equation 14where λ = [λ_{1},λ_{2}], ξ = [*c*
_{0},*c*
_{1},*d*
_{1},*c*
_{2},*d*
_{2}], and β = [β_{I},β_{C}]. The probability densities *f*(λ), *f*(ξ), and*f*(β) 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*(θ).

#### Ultradian parameter probability density: f(λ).

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 ς
, 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 λ_{1} and λ_{2}can be computed as (23)
Equation 15
We derive the covariance matrix for λ_{1} and λ_{2} from the data presented by Veldhuis and colleagues, since they provided estimates of μ_{w} and ς
for six individual subjects (53). We used *Eq. 15
* to convert each subject's estimate of μ_{w} and ς
into estimates of λ_{1} and λ_{2}. Using the six pairs of λ_{1} and λ_{2} along with estimates of the means of λ_{1}and λ_{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
Equation 16Hence, we take *f*(λ) to be the bivariate Gaussian probability density whose mean vector and covariance matrix are given in *Eq. 16
*.

#### Kinetic parameter probability density: f(β).

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 53to 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 β_{C} using the delta method formula (38)
Equation 17
Equation 18where
and *s*
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 β_{C} using *Eq. 16
*. We combined the estimates of the mean of β_{C} and the variance of β_{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 β_{C} of 0.645 min^{−1} and 0.015 min^{−2}.

The analyses of Veldhuis et al. (53) provide information on β_{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 β_{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 β_{I}because it was based on only six subjects. Therefore, because the sample variance estimate is approximately distributed as a χ^{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 β_{I}. Because we identified no reports that jointly estimated β_{C} and β_{I}, we assumed the covariance between β_{C}and β_{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
Equation 19

#### Circadian parameter probability density: f(ξ).

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
Equation 20
where the lower triangle is omitted because of symmetry. Therefore, we took *f*(ξ) to be the Gaussian probability density whose mean vector and covariance matrix are defined in *Eq. 20
*. The value of γ_{A}, the coefficient of variation, could not be determined from previous studies. Therefore, we simulate our model for the following three values of γ_{A}: 0.1, 0.3, and 0.5.

#### Joint probability density of the model parameters: f(θ).

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
Equation 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 < *t*
_{1} <*t*
_{2} <...< *t*
_{n−1} <*t _{n}
* ≤

*T*, the simulation algorithm for the cortisol model defined in

*Eqs. 1-9*proceeds through the following seven steps

All simulations were carried out using the Splus programming language.

## RESULTS

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. 2
*A*. 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.

*Step 5* of the algorithm integrates the stochastic differential model in *Eqs. 7
* and *
8
* to obtain the solution in *Eq. 9
* (Fig. 2
*C*). The impulse response function for this system, as shown in Fig. 2
*B*, is derived directly from *Eq. 9
*. That is, for a given secretory event, i.e., a *u _{k}
* and an

*A*, the infusion and clearance parameters govern the fraction of

_{k}*A*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

_{k}*A*enters the plasma. By 6 h, all of the cortisol in this secretory event has reached the plasma. In

_{k}*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. 2

*D*). 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 γ_{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 θ from *f*(θ), and for each value of γ_{A}, we simulated six 48-h cortisol series. Figures 3, 4, and 5 show the same value of θ. That is, for example, Figs. 3
*A*, 4
*A*, and 5
*A*were all simulated with the same random draw of θ. Figures 3-5differ only in their choice of γ_{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.

As expected, the variability in the plasma cortisol levels increases as γ_{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. 3
*A* the plasma levels in the second 24 h are higher. For γ_{A} = 0.1 (Fig. 3), there is less between-day variation in the plasma cortisol levels. The simulated data from each value of γ_{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 to5 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. 2
*B*). 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. 2
*B*. 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.

To appreciate better within-subject variation in plasma cortisol levels, we simulated an additional set of six cortisol series with θ fixed at θ_{0} and γ_{A} = 0.3 (Fig.7). That is, for all six draws, θ is set equal to θ_{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.

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.

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.8
*A*). 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

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.

## Acknowledgments

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.

- Copyright © 2001 the American Physiological Society