INTRODUCTION

Top Abstract Introduction Methods Discussion References

IN THE PRESENT PAPER , we implement a model for reconstructing from observed plasma hormone concentration time series the basal and pulsatile secretion rates of a protein hormone, such as luteinizing hormone (LH). Statistical methods for estimating the key parameters of the model are presented and applied illustratively to LH data for both the adult male and pre- and postmenopausal female. Then, based on the estimated parameters, we estimate (reconstruct) the unobserved secretion rate time series.

Consider the "true concentration" in vivo of a given single hormone. Let t₀ (0) represent the beginning of the observation period. Let [X(t ), t  t₀] be the hormone concentration evolving over time and [Z(t ), t  t₀] the hormone secretion rate over time. The rate of change in the concentration [X(t )] is described by the differential equation

(1)

with X(t₀ ) being some specified initial condition and  the unknown rate of elimination specific to the particular hormone of interest. All modeling of hormone dynamics starts (more or less) with this equation of conservation of mass, which merely states that the rate of change in concentration, by definition, is the difference in the rates of removal and production of the hormone.

We have chosen to model and analyze the hormone concentration [X(t )] and the rate of secretion [Z(t )] in continuous time. If one did not model with this perspective, it would be difficult to compare usual experimental data observed under different sampling schemes. We assume that, for a given subject, the observations of blood hormone concentrations are made at times t_k, k = 1, ... , n during the period [L₁, L₂], with t₀ = L₁, t_n = L₂. What is then observed is Y_k

(2)

as the concentration of the hormone plus sample measurement error (due, e.g., to assaying) at time t_k; the sampling interval t = t_k  t_k  1 is assumed to be constant, but need not be.

We postulate that an accurate description of hormone concentrations requires a stochastic formulation, which incorporates the various forms of uncertainty that occur within the different time and space scales. Here we allow for three such basic forms of variation.

First, at the cellular/glandular scale, there are variations in the instantaneous rates of synthesis of a given hormone within and among cells. For a hormone secreted in pulses, we distinguish between the instantaneous rate of synthesis (or intracellular production) and the instantaneous rate of secretion (which for a protein hormone is based on an accumulation and subsequent release of previously synthesized hormone-containing granules). Conceptually, the instantaneous rate of synthesis represents an "idealized" (or expected) rate, whereas the "realized" rate for a population of molecules and cells should be a stable random variation about this expected rate. This dispersion of variability (noise) results in the inclusion of the A_j terms in Eq. 4. Such random variations could reflect nonuniform within-cell and between-cell metabolic milieus (e.g., ATP energy stores), biochemical signals (e.g., first and second messengers), and cell biological functions (e.g., cytoskeletal apparatus state of polymerization, phosphorylation state of secretory-related proteins, etc.) (2, 5, 8, 9, 15, 18).

Second, at the level of glandular secretion of an array of hormone molecules into the circulatory system, there is nonuniform release topographically among the cells and subsequent mixing of the population of hormone molecules within the bloodstream, and this microscopic biological variability (noise) should impact the resulting concentration level. This noise is represented by the Brownian motion term _wdW(t ) in Eq. 6 (1, 11, 16, 18).

Third, at the level of the sample removal, processing, and assay from the human subject, there are additional contributions to experimental uncertainty (e.g., within-assay measurement errors). This is represented by the e_i terms in Eq. 8. Experimentally, the error in the automated measurements of protein hormones is estimated to have a typical coefficient of variation of 3-6% (5, 6, 16).

In Ref. 8, the pulsatile secretion (not concentration) of one hormone was modeled and applied to pituitary LH secretion data from a horse. This paper builds on that initial model in three ways. First, it extends the biomathematical construct and analysis to that of hormone concentration by including the modeling of hormone elimination; second, and as importantly, we here allow for greater flexibility in the modeling of the variable (LH) pulse masses via random effects, without introducing an inordinate number of parameters; and, third, we implement a maximum likelihood estimation (MLE) strategy with the calculation of appropriate secretory-parameter statistical confidence intervals. In Refs. 7 and 8 the pulse masses were assumed to be a linear function of the preceding interpulse lengths: the longer the interpulse interval, the greater the accumulation of LH mass. The former model is reasonable and appears to fit certain data well. However, such a model is rather rigid if a small interpulse interval is followed by a large mass (or vice versa), as we illustrate here. Failure to allow for more flexibility in possible burst mass values could produce spurious estimates of the true secretory variability (discussed below).

	METHODS

Top Abstract Introduction Methods Discussion References

Modeling Hormone Elimination and Secretion

Elimination. The elimination rate constant  of a biological molecule from a particular (single compartment) sampling space is related to its half-life (t_1/2 ) by: exp (t_1/2 ) = 1/2. For the components of the human male and female reproductive axes, we note that direct quantification of biexponential LH disappearance rates in hypopituitary men injected with highly purified human pituitary LH yielded a mean rapid (first) component LH half-life of ~18 min, a slow (second) component half-life of ~90 min, and an average fractional contribution of the former to total decay amplitude of 0.63 (20). Thus a monoexponential approximation of such kinetics would suggest a half-life within this broad range and its experimental uncertainties. Although the rapid and slower algebraic phases of hormone elimination do not necessarily correspond to definable anatomic compartments, the more rapid phase of hormone disappearance is thought to reflect largely hormone distribution within the vascular space after abrupt secretion or infusion, whereas the slower component may result from irreversible metabolic removal of the hormone from blood. Indeed, rate constants for the latter correlate with the sialic acid content of human LH infused into hypophysectomized rats, consistent with a view of a sialoreceptor-mediated uptake and removal of LH by metabolically relevant tissues (2). In accordance with this concept, primarily a distributional phase half-life may be evident during spontaneous LH secretion pulses, whereas irreversible metabolic removal would proceed more slowly. We discuss below that with rapidly recurring LH secretory pulses, the slower putative removal process may be less evident in the plasma LH concentration profile and could mimic a low rate of basal (nonpulsatile) LH secretion between peaks. This point will be illustrated further in a paradigm of infused LH pulses.

Pulsatile secretion. We have viewed the secretion model as arising in two stages. The first stage concerns the mechanism that governs the time occurrences of the bursts (or pulses); the second stage concerns the resulting shapes and masses of the bursts (at the various pulse times).

(3)

j  1

(4)

	MODELING HORMONE CONCENTRATION PROFILES

Given the foregoing physiological motivation, we thus propose as a biomathematical formulation of the rate of secretion Z( · ) and the concentration level X( · ) of pulsatile hormone secretion superimposed on a finite basal (time invariant) hormone secretion rate, ₀, the following

(5)

(6)

where A_j terms are IID n(0, ²_A).

To summarize, we formulated a construct to represent the intermittently observed hormone concentrations, starting at the level of the (unobserved) continuous rate of pulsatile secretion and allowing for random variation not only in the pulse times, but also via three other sources of stochastic variation (as reviewed in the introduction): 1) anticipated biological variation in the hormone mass accumulated between pulse times, which is reflected by the inclusion of the A_j terms in Eq. 4, where A_j terms are IID normal (0, ²_A); 2) nonuniform release into and subsequent mixing of the hormone molecules within the circulation, which is represented by the Brownian motion term _wdW(t ) in Eq. 6; 3) variability due to sample removal, processing, and assay (e.g., within-assay measurement errors and other technical variability), which is represented by the e_i terms IID normal (0, ²) in Eq. 8.

	RECONSTRUCTING THE PULSATILE HORMONE SECRETION RATE

Consider a discrete-time sampling of X( · )

(7)

What we actually observe is such a discrete-time sampling of X( · ) plus measurement error (e.g., due to assaying)

(8)

For simplicity, we assume that the sampling rate t is 1; our asymptotics concern n   rather than t  0, so this a natural assumption. For simplicity of notation, we will let represent what was 1  t above. Let Y = (Y₁, Y₂, ... , Y_n )' be the observed concentrations. We will condition on the (unobserved) Y₀ = y₀, pulse time T₀, and pulse mass M⁰; in practice, we estimate these three ( y₀, T₀, M⁰ ) from the data. Define where and where Above,  ~ N(0, ² I) represents the measurement error, and w ~ N(0, ²_w I) represents the increments of the Brownian motion. We can rewrite the above as Since Y has a normal distributed with mean and covariance matrix Using the above notation, the log likehood function is Let  = (₀, , ₀, ₁, ₁, ₂, ₃)',  = (², ²_A,²_w)', and  = (', ')'. Let g() =  µ and Using the above notation, one can express Y as Since A, w, and  are independent and normally distributed, we have a nonlinear random-effect model with the design matrices U₁, U₂ depending on .

In Ref. 24, the asymptotics for the MLE of   are established, with the MLE denoted by _n; it was precisely the present hormonal estimation problem that motivated the derivation of the asymptotics. The results were that

Also, we are interested in the actual values of random effects A, which embody biological variations in pulse mass (Eq. 4). Because the random effects A_j are not observable, we will use E(A_j|Y₁, Y₂, ... , Y_n); in Ref. 24 a precise formula for calculating these predicted values is presented. Also, to calculate the asymptotic standard deviations of the MLEs, one can calculate the information matrix I_n(_n ) on the basis of our model. The explicit expressions for the information matrix are given in Ref. 24.

On the basis of the parameter estimate _n and the predicted random effects we can reconstruct (or estimate) the pulsatile secretion rate, as follows

To calculate the total daily LH secretion, we integrate the reconstructed LH secretion rate, , from 0 to 1,440 min; because the normalized secretion rate ( · ) integrates to one, we can use the following very accurate approximation The above is an approximation only because of "edge effects" at the beginning and end of the day. We then obtain the daily pulsatile secretion by subtracting the daily basal secretion, ₀ × 1,440 min, from the total daily secretion.

Also, consider the (particular) likelihood equation with respect to ₀ (evaluated at the MLE _n )

(9)

That is, the sum of the predicted random effects (evaluated at _n ) is zero; this is analogous to the fitted residuals in linear least-squares regression, except that here the situation is much more complex.

We can use this approximation to calculate desired standard errors for daily LH basal and daily pulsatile LH secretion

Also, to obtain the standard error for the half-life (t_1/2 ), from the standard error for the elimination rate [ = log(2)/t_1/2 )], we use the -method (a first-order approximation).

	ILLUSTRATIVE APPLICATIONS

Adult Male and Female 24-h Serum LH Concentration Profiles

View larger version (45K): [in this window] [in a new window]   Fig. 1.   Illustrative 24-h serum luteinizing hormone (LH) concentration profiles from 2 young men, 3 premenopausal women, and 1 postmenopausal woman. Data shown by solid lines in each top paired subpanel are immunoradiometrically measured serum LH concentrations in blood collected at 10-min intervals, with the fitted (model reconstructed or predicted) curves depicted by dashed lines (fitted with random effects) and dotted lines (fitted without random effects). Bottom paired subpanels define the stochastic differential equation (SDE) model-calculated LH secretory rates, which exhibit a range of admixtures of basal and pulsatile LH secretion. Note the variable LH secretory pulse shapes, amplitudes, frequencies, and mass (integrals of the secretory bursts) and basal secretion rates illustrated in the 3 groups by gender and menopausal status. EF, early follicular; LF, late follicular; and ML, midluteal phase of the normal menstrual cycle in premenopausal women (data reanalyzed from Refs. 2, 15, 16, 22). Note different y-axes scales to accommodate the wide data range.

Table 1 gives the LH secretion statistics for the three illustrative LH data groupings: young males, premenopausal females, and postmenopausal females. SDs are given in parentheses. Table 2 shows the individual parameter estimates in two of the subjects. We also estimate parameter SD values for each of the 10 key parameters in the model.

Infused LH Pulse Profiles

Five models of combined basal and pulsatile LH release (infusion) were evaluated, and the known mass of LH injected (IU/volunteer) per pulse was regressed against the calculated mass of LH "secreted" per burst (IU/l of distribution volume) (Fig. 2A ). The inverse of the slope of this line approximates the LH distribution volume, which was measured earlier (20). In the first two models, the LH half-life was computationally estimated as a single exponential disappearance function for each LH dose infused, with or without constrained low basal rates of (residual) endogenous LH secretion. The latter reflected incomplete suppression by leuprolide and was calculated from the first measured (preinjection) basal serum LH concentration in each series. Compared with an expected LH distribution volume of 3.6-6.4 liters (75-86 kg subjects, with projected 4.5-8% of body weight as LH distribution volume), these two models overpredicted this value at 9.81 (variable basal) and 10.0 liters (constrained basal). Two other models, one with and the other without a constrained low basal rate of endogenous LH secretion, assumed known a priori measured biexponential LH disappearance kinetics, namely, a rapid and slower mean LH half-life of 18 and 90 min, respectively, with 37% of the total decay amplitude attributable to the slower component (20). Both models predicted mean LH distribution volumes similar to expectation, namely 3.96 (variable basal) and 3.60 liters (constrained basal). A fifth model also allowed the foregoing two-component LH kinetics, but assumed zero residual (basal) LH secretion and yielded an estimated LH distribution volume of 3.76 liters. Examples of the model-based fits to the serum LH concentration profiles, as well as the reconstructed LH secretion rates, after the 7.5 (low) and 50 IU (higher) dose LH infusions for three models are shown in Fig. 2A. The three models illustrated are zero basal, constrained basal, and variable basal (endogenous) LH secretion.

View larger version (19K): [in this window] [in a new window]   View larger version (38K): [in this window] [in a new window]   Fig. 2.   A: linear regression plots of relationships between injected dose of human recombinant LH (IU of First International Reference Preparation) and calculated mass of LH infused (IU/l distribution volume) for 5 proposed models of admixed basal and pulsatile LH secretion/infusion. Men were pretreated 3-4 wk earlier with leuprolide acetate (3.75 mg im) to downregulate (endogenous) LH secretion and then infused intravenously on separate days in randomly assigned order with pulses (over 0.5 to 8 min) of human LH at fixed doses of 7.5, 15, 30, 50, or 75 IU every 2 h for 4-8 consecutive injections. Blood was sampled once before and then every 10 min during the LH infusions for a total of 8-16 h. Resultant serum LH concentration profiles (assayed by 2-site immunoradiometric assay) were analyzed allowing a variable (fitted) 1-exponential decay with a constrained or freely varying basal (nonpulsatile) endogenous LH secretion rate (bottom 2 lines) and a known 2-exponential LH decay model (see Ref. 20 and METHODS) with a constrained versus freely varying versus zero basal LH secretion rate (top 3 lines). Inverse slopes of the regressions yield the mean (LH) distribution volume. Regressions and corresponding distribution volumes are constrained basal, 1 half-life: y = 0.4772 + 0.0997x and (1/0.0997) = 10 liters (A); freely varying basal, 1 half-life: y = 0.2707 + 0.1020x and (1/0.1020) = 9.8 liters (B); constrained basal, 2 constrained half-lives: y = 2.64 + 0.280x and (1/0.280) = 3.6 liters (C); freely varying basal, 2 constrained half-lives: y = 1.894 + 0.2525x and (1/0.2525) = 3.96 liters (D); and zero basal, 2 constrained half-lives: y = 0.302 + 0.266x, and (1/0.266) = 3.76 liters (E). B: illustrative fits of 2 of the intravenous injected human recombinant LH dose profiles, namely 7.5 and 50 IU LH infused every 2 h, as described in A. Each LH dose profile within a paired subpanel set shows fits for 2 half-lives (fixed) (dashed lines) and/or freely varying half-lives (monocomponent) (dotted lines). We further illustrate 3 (paired subpanel) plots for each dose schedule, namely, constrained basal, freely varying basal, and zero basal (LH) secretory rates.

	DISCUSSION

Top Abstract Introduction Methods Discussion References

Here, we develop, implement, and illustrate a stochastic differential equation (SDE) biomathematical formulation of combined basal and pulsatile LH secretion in concert with a subject-specific LH half-life and conditional pulse times to quantitatively analyze (24 h) serum LH concentration time series, e.g., as obtained earlier in healthy young men and premenopausal as well as postmenopausal women (5, 14-16). This new construct of hormone secretion and removal allows for variably shaped [e.g., skewed or asymmetric (1, 11)] LH secretion pulses superimposed upon a time-invariant basal (LH) secretion rate. In addition, secreted LH molecules are admixed in the bloodstream and subjected to a monoexponential (or higher order) elimination function (20, 21). The smaller number of essential parameters defining overall LH secretion and removal in this model [namely 10 parameters per 144-sample data set according to the present notion, versus 20-30 parameters by model-specific deconvolution analysis (18)] allows us to explore the otherwise difficult issue of joint estimation of basal and pulsatile hormone secretion rates assessed concurrently with half-life (19). Moreover, the statistical basis for the present MLE of basal and pulsatile LH secretion and LH half-life permits the construction of statistical confidence intervals for each of these parameters, as well as a reconstruction of the random-effects term defining stochastic variability in anticipated LH secretory burst mass.

The high correlations among basal and pulsatile rates of LH secretion and concurrent half-life of elimination of LH in earlier highly parameterized convolution models of neurohormone release and removal pose a formidable challenge to reliable estimation of these interdependent parameters by parametric approaches (19). The current implementation of a more parsimonious model of LH secretion and removal, including relevant stochastic contributions at several levels and a restricted resultant parameter set, begins to address some of these constraints, and might thus ultimately allow resolution not only of pulse number and mass, basal hormone release, and hormone half-life, but also of asymmetric LH secretory-pulse shape. Our use of the generalized gamma function (with 3 parameters) to mathematically depict a potentially asymmetric LH secretory burst shape seems to be suitable for evaluating the hormone secretion rate contour over time within a pulse (see preliminary analysis in Table 2 and direct sampling data in Refs. 1, 11). Such estimates should be technically accomplishable in peripheral blood with greater accuracy under conditions of sufficiently rapid blood sampling and adequately specific and reproducible hormone assays. For example, a recent clinical study sampled blood every 2.5 min in young and older men throughout the night, allowing a high density of serum LH measurements over time (13). Our formalization of pulse shape via a simple three-parameter gamma function, with a term to define the steepness of upstroke, another to depict the rapidity of declining secretion after its maximum, and a third to impart peak sharpness quantification should be useful in eventually appraising secretory event variations in different pathophysiological states. Without highly specific and precise assay methods, however, neither this nor other models would likely allow accurate discrimination of (LH) secretory pulse shape. Moreover, when pulse shape is reconstructed, we recommend comparisons of analytic results with direct secretory gland sampling, e.g., as carried out in the intact horse (1), ovariectomized sheep (11), or human (17).

By applying the present SDE nonlinear random-effects model to physiological serum LH concentration time series in men and young and older women, we could illustrate possible contrasts among healthy individuals with respect to LH secretory burst frequency and mass and basal LH secretion rate but not apparent (distributional) LH half-life. In particular, the postmenopausal woman exhibited a higher rate of calculated basal LH secretion than the young men, without any major disparity in (distributional) LH half-life. In addition, LH secretory burst mass was considerably (4-fold) higher in the older woman than that in young men. In contrast, in the young woman in the midluteal phase of the menstrual cycle, LH pulse frequency was low with an interpulse interval of ~2.5-3 h. In the luteal phase, the apparent basal LH secretion rate was intermediate and contributed ~50% of total daily LH secretion. Thus the postmenopausal woman and young men may represent extremes in basal LH secretion rates and LH secretory burst mass. The healthy young woman in the early follicular phase of the menstrual cycle seems to represent near-maximal pulsatile LH secretion (74 ± 6%) and the lowest basal LH secretion rate. During the late follicular phase, an accelerated LH pulse frequency, increased mass per burst, and augmented basal LH release rate are suggested here. We emphasize that considerable further analyses in a larger group of men and women will be important to definitively test the generality of these illustrative inferences. The basis for our estimates must be distinguished from that of earlier deconvolution estimates of LH secretion (15) in that here we implement an SDE model with random effects, jointly estimate combined pulsatile and nonpulsatile (basal) LH release, condition our analysis of secretion rates on independently estimated pulse times, and apply maximum likelihood statistical estimation (MLE).

Our estimates of LH half-life during pulsatile LH release under physiological conditions in the human male and female are consistent with independent calculations of LH distribution rates (rapid-phase elimination) after the bolus injection of highly purified human pituitary LH in hypopituitary men (20). Namely, such earlier experiments showed an initial distributional phase LH half-life of disappearance from plasma of ~18 min, a delayed (slow, second) component half-life of 90 min, and a fractional contribution of the rapid component to the total amplitude of LH decay of 0.63 (20). Here we find a range of apparent (rapid) half-lives of endogenous pulsatile LH decay toward basal of 7.4-27 min, thus suggesting that in vivo human pituitary LH secretion occurs with sufficient frequency that steady state is rarely attained in plasma and that the rapid phase of LH distribution tends to predominate. Alternatively, we point out that these data either argue against a combined model of basal (nonpulsatile) and pulsatile LH secretion or suggest a biexponential LH decay structure in these physiological states. Data by way of validation (Fig. 2A ) suggest the latter (biexponential) interpretation.

Earlier deconvolution analyses assuming purely pulsatile LH secretion (zero basal) predict an LH half-life range of 60-130 minutes in human subjects (5, 15, 16, 21). Such values approximate the directly measured second (slow) component of LH removal after bolus intravenous injection of highly purified human pituitary-derived LH in LH-deficient men (90 min) (20) or the half-lives of the metabolic removal of LH during (12, 20, 23) or after (10) steady-state LH infusions in the human (80-130 min). Accordingly, the shorter half-lives of LH disappearance estimated in our combined pulsatile-basal secretion model with a nonlinear random-effects SDE analysis reflect either primarily distributional (rapid phase) kinetics of LH and/or suggest an overestimation of the basal LH secretory rate due to (e.g., in postmenopausal women) rapidly successive LH pulse times with incomplete LH removal before the onset of the next secretory event and a coarse sampling rate (relative to pulsing rate). These issues are overcome largely via a biexponential decay structure (see below).

We explored the foregoing kinetic considerations further in experiments using bolus-injected human recombinant LH in five young men pretreated with the GnRH agonist leuprolide to produce a reversible LH deficiency state (see Fig. 2). In this clinical paradigm of experimentally reduced rates of basal LH release and known fixed exogenous LH injection doses given intravenously at 2-h intervals, a single-exponential formulation of LH decay would yield overestimates of LH distribution volume (bottom 2 curves in Fig. 2A ), underestimates of LH half-life, and overestimates of basal LH secretion rate. Accordingly, we caution that misapplication of a combined basal and pulsatile hormone secretion model with monoexponential kinetics to a known (virtually) purely pulsatile hormone release context may predict unduly short hormone half-lives and excessive basal and total hormone secretion rates. A similar inference is made using bolus testosterone injections after pharmacological testosterone depletion by ketoconazole administration in men (not shown).

Using the current SDE construct, we have not attempted to fit the human LH data to a variable two-component half-life model of elimination, which would increase the parameter set. In principle, adding a second (slower) component of LH removal, possibly reflecting irreversible tissue uptake and degradation of LH (see introduction), would reduce the estimates of basal secretion presented here, while tending to increase the apparent mass of hormone secreted within bursts. Indeed, our experiments using varying doses of infused LH corroborate this prediction (see top 3 curves in Fig. 2A ). Thus further evaluation of biexponential kinetic models of hormone disappearance will be quite relevant and instructive. In addition, we would note that independent prior knowledge of anticipated or known concurrent basal rates of (nonpulsatile) hormone release, whether zero or otherwise, would be helpful when appraising complex pulsatile hormone secretion profiles quantitatively.

Our formalization of LH secretory pulse mass includes a residual mass of LH accumulated since the last GnRH-stimulated discharge or pulse and a constant that relates the rate of pulse-mass accumulation to the prior interpulse interval length (respectively, ₀, ₁ ). We have attempted here (Table 2) preliminarily to estimate these two conceptually distinct contributions to the mass of any given burst of secreted LH. When estimating secretion from high-intensity sampling data with infrequent pulse events and/or from larger