## Abstract

Luteinizing hormone (LH) administered in pharmacological amounts downregulates Leydig cell steroidogenesis. Whether reversible downregulation of physiological gonadotropin drive operates in vivo is unknown. Most of the analytical models of dose-response functions that have been constructed are biased by the assumption that no downregulation exists. The present study employs a new analytical platform to quantify potential (but not required) pulsatile cycles of LH-testosterone (T) dose-response stimulation, desensitization, and recovery (pulse-by-pulse hysteresis) in 26 healthy men sampled every 10 min for 24 h. A sensitivity-downregulation hysteresis construct predicted marked hysteresis with a median time delay to LH dose-response inflection within individual T pulses of 23 min and with median T pulse onset and recovery LH sensitivities of 1.1 and 0.10 slope unit, respectively (*P* < 0.001). A potency-downregulation model yielded median estimates of one-half maximally stimulatory LH concentrations (EC_{50} values) of 0.66 and 7.5 IU/l for onset and recovery, respectively (*P* < 0.001). An efficacy-downregulation formulation of hysteresis forecasts median LH efficacies of 20 and 8.3 ng·dl^{−1}·min^{−1} for onset and offset of T secretory burst, respectively (*P* = 0.002). Segmentation of the LH-T data by age suggested greater sensitivity, higher EC_{50} (increased LH potency), and markedly (2.7-fold) attenuated LH efficacy in older individuals. Each of the three hysteresis models yielded a marked (*P* < 0.005) reduction in estimated model residual error compared with no hysteresis. In summary, model-based analyses allowing for (but not requiring) reversible pituitary-gonadal effector-response downregulation are consistent with a hypothesis of recurrent, brief cycles of LH-dependent stimulation, desensitization, and recovery of pulsatile T secretion in vivo and an age-associated reduction of LH efficacy. Prospective studies would be required to prove this aging effect.

- testosterone
- gonadotropin
- hysteresis
- downregulation
- dose response

biological systems commonly exhibit hysteresis,^{1} defined broadly by dependence of a response on its previous history, thus implying that nonidentical initial (onset) and delayed (recovery) pathways operate sequentially (1). Classical instances arise in muscle, lung, and gallbladder contraction and in hormonal signaling (5, 8, 14, 36). Prominent endocrine examples include sex steroid, glucose, insulin, glucagon, and parathyroid hormone stimulus-response adaptations (7, 43, 53). Unfortunately, hysteretic dose-response fluctuations have usually been ignored in analytical models or imputed to random effects (58). A random-effects model was introduced recently for the pituitary-gonadal axis and justified on statistical grounds (6, 29); however, if de facto hysteretic processes operate over time, a purely random-effects dose-response model would be overly simplistic. Accordingly, more precise understanding of endocrine-signaling systems would ultimately require allowance for and quantification of possible hysteresis, if present.

The present work investigates application of a new analytical construct of allowable dose-response hysteretic adaptations under physiological conditions in vivo. Model constructs permitted analytical estimation of possible (but not obligatory) nonidentical initial and recovery dose-response parameters,^{2,3} i.e., potency, efficacy, or sensitivity, as recently suggested for corticotropin-cortisol regulation in humans (26). The analysis targets the gonadal axis. Studies of the gonadal axis have revealed unequivocal pharmacological downregulation of gonadotropin/human chorionic gonadotropin (hCG) drive of steroidogenesis (9, 19, 20, 34, 51, 52, 60, 61). The issues are whether allowable cycles of physiological downregulation can be modeled in vivo in unstressed uninfused individuals and whether hysteretic adaptations are influenced by physiological context, such as age.

## METHODS

#### Overview.

Three distinct hysteretic constructs (potency, sensitivity, and efficacy) were evaluated using archival data in 13 young (ages 18–30 yr) and 13 older (ages 60–78 yr) men. Healthy inclusion criteria, along with diet, activity, and weight status, are described elsewhere (32, 56). Body mass index ranges were similar: 20–26 and 22–28 kg/m^{2} for young and older groups, respectively. Upon informed consent and with Institutional Review Board approval, 26 healthy subjects underwent blood sampling every 10 min for 24 h in a clinical research unit for later measurements of serum luteinizing hormone [LH (effector)] and testosterone [T (response)] concentrations. These data were subjected previously to linear cross-correlation and conventional dose-response analyses. They are employed here in an entirely new analytical context.

#### Analytical construct.

The goal was to relate time-varying LH concentrations (input and effector) to time-varying T secretion rates (output and response) via a new hysteresis dose-response model in individual young and older men. This requires construction of a five-parameter nonlinear (logistic) dose-response function with allowance for inhibitory (negative) hysteresis after an estimable time lag. Figure 1 illustrates the concept, which was first outlined by Keenan et al. (26) for corticotropin's feedforward onto cortisol secretion. For the gonadal axis, LH concentration profiles (145 samples in each subject over 24 h) are reconstructed by a modified deconvolution model to yield fitted (reconvolution) concentration values. LH secretion is represented by a train of variable-amplitude pulses superimposed on basal (time-invariant) nonpulsatile secretion (39), whereas LH elimination proceeds via a biexponential function (28). The modified deconvolution procedure incorporates an integral equation model, rather than the usual difference-based model (see appendix b. *Estimation Algorithm. Stages I–III*). An advantage of the integral form is that elimination and secretion are viewed as occurring continuously between consecutive samples, which corresponds to biological reality. In addition, the modified methodology presented here includes *1*) permissible addback of initially deleted pulses if significant via the Akaike information criterion (AIC) and *2*) allowance for two possible LH secretory burst shapes (waveforms), one during the day and the other at night, as suggested recently for ACTH (27). T concentrations are deconvolved analogously to yield sample T secretion rates for subsequent use in construction of dose-response hysteresis models. Then, after alignment of pulsatile LH concentrations and T secretory responses, the protocol estimates five-parameter nonlinear LH → T feedforward (logistic-like) dose-response functions. The unique aspect is estimation of a classical four-parameter (basal, potency, sensitivity, and efficacy) dose-response model simultaneously with a fifth downregulation parameter and a hysteresis-inflection time point. The inflection point is the time delay before allowable downregulation of the LH-T dose-response function optimized over the collection of T pulses. Downregulation is implemented as an allowable decrease in potency, sensitivity, or efficacy of LH-T action after the estimable inflection point from the onset of the T secretory pulse.

Figure 1 schematizes the pulsatile downregulation hysteresis concept. In particular, LH concentrations (reconvolved) are related to T secretion rates (deconvolved) via five simultaneously estimable dose-response parameters and a time-delayed hysteresis inflection point. Random effects on T secretory-burst amplitude and the standard deviation (SD) of residual model error are calculated concurrently, as described previously (23, 24).

The existence of parameter shifts was tested by the signed-ranks test of hysteretic onset and recovery parameters. Justification for each model was evaluated by likelihood-ratio testing using a χ^{2} statistic: each hysteresis model was compared with a no-hysteresis construct on a subject-by-subject (*n* = 26) basis using twice the (negative) log-likelihood function differences at two degrees of freedom. The final step involved assessment of statistical contrasts by age and also by type of downregulation model using Wilcoxon's rank-sum (2-sample) test (young vs. older comparison) and two-way analysis of covariance (ANCOVA) of logarithmically transformed model error SD values (effect of age and model type), respectively. Multiple means were compared by post hoc Tukey's honestly significantly different test (45, 46).

## RESULTS

The modified integrative-deconvolution methodology (see appendixes a–c) has never been applied to LH concentration-time series. Hence, significant results are presented in Fig. 2. Statistical comparisons disclosed prominent effects of age (13 older vs. 13 young men). Confirmatory of single-waveform analyses (30), LH secretory bursts per 24 h were more frequent (*P* = 0.008) and LH secretory bursts were smaller (lower mass of LH released per pulse) in older than young men (*P* = 0.007). New findings were that the daytime LH secretory-burst mode (time delay from objective secretory-burst onset to maximum) was 45% shorter (*P* = 0.005), whereas basal (nonpulsatile) LH secretion was 58% higher (*P* = 0.035), in older than young men. The nighttime modes of LH secretory bursts also differed by age (*P* = 0.043): median 16 (range 3–18.4) min for young men and 12 (range 3–17.7) min for older men.

Figure 3 shows estimates of initial and recovery (hysteretic) dose-response parameter values for models of sensitivity, potency, and efficacy hystereses. *P* values reflect paired nonparametric comparisons of onset vs. recovery LH-T dose-response parameter estimates in the group of 26 subjects. All three models predicted strong downregulation (*P* < 0.005). In the sensitivity model (Fig. 3*A*), 25 of 26 men exhibited a hysteresis-like downregulating shift. Median onset and recovery sensitivities were 1.1 and 0.10 slope unit, defining an 11-fold shift (*P* < 0.001). In the potency model (Fig. 3*B*), LH-T potency shifts were inferable in all 26 men. Median onset and offset EC_{50} values for the 26 subjects were 0.66 and 7.5 IU/l, respectively (*P* < 0.001, paired comparison). Figure 3*C* presents estimates of LH-T efficacy (asymptotically maximal T secretion). The median onset and offset efficacy values were 20 and 8.3 ng·dl^{−1}·min^{−1}, respectively (*n* = 26 individuals, *P* = 0.002). Downregulation of efficacy was inferable in all subjects.

Table 1 shows ANCOVA results of model error (SD) comparisons, with the no-hysteresis model SD as the covariate and segmentation by age (2 factors) and model type (3 factors). Each hysteretic model reduced mean model residual error significantly (all *P* < 0.005). Among the three hysteresis models, residual model error was lower for the potency model than for the efficacy construct (*P* = 0.011). Main effects of age and model type were significant at *P* < 0.001, but there was no significant interaction (*P* = 0.32). Residual model error (SD) was lower by age only in the sensitivity model (*P* = 0.025). Compared with no hysteresis, model justification was confirmed by likelihood-ratio testing using a χ^{2} comparison of twice the difference of (negative) log-likelihood functions (Table 1).

Individual outcomes for the three hysteresis models in a young and older man are illustrated in Fig. 4. Sensitivity, efficacy, and potency model estimates were obtained in a young and an older man. Figure 4 shows the deconvolved T secretion estimate and hysteresis-predicted T secretion fit, maximum-likelihood estimates (MLE) of logistic dose-response parameters, and the LH concentration (input) signal, with and without the allowable LH-T alignment shift. Estimated initial (onset) mean dose-response curve and delayed (recovery) mean dose-response curve (after the hysteretic time delay from the T pulse onset) are shown. Each set of dose-response curves reflects multiple pulse-by-pulse dose-response hysteretic cycles, which differ only by allowable random effects in efficacy.

Median sensitivity shifts were 8.6-fold larger in the older (1.8 slope units, *n* = 13) than the young (0.21 slope units, *n* = 13) cohort (*P* < 0.001), indicating greater downregulation (Table 2). Figure 5 show that absolute values of onset (initial) and offset (recovery phase) exponential sensitivity estimates were higher in older than young individuals (*P* < 0.001 and *P* = 0.004, respectively). Efficacy shifts did not differ by age (*P* = 0.34). These data (median and range) are summarized by age and model in Table 3. Analysis of data from an earlier study wherein ketoconazole was used to block T synthesis partially (59) showed that hysteresis was abolished (*n* = 6 men, plots not shown).

Potency (see appendix) also is an exponential dose-response term that yields an EC_{50} (estimated LH concentration stimulating one-half maximal T secretion) when the potency value is divided by the sensitivity coefficient (23, 27). Because potency in the potency-hysteresis model and sensitivity in the sensitivity-hysteresis model are estimated before (onset) and after (recovery) LH-associated downregulation, two EC_{50} values can be calculated for each subject. There is only a single EC_{50} estimate in the efficacy-shift model. For the potency- and sensitivity-shift models in the combined cohorts (*n* = 26), onset and recovery EC_{50} values were 0.66 and 7.5 IU/l and 2.7 and 15 IU/l, respectively (both *P* < 0.001). Median EC_{50} shifts (IU/l), i.e., the difference between onset and recovery values, in the potency-hysteresis model were similar in young and older men (*P* = 0.42, age effect; Tables 2 and 3). In contradistinction, the downregulating shift in LH EC_{50} was 3.9-fold smaller in older than young men in the sensitivity-hysteresis model: 4.5 vs. 13 IU/l (*P* = 0.015). At the same time, individual onset and recovery EC_{50} values in the sensitivity model were reduced (apparent potencies were enhanced) in the older men (*P* = 0.001 for onset and *P* = 0.011 for recovery; Fig. 6). Assessment of efficacy by age stratum predicted 2.7-fold higher LH efficacy in the sensitivity-shift model in young than older men: 43 (range 20–247) vs. 16 (range 5.9–96) ng·dl^{−1}·min^{−1} (*P* = 0.004; Fig. 7).

Estimated time delays (min) from onset of a T secretory burst to downregulation (hysteretic-inflection delays) were similar in all three models and invariant of age (*P* ≥ 0.21, ANCOVA). The median (absolute range) across the 26 subjects and 3 models was 23 (10–43) min. Model error (SD) was lower in the older cohort.

## DISCUSSION

The present analyses demonstrate that hysteresis models of allowable reversible quenching of pulsatile gonadotropin (LH)-T stimulation can be constructed and applied to physiological data obtained in vivo. According to such constructs, paired LH-T time series in healthy men exhibit evidence of dose-response downregulation within individual T pulses. Statistical model justification was highly significant (*P* < 10^{−4}) compared with a no-hysteresis construct. In addition, the sensitivity-downregulation construct revealed distinct age-associated adaptations in the male gonadal axis. Foremost predictions were lower LH efficacy [asymptotically maximal T secretion due to LH action (*P* = 0.004)] and greater absolute downregulation of T response sensitivity to LH [sensitivity-shift size (*P* < 0.001)] in older than young men. If verified in larger and prospective cohorts, these outcomes would point to dynamic pituitary-gonadal effector-response coupling in vivo, which may be influenced by one or more factors associated with healthy aging.

The accompanying studies differ from earlier model analyses in several important ways. *1*) The methodology extends the recent concept of allowable (estimable) dose-response hysteresis (26) to a major regulatory interface in the male gonadal axis. *2*) Feedforward LH → T dose-response functions are estimated conditional on deconvolution-based reconstruction of LH concentrations and T secretion rates, as respective input (effector) and output (response) signals (6). *3*) The deconvolution methodology is integrative and permits two waveforms (shapes) of LH and T secretory bursts. The latter allowance is based on an idea presented recently for the ACTH-cortisol axis (27). *4*) According to such analyses, older age is associated with diminished LH efficacy and, apparently, heightened T sensitivity and LH potency. The last two preliminary findings could be consistent with the broader concept of target organ denervation hypersensitivity, which is recognized in neuronal, vascular, renal, muscular, and other adaptive physiological systems (15, 22, 38, 41, 44, 49, 50). A relative LH-deficiency state was inferred here in older men by their lower mass of LH secreted per burst (Fig. 2; *P* = 0.007). In this regard, increased Leydig cell responsivity (decreased EC_{50} of exogenous LH stimulation) has been observed in vitro after short-term in vivo LH deficiency in the adult male rat (17) and in the presence of adenosine in Leydig tumor cells (11). Selective enhancement of target tissue responses is also evident in glucose-insulin (48), calcium-parathyroid hormone (53), and other endocrine interfaces (15, 44, 49, 60).

Whether the primary mechanism mediating reduced LH efficacy/lower T secretory responses to inferably maximal concentrations of LH in vivo in aging volunteers is decreased LH biopotency or diminished Leydig cell responsivity is not known definitively. Indeed, older men have normal or decreased in vitro LH bioactivity (13, 40). Attenuated LH efficacy would be consistent with diminished T responses to pharmacological hCG stimulation (54, 55) and with inferences in the aged male rat (18, 42). However, infused hCG is distinctly nonphysiological in the male because of lutropic dose and binding kinetics (12, 34). In the present uninjected physiological state, the sensitivity-hysteresis model predicted a 2.7-fold reduction in median LH efficacy under endogenous LH-pulse drive: T was secreted at 43 and 16 ng·dl^{−1}·min^{−1} in young and older men, respectively (*P* = 0.004). Earlier methods of dose-response reconstruction inferred lower absolute LH efficacies in men (e.g., 3.8–5.8 ng·dl^{−1}·min^{−1}) (32, 33). On the basis of available direct measures of T kinetics (57), the nominal T distribution volume is 25 liters (250 dl) per square meter of surface area. If so, calculated initial (onset) LH efficacy (milligrams of T secreted per day per 1.73 m^{2}) by age stratum would be 26 and 9.9 in young and older men, respectively. These asymptotic estimates are consistent with the severalfold augmentative effect of pharmacological doses of hCG on unstimulated daily T secretion rates of 4–8 mg in healthy young individuals (40). However, no studies have used randomly ordered pulses of recombinant human LH to estimate directly LH-T dose-response properties or maximal T secretion rates in men.

The median time delay of 23 min from T secretory-burst onset to hysteretic downregulation would be consistent with rapid (within 20–30 min) pharmacological desensitization observed in several other endocrine systems in vivo and in vitro (2, 3, 16, 21, 35). Here, reversible cycles of desensitization were consistently inferable (in 25 of 26 subjects) under physiological conditions for the endogenous LH-T system. Consistently inferable downregulation in all three model forms would be consonant with the capability of submaximally stimulatory LH concentrations to induce Leydig cell response refraction in some in vitro studies (19, 20).

In summary, the aptness of the three complementary analytical models of reversible hysteresis in framing the dynamics of pituitary gonadotropin-gonadal sex steroid coupling suggests that this type of quantitative methodology should find application in other physiological contexts, in which brief, recurrent, and reversible agonist-response adaptations are postulated.

### Perspectives and Significance

Demonstration of unperturbed effector-response cycles under physiological conditions in vivo is inherently difficult. Conventional approaches include administration of a pharmacological effector at two separate times so as to compare initial and delayed responses and/or mathematical simulation of potential mechanisms of tachyphylaxis, hypersensitivity, or tolerance (37, 47). Neither of these strategies necessarily recapitulates underlying physiological effector-response adaptations. If time-dependent response adaptations are assumed to reflect altered dose-response dynamics, then the presently proposed analytical estimation model may be appropriate. The objective is to obtain noninvasive, uninjected, nonpharmacological estimates of recurrently adaptive effector-response properties. Because response desensitization and response hypersensitization accompany G protein-coupled signaling (4, 10), quantification of successive effector-response cycles in such systems may be especially informative to physiologists. Accordingly, the concept of rapid in vivo hysteresis should provide a basis for comparable analyses in many other regulatory and integrative physiological systems.

## GRANTS

This work was supported in part by National Center for Research Resources Center for Translational Science Activities Grant 1 UL 1 RR-024150 and National Institute on Aging Grant AG-031763.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## ACKNOWLEDGMENTS

We thank Donna Scott and Samuel Rudolph for manuscript preparation, Ashley Bryant for data analysis and graphics, the Mayo Immunochemical Laboratory for assay assistance, and the Mayo research nursing staff for protocol implementation.

## APPENDIX A: Three Estimation Models

In many biological systems, in vivo elimination is not well represented by a single-exponential decay. Rather, there are often two components: *1*) a fast term reflecting advection, diffusion, and mixing in the blood and *2*) a slow component embodying irreversible metabolism or clearance (58). The combined kinetics can be described by two coupled differential equations
*X*(*t*) denotes the concentration at time *t*, with *X*^{(1)}(*t*) for the fast component and *X*^{(2)}(*t*) for the slow component. The fast and slow rates of elimination are α^{(1)} and α^{(2)}, and the fast and slow fractions are *a* and (1 − *a*), respectively, while *s* is dummy variable. *Z* denotes the secretion rate. Because data are not observed continuously, but at some sampling rate (Δ*t*), observations are made at times (*t*_{1}, *t*_{2},…, *t*_{n}), *t*_{i + 1} − *t*_{i} = Δ*t*. The two models are discretized (can be evaluated at the sampling times) as follows: by a difference equation
^{(1)} = exp[−α^{(1)}Δ*t*] and α̃^{(2)} = exp[−α^{(2)}Δ*t*]. Virtually all endocrine modeling has utilized the difference equation (*Eq. 3*), rather than the integral equation (*Eq. 4*), in estimation. The integral equation model subsumes the difference equation model; hence, it is more comprehensive.

In the difference equation approximation, the fractions [1 − α^{(1)}Δ*t*] and [1 − α^{(2)}Δ*t*] are often used in place of α̃^{(1)} and α̃^{(2)}, representing approximations to derivatives at zero; such use is appropriate when Δ*t* is small. An advantage of using α̃^{(1)} and α̃^{(2)} is that the difference and differential equation models are constrained to have the same half-lives of elimination (or as close as possible for Δ*t* sampling).

When Δ*t* is small compared with the half-lives and/or the rapidity of change in *Z*, there is not much difference between the difference and integral models. However, when this is not the case, inaccurate estimates can occur. Specifically, *Eq. 3* compared with *Eq. 4* does not account for time-varying *Z* and elimination that are, in fact, taking place within the time interval (*t*_{i}*,t*_{i + 1}) defining consecutive observation times. It is as if elimination were occurring only at time *t* and secretion were constant at the initial value *Z*(*t*_{i}) over the interval. When the sampling interval is similar to or longer than the half-life of elimination, inaccurate estimates occur, since time-varying secretion and elimination within the interval are ignored.

The difference equation model (*Eq. 3*) can be derived from the discrete time integral equation model (*Eq. 4*) by assuming a constant secretion inside the integral over the interval at the initial value *Z*(*t*_{i}) and ignoring elimination (i.e., α is taken as 0 in the exponent). An intermediate model (*Eq. 5*), which lies between *Eqs. 3* and *4*, can be constructed by assuming constant secretion inside the integral at the initial value *Z*(*t*_{i}), but allowing constant-rate elimination over the interval in the following intermediate model

One might propose using the integral equation model (*Eq. 4*) exclusively. The difficulty is that the numerical integrations involved in *Eq. 4* make the computations beyond the reach of the resources of most investigators, especially when one factors in the selection of pulse times. An alternative is to frame the basic procedure as a sequence of three stages, as done here. *Stage I* employs the difference equation model (*Eq. 3*); the parameters are estimated, and the AIC-penalized maximum-likelihood pulse times are chosen. *Stage II* utilizes intermediate *Eq. 5* to create more accurate slow half-life estimates and serve as a transition to *stage III*. *Stage III* uses the integral equation model (*Eq. 4*), which is the most accurate for estimating a variable secretion rate function (*Z*).

## APPENDIX B: ESTIMATION

Since *Z* is predicated on release occurring at certain (unobserved) pulse times in a given hormone concentration profile, a previously published pulse-detection algorithm was applied (24). The algorithm is a nonlinear diffusion equation that selectively smoothes the profile using a diffusion coefficient that is inversely related to the degree of positivity of the local gradient in the concentration profile. One begins with all local minima: those followed by a rapid or marked increase are smoothed minimally, whereas those of slow or little increase are smoothed more. The smoothed profile is not used in parameter estimation; it is used only in construction of possible pulse times. As the algorithm proceeds, a sequence of diminishing numbers of pulse times is constructed, denoting increasing susceptibility of the pulse time to smoothing: P = [*P*^{(1)}, *P*^{(2)},…, *P*^{(M − 1)}, *P*^{(M)}], where *P*^{(1)} is the most robust and *P*^{(M)} is the least robust (most diminutive and easily smoothed) pulse time.

Putative pulse sets are constructed by including pulse times with indexes up to *m,* P_{m} = [*P*^{(1)}, …, *P*^{(m − 1)}, *P*^{(m)}], 1 ≤ *m* ≤ *M*. There are *M* such sets, wherein pulse times in any given set *P*_{m} are not arranged in increasing order of sampling time, but in order of susceptibility to smoothing (removal). Thus, denote the pulse times of *P*_{m}, rearranged in increasing order of time, as T^{m} = [*T*^{(1)},…, *T*^{(m − 1)}, *T*^{(m)}]. Parameter estimation then proceeds by penalized MLE (PMLE) conditional on each pulse-time set, where penalization is on the number of pulse times *m*. The AIC and Bayesian information criterion can be used as model comparators.

Consider a given fixed pulse set T^{m} = [*T*^{(1)},…, *T*^{(m − 1)}, *T*^{(m)}], and describe the waveform of mass release (realization of secretion rate over time) by a three-parameter generalized gamma density (24). For thyroid-stimulating hormone and corticotropin, there is often a day-night difference in the waveform, one for day (D) and one for night (N) (25, 31). If daytime (to be estimated) is defined by the interval (ϕ_{1},ϕ_{2}), then the waveform at pulse time *T*^{(k)}, depending on daytime or nighttime, is

The secretion rate is given as the sum of a nonpulsatile (basal) secretion rate (β_{0}) and a train of secretory bursts (pulses). Each secretory-burst mass is described as a weak linear function of the preceding interpulse interval [*T*^{(k)}–*T*^{(k − 1)}] plus a random effect [*A*^{(k)}]. The latter allows for variations in burst size, which are not explicitly modeled by the linear function (e.g., variability due to nonuniform access of secretory granules to capillaries from burst to burst). Thus, overall secretion is

Let θ denote the secretion and elimination parameters for the model, wherein fast and slow rates of elimination α^{(1)} and α^{(2)} are to be estimated, and their fractions [*a* and (1 − *a*)] are taken as literature-based population values. The *m* random effects (one for each burst) are assumed to be independently and normally distributed (IID normal) with mean zero and variance σ_{A}^{2}. For a given pulse set T_{m} = [*T*^{(1)}, …, *T*^{(m − 1)}, *T*^{(m)}], a parameter choice θ, and the secretion rate function *Z*, any one (or all) of the models (*Eqs. 3*–*5*) can be applied, resulting in the (true) concentrations with measurement error
_{ε}^{2}. Because the random effects (*A* values) in the secretion rate (*Z*) enter linearly, the resulting likelihood for the observed concentrations (*Y* values) is Gaussian. An AIC (or Bayesian information criterion) penalty term is appended, penalizing the number of pulse times *m*, resulting in a penalized log-likelihood function that is maximized over θ and *m* as

The optimization is constrained, since any estimate θ̂ must result in a nonnegative *Z*; i.e., the conditional expectation of the secretion rate evaluated at θ̂ must be nonnegative: E_{θ̂}[*Z*(*t*_{i}), *i* = 1, …, *n*|*Y*_{i}, *i* = 1, …, *n*] ≥ 0. This requires calculating the conditional expectations of the *A* values, conditioned on the *Y* values: [*Â*^{(1)}, *Â*^{(2)}, …, *Â*^{(m)}] = E_{θ̂}[*A*^{(k)}, *k* = 1,…, *m*|*Y*_{i}, *i* = 1,…, *n*]. Values are estimated from the definition of *Z* (*Eq. 7*), and nonnegativity is assessed. The results of the PMLE are θ̂, *m*̂, and the AIC-optimal pulse-time set P_{m̂} = [*P*^{(1)},…,…,*P*^{(m̂ − 1)}, *P*^{(m̂)]}.

In previous analyses (39), allowable pulse-time sets included only the *M*-specific collection of decreasing sets. The possibility arises that one or more of the *m*̂ pulse times removed could enhance the fit if added to the resulting AIC-optimal set, *P**m*̂. Such sets were excluded, because smoothing yielded the specific decreasing collection P_{j}, *j* = 1,…, *M*. To address this possibility, one may evaluate the effect of addback of previously removed pulse times, i.e., by creating new pulse-time sets {P̃*m*̂, *P*^{(k)}}, *k* = *m*̂ + 1,…,*M*. For each of the *M*-*m*̂ such sets, one calculates the PMLE, pl[θ|*Y*_{i}, *i* = 1,…, n; {P_{m}, *P*^{(k)}}], which is compared with that achieved by *P**m*̂via AIC. The best of such qualifying single pulse times is added, if one exists. The procedure is then repeated for the remaining *M*-(*m*̂ − 1) pulse times, until none of the potential addback pulse times is preferred, allowing for the AIC penalty. The aggregate of these procedures constitutes the estimation algorithm.

#### Estimation Algorithm

Observed data: *Y*_{i} = *X*(*t*_{i}) + ε(*i*), *i* = 1,…, *n*.

##### Stage I.

Starting parameter θ_{*} and putative pulse sets P_{m}, *m* = 1,…, *M.*

Difference equation model: *X*(*t*_{i} + 1) = [*a*α̃^{(1)} + (1 − *a*)α̃^{(2)}]*X*(*t*_{i}) + Δ*tZ*(*t*_{i}).

Penalized log-likelihood function: pl(θ|*Y*_{i}, *i* = 1,…, *n*; P_{m}).

Result: PMLE θ̂^{(stage I)} and AIC-optimal pulse set: *P*̃_{m̂ + s} = { *P*_{m̂}, *P*^{(k1)}, *P*^{(k2)},…, *P*^{(ks)}}, including addback of *s* pulse times.

##### Stage II.

Start with results from *stage I*, θ̂^{(stage I)}, P̃_{m̂ + s}.

Intermediate model replaces the difference equation model in penalized log-likelihood function

Penalized log-likelihood function: pl(θ|*Y*_{i}, *i* = 1,…, *n*; P̃_{m̂ + s)}

Result: PMLE θ̂^{(stage II)}; pulse times in P̃_{m̂+s} allowed to shift (slightly).

##### Stage III.

Start with results from *stage II*, θ̂^{(stage II)}, P̃_{m̂+s}. The integral equation model replaces the intermediate equation model in penalized log-likelihood function

Penalized log-likelihood function: pl(θ|*Y*_{i}, *i* = 1,…, *n*; P̃_{m̂+s)}.

Result: PMLE θ̂ [=θ̂^{(stage III)}]; pulse times in P̃_{m̂+s} allowed to shift (slightly).

##### End of Algorithm

When the pulse times P̃_{m̂+s} and the final parameter estimates θ̂ are obtained, the secretion rate function *Z* can be calculated (at the observed sample times) as conditional expectations evaluated at θ̂

One can now calculate the model fits, or the predicted concentrations. Fits are obtained by a convolution (*Eq. 4*) using the estimated secretion rate (*Eq. 7*) and estimated biexponential kinetics [α̂^{(1)} and α̂^{(2)}], two components of θ̂. Subscript r denotes reconvolved concentrations

## APPENDIX C: APPLICATION TO LH AND T ANALYSES

To illustrate the sequential three-stage estimation algorithm, Fig. 8 depicts the three main equations used in the program. Figure 9 applies these analytical steps using an illustrative LH time series obtained in a young male subject. Figure 10 presents corresponding stepwise deconvolution data for T obtained in an older man. These two profiles typify the spectrum of profiles observed here. Figure 11 shows more detailed results from the integral-deconvolution stage in the same young man (LH) and older man (T).

## APPENDIX D: LH FEEDFORWARD ON T: DOSE RESPONSE WITH HYSTERESIS

Methods described in appendixes a–c yield *1*) estimated (integral model) reconvolved LH concentrations: *Ŷ*_{L,i}, *i* = 1,…, *n* (*Eq. 11*), which will be the basis for the LH feedforward signal on T secretion, and *2*) estimated T secretion rates: *Ẑ*_{T,i}, *i* = 1,…, *n* (*Eq. 10*), which will be the output of the dose-response function. To estimate the dose-response function, there are two further considerations: *1*) a varying time delay between the onset of an LH pulse and the observed T secretory response and *2*) potentially variable loss of responsiveness of T to the LH feedforward signal. If one wishes to accurately recover the dose-response structure, both of these considerations must be addressed.

Utilizing the fitted (reconvolved) LH concentrations: *Ŷ*_{L,i}, *i* = 1,…, *n,* the LH feedforward signal F_{L}(*t*) was constructed piecewise, from one T pulse-onset time to the next. This allows one to account for varying time lags between the onsets of LH and T pulses. For each T pulse time *T*_{T}^{(k)}, the LH pulse *T*_{L}^{(j)} nearest within the allowable time, [−60,10] min, was identified. If an LH pulse exists, it is shifted to the T onset point, so that the two onset points are aligned. If no such pulse exists within the time interval, then a time lag of 40 min was applied between the LH and T pulse-onset times (32). The results of this phase of the dose-response analysis are shown in Fig. 4, *bottom*.

Estimation proceeds by using estimated T secretion rates *Ẑ*_{T,i}, *i* = 1,…, *n* assumed to be *Z*_{T}(*t*_{i}) + error. Random effects in efficacy (*A* values) are included to accommodate pulse-by-pulse size variability. To allow for possible loss of responsiveness to LH drive, a time-delayed shift in the dose-response function is permitted. The shift is assumed to occur at time *M*_{LonT} (to be estimated) following the T secretion pulse onset. From pulse onset until time *M*_{LonT}, the LH-T relationship is defined by one dose-response curve and, thereafter (until the next T pulse onset), by a shifted curve. Hence, a hysteresis-like phenomenon occurs: the recovery phase of the dose response shifts within a T pulse after an estimable delay and then resets to the initial curve at the next secretory pulse (Fig. 1). Dose-response hysteresis types can be evaluated by three models: *model 1* [one-half maximally effective stimulus concentration (LH potency)]

*model 2* [dose-response slope (testis sensitivity)]

*model 3* [asymptotic maximum (LH efficacy)]

Each model is fit separately, allowing subsequent comparisons of model residual error (see *Statistics*). In addition, generalized likelihood-ratio testing was used to evaluate each hysteresis model compared with the no-hysteresis construct.

## Footnotes

↵1 Hysteresis is defined as dependence of a dose-response function on the time evolution of effector concentrations, such that, after (compared with before) a hysteresis inflection point, the magnitude of the response decreases or increases at any given effector concentration.

↵2 Onset (initial) dose-response parameter is defined as potency, sensitivity, or efficacy dose-response parameters estimated using paired effector-response data beginning at the onset of a response and continuing until the time of hysteretic inflection.

↵3 Recovery (delayed) dose-response parameter is analagous to onset parameter, except estimation is based on the time window beginning with the hysteretic inflection and continuing until the next response onset.