## Abstract

Neuroendocrine axes are feedback- and feedforward-coupled dynamic ensembles. Disruption of selected pathways in such networklike organizations may explicate loss of orderly hormonal output as observed in aging. To test this notion more explicitly, we implemented an earlier computer-assisted biomathematical model of the interlinked male hypothalamo [gonadotropin-releasing hormone (GnRH)]-pituitary [luteinizing hormone, (LH)]-testicular [Leydig cell testosterone (Te)] axis (*Am J Physiol Endocrinol Metab Physiol* 275: E157–E176, 1988; Keenan D., W. Sun, and J. D. Veldhuis, *SIAM J Appl Math* 61: 934–965, 2000). Thereby, we appraise mechanistic hypotheses for more disorderly LH and Te secretion in aging men. We compare model predictions with monitored abnormalities in the older male, namely, irregular patterns of individual and synchronous LH and Te release, reduced 24-h rhythmic Te output, and variably elevated LH secretion. Among the mechanisms examined, the most parsimonious aging hypothesis would entail impaired LH feedforward on Te without or with attenuated Te feedback on GnRH/LH secretion. This investigative strategy should aid in exploring new postulates of disrupted feedback networks in pathophysiology.

- gonadotropin-releasing hormone
- luteinizing hormone
- Leydig cell testosterone

in aging primates,reproductive quiescence emerges in midlife in the female, and reproductive potential declines during the later lifetime in the male (1, 4, 5, 8, 26, 31). Clinical analyses in aging men have disclosed progressive age-related decrements in testosterone's bioavailability (so called andropause) and somatotropic-axis output (so called somatopause) (4, 5, 30, 35). The pathophysiological mechanisms underlying waning activity of the gonadotropic and somatotropic axes are not known but likely contribute to the frailty that accompanies aging (21-24).

Dynamic properties of the aging male reproductive axis show evident disruption before any reproductive loss occurs (36). For example, the coordinate release of luteinizing hormone (LH) and testostrone becomes markedly asynchronous in healthy older men despite normal (young adult) serum concentrations of LH and testostrone (12, 16, 27). The joint secretion patterns of ACTH and cortisol are likewise disturbed in aging humans despite normal daily secretion rates of both hormones (19). These data suggest the postulate that intra-axis feedback-dependent (network level) control of neurohormone secretion is impaired in the earlier stages of aging.

Deterioration of bihormonal secretory synchrony in a neuroendocrine ensemble (above) could denote regulatory defects in interglandular pathway communication (3-5, 18, 20, 32, 34, 36). However, interpreting or predicting a priori the precise site, nature, relevance, and/or magnitude of a postulated linkage defect in a complex interconnected axis (6, 7, 26) is difficult on intuitive grounds alone. Here, we implement an earlier computer-assisted biomathematical model (6, 7, 26) to test selected hypotheses for disrupted control of the GnRH-LH-testesterone feedback axis in the healthy aging male.

## METHODS

### Introduction

The male reproductive axis in a simplified construct consists of three primary control signals: hypothalamic gonadotropin-releasing hormone (GnRH), pituitary LH, and gonadal testosterone (Te). Dose-response relationships in the male system have been largely defined for individual nodes acting in isolation [e.g., GnRH's stimulation of LH secretion, LH's dose-dependent stimulation of Te secretion (26), etc.]. However, the implicitly dynamic nature of this network arises from time-lagged feedback and feedforward interactions among all three interconnected loci (6,7). These connections are illustrated schematically in Fig.1. Here, we use this biomathematical notion to compare various hypotheses of feedback and/or feedforward modifications due to aging.

As formalized in detail in Refs. 6 and 7, we will denote the three primary hormone concentrations [GnRH (G), LH (L), and Te] in blood as: X(t) = [X_{G}(t), X_{L}(t), X_{Te}(t)] and their respective rates of secretion as: Z(t) = [Z_{G}(t), Z_{L}(t), Z_{Te}(t)]. Assuming simple conservation of mass, we begin with
Equation 1The rate of hormone elimination is represented as proportional to its concentration [αX(t)]. This has been experimentally observed by injecting purified hormone into the blood and monitoring its exponential disappearance rates in hormone-deficient volunteers (26).

An earlier more comprehensive formalism (6, 7) relates how each rate of secretion: Z_{G}(t), Z_{L}(t), Z_{Te}(t) depends on the history of relevant hormone concentrations within the interconnected feedback system {[X_{G}(s), X_{L}(s), X_{Te}(s)], s ≤ t}.
We assume that GnRH feedforward on LH, LH feedforward on Te, and Te feedback on GnRH or LH are exerted via a time-delayed (time delay = [l_{1}, l_{2}]) and time-averaged effect of the relevant hormone concentration (6, 7). Thus, at time t and 0 < t − l_{2}≤ t − l_{1} < ∞, the feedback signal intensity (i.e., ƒ
denotes a time average) is
We consider that feedback is exerted on the relevant control nodes (system components) in Fig. 1 via putatively monotonic dose-response logistic functions (denoted by H(·)′s); see Fig. 2. For the simplified male GnRH-LH-Te axis, we allow four major feedback/feedforward dose-response functions: H_{1,2}(·), which describes the GnRH pulse firing rate as a joint function of Te (feedback) and GnRH (autofeedback) concentration, H_{3}(·), which gives the rate of GnRH pulse-mass accumulation as a function of Te (feedback) concentration, H_{4}(·), which defines the rate of Te secretion as a function of LH (feedforward) concentration, and, H_{5,6}( · , · ) for the rate of accumulation of the LH pulse mass as a function of Te (feedback) and GnRH (feedforward); see Figs. 1 and 2. The subscripts correspond to feedback/feedforward relationships [(1)–(7)] given in below; the time delays for the j th relationship are denoted as (l_{j,1}, l_{j,2}).

### GnRH-LH-Te Feedback/Feedforward Interactions

#### Hypothalamus: Te and GnRH Feedback [interactions (1)–(3)] on GnRH.

The blood Te concentration (ng/dl) is assumed to exert a negative time-delayed (nominally 80–90 min) feedback (11, 12), and GnRH concentration (pg/ml) a (auto) negative time-delayed (approximately 1–1.5 min) feedback (2), on the GnRH pulse firing rate (given in pulses/day) (see Figs. 1 and 2)
Equation 2The Te concentration (ng/dl) is allowed to exert negative and time-delayed (e.g. approximately 80–90 min) (26) feedback on the rate of GnRH hormone synthesis (pg · ml^{−1} · h^{−1}) (see Figs.1 and 2)
Equation 3

#### Testes: LH feedforward [interaction (4)] on Te secretion.

The blood LH concentration (IU/l) is viewed as producing a positive time-delayed (∼20–30 min) (12, 26) feedforward action on Te rate of hormone synthesis (ng · dl^{−1} · h^{−1}) (see Figs.1 and 2)
Equation 4

#### Pituitary: GnRH feedforward and Te feedback [interactions (5)–(7)] on LH secretion.

The portal GnRH concentration (pg/ml) is presented as a positive time-delayed (∼0.5–1.5 min) feedforward signal driving the rate of pituitary LH hormone (mass) synthesis (IU · l^{−1} · h^{−1}) (2) (see Figs. 1 and 2). Conversely, the plasma Te concentration (ng/dl) exerts a negative time-delayed (nominally 80–90 min) feedback on rate of LH hormone synthesis (IU · l^{−1} · h^{−1}) (26).
Equation 5where T
, T
, … are the LH pulse times, and
The function Γ(·) is a sum of exponentials, which thus describes the expected cascade of cellular reactions (2nd and 3rd messengers, etc.) initiated by the GnRH signal acting on LH synthesis and storage for release (7).

The times of pituitary LH pulses are envisioned to mirror hypothalamic GnRH signals with a slight time delay τ_{L} (e.g., 1–2 min) and with a brief refractory period r_{L}(e.g., ∼10–15 min) after a GnRH pulse time (2). This refractory period is the interaction (7).

We model these feedback interactions [except for (7)] via monotonic logistic dose-response functions
Equation 6
or where the coefficients themselves are described by logistic functions, e.g.,
If B_{i} > 0, the feedback is positive (i.e., feedforward effect); if B_{i} < 0, the feedback is negative. For each hormone, the resulting nonbasal rate of synthesis S(·) will depend on time-delayed feedback from the system through such interface functions H(·), plus noise ξ(·). The H(·) defines a “conditional expectation,” or an average rate for all the cells in the endocrine gland (see Figs. 1and 2); the actual realized rate of synthesis S(·) varies about H(·). Experimental evidence suggests that the noise varies about the conditional expectation in a stationary manner; we have formulated the noise ξ(·) similar to an Ornstein-Uhlenbeck process, except that we have constrained the resulting rate of synthesis S(·) to be nonnegative and bounded (7). It is important to allow for such noise as ξ(·) in the realized rates of synthesis, because in this specific context, one could explore whether the young and older male reproductive systems differ because: *1*) the average rate of hormone synthesis H(·) is altered by age and/or *2*) with age, there is greater variability about this average rate H(·). Other pathway-specific postulates of aging are evaluated below.

### Circadian Rhythm

The mean serum testosterone concentration in young men typically decreases during the late day and rises in the early morning (9-12, 24, 26, 28, 30-32). The mechanisms subserving this diurnal rhythmicity are not well understood. In a simplified construct, the nyctohemeral rhythm could be viewed as a modulation of LH's feedforward on the rate of secretion of Te, H_{4}(·), by a strictly positive, 24-h periodic function μ(·). In the present simulations, the μ(·) function was assumed among other possibilities to consist of one harmonic with amplitudes B_{0}, B_{1}, and phase shift φ_{1} being chosen so μ(·) would vary between a high at 4 AM and a low at 4 PM (e.g., B_{0} = 1.0, B_{1} = 0.35, φ_{1} = 4.0, starting at 8 AM; individual variability allowed, as discussed below
Equation 7

### Construct of the Male GnRH-LH-Te Axis

According to the foregoing background and detailed motivation given in (6, 7, 26), the ensemble formulation of the male axis can be given by the following system of (stochastic differential) equations.

#### Pulse generator.

The pulse times of GnRH: T
, T
, T
, … are the result of a (point) process with an intensity
Equation 8
where λ(t)dt is interpretable as the instantaneous probability of a pulse in the infinitesimal time increment (t, t + dt). The conditional density for T
given T
and λ(·) will be required to satisfy
Equation 9resulting in GnRH pulse times T
, T
, … , and LH pulse times (due to the refractory period r_{L} and time-delay τ_{L}):
and the associated counting processes:
and
This model for the pulse generator, based on assumed loosely coupled GnRH neuronal firing within an ensemble of excitable hypothalamic neurons, is: *1*) dynamical, in the sense that the probabilistic structure of the pulse times depends on feedback within the system, and *2*) flexible, wherein a Weibull renewal process (corresponding to a constant λ) allows for variability of the interpulse interval lengths, independent of the mean, via the parameter γ ≥ 1 (6, 7).

#### Agonist-driven (nonbasal) rates of synthesis [S(·)] for GnRH, LH, and Te.

(feedback from Te, plus allowable variability) (feedforward from LH, plus allowable variability) (feedforward/feedback from G and Te, plus allowable variability) (accumulation of newly synthesized GnRH and LH in secretory granules, and where η represents the time rate of release of newly synthesized granules.)

#### Pulse masses M^{j} for GnRH and LH.

(normalized rates of secretion) (fractional mass M remaining at time T ) (Available mass equals that not released from previous pulse plus newly synthesized granules.)

#### The male axis (GnRH, LH, Te) core equations.

Combining the above (below, x_{+} denotes max {x, 0})
Equation 10
where the W_{i}(·) 's are independent standard Brownian motions and ς_{i}dW_{i}(t), subscripts i = G, L, Te, describe the in vivo biological variability due to the diffusive behavior (mixing, advection, turbulence, etc.) of the neurohormone within the blood. The values of ς_{i}, i = G, L, Te, were chosen so that the resulting coefficients of variation were approximately 6%. Te synthesis and secretion are assumed to be coupled (without an intervening response cascade) by simple diffusion of the synthesized steroid, because no granular storage form exists for Te (unlike LH) (26).

#### Allowance for variation among individuals.

To allow for the variation among individuals, certain parameters may vary randomly about a population value. In particular, we allow the Te basal secretion β_{Te} (*Eq. 10
*) to be uniformly distributed between 1,000 and 1,500 ng · dl^{−1} · h^{−1}, the circadian rhythm amplitude B_{1} (*Eq. 7
*) to be uniformly distributed between 0.2 and 0.5 with B_{0} = 1, the pulse regularity parameter γ (*Eq. 9
*) to be uniformly distributed between 2 and 6, and the slope of the Te feedback on the instantaneous rate of GnRH pulsing (*Eqs. 6, 8,* and *
9
*) to be uniformly distributed between −0.005 and −0.009 (6, 7, 26).

What one would actually observe in most experiments is a discrete-time sampling of LH and Te (ordinarily not GnRH) with measurement error due to sample collection, processing, and assay Equation 11

### Overview of Approximate Entropy (and Cross-Approximate Entropy)

The approximate entropy (ApEn) statistic was developed to quantify the degree of irregularity, or disorderliness, of a single time series. Analogously, cross-ApEn quantifies joint pattern synchrony between two time series. The methods have broadly applied to both medical and nonmedical data (14, 15, 17).

For a given data series Y_{k}, k = 1, … , N, ApEn is calculated as follows: the (sample) standard deviation of Y, SD (Y), is calculated, and two values, m (pattern length) and r (discriminating threshold), are specified. The validated choices for m and r for neurohormone series of the length sampled here are m = 1, r = 0.2 SD (Y). One then calculates
Equation 12

#### Removing trends before ApEn calculation.

The formulations of ApEn and cross-ApEn apply to times series that are at least asymptotically stationary (14, 15, 17). One direct consequence of this assumption is an approximately constant marginal standard deviation SD (Y), which is used as part of the yardstick [r = 0.2 SD (Y)] to monitor significant pattern differences. If a series is not stationary in its mean, then one would need to remove an estimated mean function, so that the resulting estimated standard deviation is not inflated. In the case of Te, which can exhibit prominent rhythmicity over 24 h (and to a lesser extent, LH), any major effects of a circadian rhythm in the concentration should thus be minimized before ApEn and cross-ApEn are computed. In *Application to Clinical Data,* ApEn and cross-ApEn are calculated for clinical data from older men and young men; Fig. 3 (*left*) displays a Te concentration (solid line) profile of a young male sampled every 10 min over 24 h (starting at 8 AM). The effect of the circadian rhythm on Te concentrations and its (sample) SD are apparent, which in part is a consequence of the nonstationarity in the mean.

#### Periodic rhythmicity.

If the effect of the circadian rhythm were exerted purely via a cosine function, then its removal by an estimated cosine function would probably suffice. We tested this notion via trigonometric regression at the first discrete Fourier frequency λ = (2π/n), which will provide the least-squares fitted cosine function, Acos(λt) + φ. The periodogram I^{(n)}(·), a normalized version of the squared modulus of the Fourier transform of a series {Y_{k}} provides the squared amplitude at this first frequency and all the other discrete Fourier frequencies
Figure 3 displays the fitted cosine function superimposed (interrupted line) on the Te concentrations (Fig. 3,*top left*). The cosine-detrended series is shown for comparison (Fig. 3, *middle left*). The periodogram of the original Te concentrations is compared on Fig. 3, *top right.*The sample SD of the cosine-detrended series is 65, which is a reduction from 79 of the original series. There appear to be residual nonstationarities after cosine detrending: see Figure 3, *middle right,* wherein the periodogram of the cosine-detrended series is superimposed (interrupted line) on the periodogram of the original series. One could attempt to remove these additional low frequencies more systematically, as suggested in the following subsection.

#### Mean function estimation via smoothing by the heat (or diffusion) equation.

If modulation of the 24-h neurohormone profile should occur in a more complicated nonlinear manner, then the effects would probably be smeared over a low-frequency band. One way to remove a complicated nonlinear trend, without necessary knowledge of its form a priori, is by implementing the heat equation. In this strategy, the data serve as the initial conditions with Dirichlet, Neumann, or mixed boundary conditions. We have used Dirichlet conditions and subtracted the resulting smoothed version from the data Y_{k}, with
(the mean) then added back. This procedure removes the low-frequency band and thus achieves a “straightening out,” without making a priori assumptions about its form (see more detailed motivation in
).

Suppose that the data Y_{k} are obtained by sampling a function {Y(t), 0 ≤ t ≤ 1}, with t being observational time (units of days): Y_{k} = Y((k − 1)/(N − 1)), k = 1, 2, … , N. Let U be the solution of the following equation, with accompanying initial and boundary conditions
The variable s represents algorithmic time. Except for the boundary conditions, one can envision the algorithm (at algorithmic *time *s) as convolving Y(·) with a Gaussian kernel whose variance is C^{2}s. The equation is implemented for 0 ≤ s ≤ S (see
), where the terminal function is [U(t, S), 0 ≤ t ≤ 1], and the resulting smoothed (and discrete) version of the data Y_{k} isŶ_{k} = U((k − 1)/(N − 1), S), k = 1, 2, … , N. This procedure is minimally dependent on the precise nature of the underlying slow trends. Figure 3,*top left,* shows serial serum testosterone concentrations with the superimposed heat equation-smoothed series (dashed-dotted line). Figure 2, *bottom left,* gives the resulting series with the heat equation-smoothed component removed. For comparison, Figure 3, *bottom right,* presents the periodograms of the original (solid line) and heat equation-smoothed Te series (interrupted line). Note that only low-frequency terms are removed. The sample SD of the heat equation-detrended series fell to 43 in this case (compared with a value of 79 in the original and 65 in the cosine-detrended series). The resulting ApEn values for these three Te series were 1.5 (original), 1.6 (cosine detrended), and 1.8 (heat equation detrended).

## COMPUTER-ASSISTED EXPERIMENTS

On the basis of the mathematical construct of the GnRH-LH-Te axis given in *Construct of the Male GnRH-LH-Te Axis* (above and Refs. 6, 7, 26), we performed selected computer-assisted experiments, as described below. Each consisted of 500 simulations implemented in Matlab.

#### Normal young male (experiment 0).

Published young-adult male 24-h serum LH and Te concentration time series (12, 18) are emulated realistically by the foregoing biomathematical construct (6, 7). The H dose-response functions for the normal young male are displayed in Fig. 2
*A.* Figure 4displays one realization from each of the six (*numbers 0–5*) experiments, with the randomness in each realization being initiated with the same seed. These profiles offer a visual impression of the dynamics produced by modifying a given dose-response H function. Figure 5 gives the resulting histograms of the 500 values of (heat equation detrended) LH ApEn, Te ApEn, and the cross-ApEn. And, Fig.6 presents the histograms of the 500 LH and Te means, and the ratio of the 24-h amplitude/mean of the simulated Te profiles for each of the experiments. We obtained comparable histograms in the case where there was no cosine modulation of LH's feedforward on Te (i.e., B_{1} = 0 in*Eq. 7
*) (not shown). Thus comparable model-specific predictions emerge both before and after eliminating all cosine rhythmicity. Simulations of joint LH and Te time series yielded LH and Te ApEn and LH-Te cross-ApEn histograms with median values of 1.2, 1.8, and 2.1, respectively. These values in the young adult male served as baseline model parameters for further comparisons with simulated disturbances of H functions in various aging models (*experiments 1–7,* below).

#### Heightened Te feedback in aging (experiment 1).

As reported experimentally by Winters et al. (34) and Deslypere et al. (3) parenteral delivery of androgen seems to impose greater feedback inhibition of LH release in aging than in young men. To simulate the impact of this inference mathematically, we left-shifted the dose-response relationship of Te feedback on both GnRH and LH (Fig. 2
*A*). The resultant histograms of uni- and bivariate ApEn of LH and Te concentrations and of the fractional (24 h) Te cosine amplitude are given in Fig. 6. Te ApEn rose, whereas Te concentration and Te cosine rhythmicity fell, consistent with clinical observations in older men. However, LH ApEn, LH-Te cross-ApEn, and the mean 24-h LH concentration fell, which would be inappropriate for known aging data (13).

#### Attenuated LH feedforward on Te (experiment 2).

Human experiments using injected human chorionic gonadotropin as a surrogate for an LH stimulus have revealed a significant decline in maximally stimulated Te production in older men (26, 30,31). Other clinical studies concur with this notion (12). To simulate this putative physiological decrease in testis responsiveness to LH, we reduced the maximum of the LH-Te dose-response curve, thus limiting LH actions on Leydig cells (Fig.2
*A*). This paradigm recapitulated all six observed changes in aging: higher uni- and bivariate ApEn for LH, Te, and LH-Te; reduced mean Te and elevated mean LH concentrations; and, lower 24-hour Te rhythmicity (Fig. 5 and 6).

#### Heightened GnRH feedforward on LH (experiment 3).

A recent detailed intravenous GnRH dose-responsive study of LH release in older versus young men revealed greater maximal LH secretion after GnRH in older individuals (37). To represent this possible primary change in the aging male axis, the maximum of the GnRH-LH dose-response function was increased (Fig. 2
*B*). This simulation mimicked expected age-related changes in LH-Te only by way of elevated 24-h serum LH concentrations (Fig. 6).

#### Accelerated GnRH pulse frequency and decreased GnRH pulse mass (experiment 4).

Deconvolution analysis of intensively sampled (every 2.5 to 10 min) serum LH concentrations in young and older men has disclosed reduced LH secretory burst amplitude, possibly associated with a rise in LH pulse frequency (11, 26, 31). The latter conjecture of accelerated GnRH pulse frequency as a primary mechanism in aging was evaluated by augmenting the mean frequency (intensity) of the point process defining GnRH pulsatility (reduced mean waiting times) (Fig.2
*B*). All three ApEn measures increased (univariate LH and Te, and bivariate LH-Te), but serum LH and Te concentrations did not rise and fall, respectively, as recognized in older individuals.

#### Combined increased Te feedback (experiment 1) and decreased LH feedforward (experiments 2 and 5).

The results for *experiment 5* (combined *experiments 1–2;* see Fig. 2
*B*) are shown in Figs. 5 and 6. This model predicted only that ApEn of LH would rise and that concentrations of Te would decline, as recognized in aging men.

## APPLICATION TO CLINICAL DATA

We next evaluated the impact of detrending of the LH and Te time series using clinical data from healthy older and young men (described in Refs. 11, 12, 18). As shown in Fig.7, both the cosine and heat equation detrending procedures lowered the SDs of the LH and Te time series, thus reducing r [e.g., 0.2SD(Y)] in *Eq. 12
* and increasing ApEn, preserving an age contrast. The latter was also prominent in analysis of the untransformed (nondetrended) data (Fig.7).

## NEW HYPOTHESES SUGGESTED BY COMPUTER-ASSISTED SIMULATIONS

Earlier clinical experiments in young men have disclosed that LH ApEn rises consistently on pharmacological withdrawal of Te negative feedback, e.g., via treatment with the antiandrogenic drug, flutamide (25), or the steroidogenesis inhibitor, ketoconazole (29). Consequently, we also tested the network-predicted response to simulated attenuation of androgen negative feedback on GnRH and LH (*experiment 6*). Such a postulated age-related reduction in Te negative feedback (right shifted response curves) yielded an anomalous rise in serum Te concentrations, which is not evident in older individuals, but predicted appropriate increases in ApEn of LH and Te individually and jointly, as well as a rise in LH concentrations (Figs. 8, *A* and*B*).

The foregoing analysis was also combined with *experiment 2*(reduced LH feedforward on Te) to test the hypothesis of dually impaired Te feedback and LH feedforward [see Fig. 8, *A* and*B* (*experiment 7*)]. This bipartite postulate forecasted increased LH and/or Te ApEn and cross-ApEn values and reduced Te and elevated LH concentrations, as reported in older men (see introduction).

## DISCUSSION

To test specific hypotheses of interglandular pathway disruption within the aging male GnRH-LH-Te axis, the present experiments implement an earlier stochastic differential equation-based biomathematical model that embodies the dynamic time-delayed and dose-responsive feedback and feedforward linkages known within this interactive system (6, 7, 26). We here apply this more formal computer-assisted construct to explore selected hypotheses of altered regulatory behavior due to aging. Several primary hypotheses for more disorderly LH and Te secretion patterns in healthy older men originated from published clinical studies (1, 3, 8, 11, 12, 16,18, 20, 25, 26, 31, 32, 34). In addition, new (secondary) hypotheses emerged on the basis of initial predictions of the ensemble model. Available clinical literature points to multiple plausible pathway (interglandular) or nodal (glandular)-level disturbances in the aging male; e.g., hypothalamic GnRH release, reduced GnRH feedforward on LH secretion, impaired LH feedforward on Te production, and altered Te feedback on GnRH and LH output. However, most earlier experimental investigations of necessity have addressed individual nodal alterations in isolation; i.e., by analyzing GnRH's dose-responsive stimulation of pituitary LH release (37), hCG's feedforward drive of gonadal Te secretion (31), and Te's feedback restraint of GnRH/LH production (22, 24). The present networklike analyses show that single-node studies cannot readily distinguish between partial failure of a control node and disruption of internodal signaling (16). Indeed, a crucial prediction of the present computer-assisted experiments is that disturbing the function of any one regulatory locus disrupts the networklike behavior of this axis, due to the strong feedback and feedforward dependencies that interlink GnRH, LH, and Te secretion (6, 7, 26). Because dynamic interfaces in neuroendocrine axes are multiple, nonlinear, time-lagged, and stochastically variable, a priori intuitive predictions of the impact of nodal disruption and/or pathway failure on overall system performance are very difficult. Indeed, the current formalism highlights the unexpected implications of some earlier intuitions (below).

As a sensitive new statistical tool to quantify the network-level consequences of selected disruption of any given regulatory site and/or its interconnections, we applied the univariate and bivariate (joint) regularity measures, ApEn and cross-ApEn (12, 14-18,27). This family of measures appears to capture among the earliest quantifiable changes in the pattern orderliness of reproductive hormone secretion in the healthy aging male (16) and female (18). As complementary outcome measures, we also used mean and 24-h rhythmic LH and Te output to monitor aging hypothesis-specific predictions (26). Among the five primary (single node) hypotheses explored, we observed that only partial failure of LH's feedforward drive of Te secretion forecasts the principal changes in LH and Te reported in aging men, namely: *1*) elevated LH concentrations, *2*) more disorderly LH and Te release, considered both individually and jointly,*3*) reduced Te concentrations, and *4*) blunted 24-hr rhythmicity of Te secretion. Accordingly, from a feedback/feedforward model perspective, muted feedforward drive by LH of Te secretion appears to be necessary and sufficient to predict the multifold GnRH-LH-Te disturbances documented consistently in clinical experiments (26).

In complementary computer-assisted simulations, we imposed accentuated versus reduced Te negative feedback on GnRH-LH release. These putative pathway defects recapitulated some, but not all, the expected age-related changes in LH and Te secretion recognized in the aging male. In particular, model-based enhancement of Te's negative feedback on GnRH and LH secretion, as suggested by reports of heightened suppression of LH production by pharmacological delivery of androgens in older men (3, 33, 34), forecasted a rise in Te ApEn and fall in Te concentrations, but not an increase in LH ApEn and LH-Te cross-ApEn, as required by clinical data (16). Conversely, computer-simulated attenuation of Te's negative feedback on GnRH and LH secretion, as intimated in other clinical studies (26), predicted higher LH concentrations and LH and Te ApEn and LH-Te cross-ApEn, but an anomalous rise in serum Te concentrations. Thus the foregoing computer-assisted hypothesis testing argues against excessive or deficient Te feedback on GnRH/LH as an exclusive mechanism driving the constellation of reported neuroendocrine changes in the human male reproductive axis in aging.

Unlike the foregoing hypothesis of an isolated reduction in Te's negative feedback on GnRH-LH output in the aging male, a bipartite hypothesis, which combined the latter defect with impaired LH feedforward drive of Te secretion, predicted the expected (older male) neuroendocrine phenotype of elevated ApEn of LH and Te, higher cross-ApEn of LH-Te, increased LH and reduced Te secretion, and blunted 24-h rhythmic Te production. Accordingly, these biostatistical experiments allow for, but do not prove, a more complex gonadal-axis pathophysiology in older men, consisting of dual defects that involve both positive and negative control; i.e.: *1*) an impairment of LH's feedforward on Te, and *2*) a reduction in Te's feedback on GnRH and LH. To affirm or refute this bipartite prediction will require further clinical and laboratory experiments.

Further analyses were carried out to evaluate the network implications of heightened endogenous GnRH feedforward on LH in the aging male. This clinical hypothesis reflects the empirical finding that graded GnRH injections stimulate greater LH secretion in older men than young men (37). However, computer-simulated feedforward enhancement of GnRH on LH, although emulating increased LH output, failed to drive elevated ApEn or cross-ApEn of LH and Te or lower Te concentrations. According to these data, we posit that the reported heightening of GnRH potency and/or efficacy in older men may be a consequence (rather than a cause) of diminished Te feedback. Indeed, reduced Te feedback is inferrable from the well established 30–50% fall in bioavailable (nonsex hormone binding globulin associated) testosterone in older individuals (9, 10, 13, 26, 31, 33).

Accelerated LH pulse frequency and reciprocally reduced LH pulse amplitude have been reported in some aging men (3, 11, 26, 28,31). We reasoned that both neuroendocrine features might be secondary consequences of restricted Te availability and hence diminished negative feedback on GnRH/LH. In particular, on the basis of the present biomathematical systems model, restraining Te negative feedback (by right shifting the corresponding negative dose response functions of Te on GnRH and LH secretion) projected an increase in GnRH/LH pulse frequency and a reciprocal decrease in LH peak amplitude. The latter inverse relationship is also consistent with reportedly damped gonadotrope responsiveness to higher GnRH pulse frequencies in experimental animals (2). Thus the foregoing model-predicted considerations highlight the thesis that examining any one control node (e.g., here GnRH's feedforward action on LH) in isolation may not illuminate full network operation (see introduction).

The bipartite postulate of age-related defects in LH drive of Te secretion along with augmented Te feedback on GnRH and LH output nullified three of the age-appropriate outcomes otherwise predicted by the single (uncombined) model of decreased LH feedforward on Te. In contrast, combining impaired LH stimulation of Te secretion with reduced Te feedback on GnRH/LH predicted relevant age-related entropic, mean and 24-h rhythmic changes in LH and Te output. These feedback distinctions, if validated independently by further in vivo physiological experiments, illustrate the challenge of intuitively based predictions and the utility of more formalized biomathematical tools to explore combinatorial defects within an axis.

The present study did not attempt to appraise all possible paired (or multiply combined) mechanisms of pathway disruption in the aging male. Thus other plurinodal or multipathway models might also explicate reported disturbances in the older male reproductive axis. Further studies will be needed to explore this interesting issue. In addition, it will be important to investigate further model sensitivity to different renditions of circadian linkages (e.g. imparting diurnal periodicity to GnRH pulse generation and/or GnRH's action on gonadotropes), of the interface functions for LH-Te and Te-GnRH-LH, of time-lagged and time-averaged inputs, and stochastic variability embodied in feedback and/or feedforward interfaces. The current model-based perspective on hypothesis testing will be important to explore further in relation to the more complex menstrual and aging-related dynamics of the female reproductive axis, as well as the integrative control of the thyrotropic, lactotropic, somatotropic, and corticotropic axes.

### Perspectives

The time-evolving behavior of an integrated neuroendocrine ensemble is difficult to envision intuitively, because of stochastic features, time delays, joint feedback and feedforward activities, and nonlinear dose dependencies mediating the several interactions. Analytically tractable computer-assisted biomathematical models seem critical to amalgamate the consequences of such complex linkages conceptually. Here, we illustrate an application of one published model structure of the male GnRH-LH-Te axis to specific hypothesis testing in aging. In physiological applications of this kind, the biostatistical construct should enhance understanding of and insights into potentially confounding, costly, invasive, and/or extended experimental measurements of observable output (e.g., LH and testosterone). As importantly, the model should ultimately aid in the accurate prediction of unobservable output (e.g., GnRH) and internodal connections (e.g. Te's feedback on GnRH LH). Assumptions inherent in the primary model structure must be made explicit, remain faithful to the extant body of knowledge, and be tractable to amendment by evolving scientific data. In the last regard, ad hoc elements of the model thereby will be refined (e.g., the stochastic properties of the pulse generator system) and fundamental connections refuted or corroborated (e.g., the putative sinusoidal modulation of the LH Te feedforward pathway). Accordingly, structural model progress is coupled to experimental developments and vice versa. Robust model forms should therein foster more insightful data reconstruction, enhance experimental strategies, and stimulate new mechanistic hypotheses.

## Acknowledgments

Supported in part by the University of Virginia Center for Biomathematical Technology, National Institutes of Health (NIH) General Clinical Research Center Grant M01-RR-00847, NIH RCDA 1K04-HD-00634, and NIH U54 Center for Specialized Studies in Reproduction P30-HD-28934.

## Appendix

### Removing the Circadian Rhythm via the Heat Equation

The heat (or diffusion) equation is one of a triad of fundamental partial differential equations in physical mathematics. Historically, it was used to describe the time evolution of such spatially defined processes as the temperature distribution in a heat-conducting medium and dispersion of a chemical concentration. Today, however, it is also being used in such areas as image processing for noise removal and image restoration. Within the later context, it is typically used as a low-pass filter. However, as described below, other special properties make it appropriate to the present intent of removing the effects of the circadian rhythm. To this end, we construct its complementary high-pass filter, which is analogous to image deblurring or enhancement.

For the present problem, this procedure functions essentially as Gaussian smoothing, wherein the degree of smoothing is not strictly fixed. Rather, the oscillatory nature of the data influences the amount of smoothing. One must slightly qualify the preceding expectation, because there is the issue of what to do at the boundary of the data series. For the present description, we illustrate that the algorithm handles the boundary in an intelligent manner.

Let F(r) denote a signal at time r, and let μ(r) be a periodic function, which modulates the signal F. In the present construct, F is the feedforward signal of LH on testosterone, which is a time-delayed and time-averaging of LH concentrations. Because of the pulsatile nature of LH secretion, F will tend to be composed of higher than circadian frequencies with periods on the order of 1–3 h, along with even higher frequency sample-by-sample variability. The circadian rhythm μ, on the other hand, primarily represents low frequencies with periods on the order of 24 h. In methods, the rate of secretion for testosterone was defined, at *time*r, as a sum of the basal and nonbasal components: β + H[μ(r)F(r)], where H is a dose-response function. The resulting Te concentration is then given by
Equation 13where α is the rate of elimination. If one ignores the first term due to the initial Te concentration, then X^{(μ)} is a convolution integral. Let X represent the solution where μ ≡ 1, so there is no circadian rhythm. The goal is to recover X approximately, given only observed X^{(μ)} and without making any rigid assumptions concerning μ. We thus have
Equation 14
Equation 15where E is the one-sided exponential (e^{−αx}, x ≥ 0; 0 for x < 0), and * denotes convolution.

Removing the circadian input requires eliminating the effects of μ from the right-hand side of *Eq. 14,* without altering F in the process. The output of the heat equation algorithm after *time *a is essentially a convolution of the input X^{(μ)} with the Gaussian density, G_{a}, of mean 0 and variance a^{2}: G_{a} ∗ X^{(μ)}. The foregoing is an appropriate filter in the present setting, because it is a convolution, and, more importantly, because convolving G_{a}with a cosine function of frequency λ, merely multiplies the cosine function by e^{−(λ2}
^{a2}
^{/2)}.

For the sake of exposition, let μ and F each consist of a single frequency (low and high, respectively): μ(r) = 1 + B_{1}cos (λ_{0}r + φ_{0}), F(r) = D_{0} + D_{1}cos (λ_{1}r + φ_{1}). If μ and F are respectively a sum of low-frequency terms and a sum of high-frequency terms, the result will follow by linearity. What is key is that when a^{2} is not terribly large, low-frequency sinusoidals remain virtually unaffected, whereas high-frequency sinusoidals are dampened. The properties of linearity, commutivity, and associativity of convolutions are also important, as well as the fact that since G_{a} integrates to one, its convolution with a constant results in that constant: G_{a} commutes with E, E and G_{a} each commute with the operation of averaging, and G_{a} ∗ β = β.

Because μF equals D_{0}μ + F − D_{0} plus any terms of reduced amplitudes at frequencies λ_{1} ± λ_{0} (which are high frequencies), then G_{a} ∗ μF ≈ D_{0}μ. Also, because μF = D_{0} + a sum of cosines at λ_{0}, λ_{1}, λ_{1} ± λ_{0}, one can reasonably assume that the average of μF: μ̅F̅ ≈ D_{0}. In the region of the dose-response function H that is approximately linear
and hence
If the domain of the spatial variable (t) were unbounded (−∞ < t < +∞), then the solution of the heat equation U(t, s) in *Mean function estimation via smoothing by the heat (or diffusion) equation* is given by the convolution of X^{(μ)} with G_{a}, where X^{(μ)} = Y (the data) and a^{2} equals C^{2}s. The variable t is observational time, and s is (smoothing or) algorithmic time. If the spatial domain is bounded, as in most applications, the solution is slightly different, but one can still solve the heat equation numerically by way of the algorithm described in methods.The remaining issue is how long to implement the algorithm (i.e., what is the choice of a) in the present context? Here, because μ and F will tend to have frequencies at the two extremes, the choice of a is less critical. In *Mean function estimation via smoothing by the heat (or diffusion) equation* the algorithm is implemented for an amount of algorithmic *time*S (i.e., a^{2} above equals C^{2}S). Certainly if S were taken to be enormous, the algorithm would smooth the input to a constant. In the present application, where there are few midrange frequencies of interest for ApEn analysis, the choice does not require subtlety.

## Footnotes

Address for reprint requests and other correspondence: J. D. Veldhuis, PO Box 800202, Univ. of VA, Charlottesville, VA 22908-0202 (E-mail: jdv{at}virginia.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