INTRODUCTION

Top Abstract Introduction Methods Discussion References

THE MALE REPRODUCTIVE AXIS consists of three physiologically distinct but interacting functional control nodes: a gonadotropin-releasing hormone (GnRH) pulse generator, endowed by hypothalamic neurons; luteinizing hormone (LH), produced in the anterior pituitary gland; and testosterone (Te), secreted by Leydig cells in the testes. In health, this multinodal feedback and feedforward interactive system yields a pseudo-steady-state output of pulsatile (episodic) neurohormone release that shows circadian modulation, e.g., 24-h variations in serum Te concentrations and, to a lesser extent, in LH (30). Dose-response relationships have been largely defined for individual nodes acting in isolation, e.g., GnRH's stimulation of LH secretion (1) and LH's dose-dependent stimulation of Te secretion (10). Similarly, earlier simulation modeling of neurohormone release has typically encompassed a single hormone, not the entire interacting feedback system. A major limitation inherent in thus isolating system components that are so highly coupled physiologically via time-lagged feedforward (e.g., LH-Te) and feedback (e.g., Te-LH) mechanisms is omission of the influences due to communications among the component(s). Moreover, artificial isolation of functional elements can make inference of system behavior more difficult. Given these issues, we present an initial biomathematical model of an entire interconnected three-nodal system, i.e., the (male) hypothalamic (GnRH)-pituitary (LH)-gonadal (Te) axis, and test its basal pulsatile output, modulated circadian responses, and predicted performance in selected simulations and prior human experiments.

Glossary

i	Elimination rate constant for hormone i
i	Basal secretion rate of hormone i
t	Discretization step size used in computer simulations
dS(t)	Incremental secretion in interval (t, t + dt) (from Ref. 16)
dWi(t)	Stochastic noise term (via differential of Brownian motion) (see section VIII)
dXi(t)	Incremental change in concentration of hormone i at time t
Zi(t)dt	Incremental change in secretion of hormone i at time t
	Parameter that (probabilistically) controls interpulse lengths
GnRH	Gonadotropin-releasing hormone
H( · )	Dose-response (interface) functions
IID	Independent and identically distributed
(lj,1, lj,2)	Time-delayed interval for jth feedback interaction
( · )	Cosine function specifying periodic (circadian) input
( · )	(Stochastic) pulse generator intensity function
LH	Luteinizing hormone
Mji	Pulse mass j for hormone i
p(t)	Pulse generator (stochastic) intensity
i( · )	Pulse shape for hormone i
SDE	Stochastic differential equation
T ij	Pulse time j for hormone i
Te	Testosterone
Wi( · )	ith Brownian motion process (one of "noise" processes) (see section VIII)
Xi(t)	Concentration of hormone i at time t
Yik	Observed concentration of hormone i at time tk
Zi(t)	Secretory rate for hormone i at time t

	I. GENERAL METHODS

A. Background Physiology

Our thesis is that the complex dynamic of the male hypothalamic-pituitary-gonadal feedback system occurs because the rate of secretion of any given hormone within the system (GnRH, LH, and Te) depends on relevant time-delayed, nonlinear feedback and feedforward signals derived from and acting on all or some of the components of the system. We further assume that pertinent dose-responsive interfaces (either inhibitory or stimulatory) connect hormone signal input to nodal output; e.g., a GnRH concentration-LH secretion rate dose-response relationship functionally links the hypothalamic GnRH signal and the time-delayed pituitary LH secretory output (8, 9, 11-14, 25, 30, 35). Similarly, a dose-response relationship exists for LH concentrations and Leydig cell Te secretion, as established by in vitro and in vivo experiments (2, 7, 10, 30). By formulating a physiologically linked network using the three primary and interacting nodes of hypothalamus, pituitary gland, and testes (Fig. 1), we can examine basal hormone output, test system performance in specific computer simulations compared with prior clinical experiments, evaluate relevant internodal linkages, and later predict possible axis dysregulation.

View larger version (16K): [in this window] [in a new window]   Fig. 1.   Schematic illustration of time-delayed negative feedback () and positive feedforward (+) within human male gonadotropin-releasing hormone-luteinizing hormone-testosterone (GnRH-LH-Te) axis. Broad arrows, feedforward (+) stimulus-secretion linkages; narrow arrows, feedback () inhibition. "H" functions are developed further in section I and Fig. 2 and define dose-response relationships at each feedback interface within axis (see section VIIIA2).

1. Core equations. Considering the foregoing background, we can relate the "true blood hormone (GnRH, LH, or Te) concentrations" in vivo to corresponding secretion rates. To this end, let time 0 represent the onset of the observation period, X_G(t), X_L(t), and X_Te(t) (t  0) the hormone concentration values, and Z_G(t), Z_L(t), and Z_Te(t) (t  0) the corresponding hormone secretion rates for GnRH, LH, and Te, respectively, over time. The core model of a dynamic GnRH-LH-Te feedback system will then encompass the following general equations. These equations state that, by definition, for each hormone signal involved (GnRH, LH, or Te), the rate of change of its concentration in blood is the difference between its rates of elimination and production

(1)

where X_G(0), X_L(0), and X_Te(0) are specified initial GnRH, LH, and Te concentrations and _G, _L, and _Te are the rates of elimination of the respective hormones; the rates could be allowed to be concentration dependent, as may be the case for higher levels of LH (32). Concentration (and secretion rate) units are mass of neurohormone (secreted) per unit distribution volume (per unit time).

	II. FEEDBACK BEHAVIOR ANTICIPATED IN THE MALE GONADOTROPIN-RELEASING HORMONE-LUTEINIZING HORMONE-TESTOSTERONE AXIS

Pulses of GnRH from the hypothalamus drive corresponding episodes of pituitary LH secretion in a nearly uniformly 1:1 ratio (1, 8, 9, 14, 25, 27, 30, 35). Pulsatile LH concentrations bathing Leydig cells in the testes in turn stimulate Te secretion in a dose-dependent manner (7, 10). Thus serum Te and, to a lesser extent, serum LH concentrations vary episodically over 24 h. In addition, available data are consistent with circadian variations, with maximal hormone concentrations of LH and Te occurring in the human during the later portion of nighttime sleep (34). In general, the amplitude of LH pulses tend to vary inversely with event frequency and to be maximal at night (2, 9, 24, 25, 33). Despite these variations, mean daily serum LH and Te concentrations remain within a relatively narrow (0.5- to 1.5-fold) physiological range. This is thought to reflect homeostatic feedback control, which we embody in the coupled equations above (Eq. 1). More explicitly, in young men, serum LH and Te concentrations are positively cross-correlated at a +20- to 50-min lag (Te following LH) and negatively cross-correlated at a 80- to 100-min lag (testosterone preceding LH) (34). The former reflects LH's dose-dependent feedforward action on Leydig cell Te biosynthesis. The latter presumably reflects Te's negative feedback on GnRH-LH secretion (30). Reduction or removal of Te's negative-feedback signal via an androgen-receptor antagonist or an inhibitor of Leydig cell steroidogenesis increases the frequency and amplitude of pulsatile LH release as well as, in the former case, the mean serum Te concentration (17, 29, 36). Conversely, continuous intravenous infusion of Te in steroidogenesis-inhibited men suppresses pulsatile LH release by reduced LH (and presumptively GnRH) pulse frequency with escape of occasional higher-amplitude LH pulses (36). These published clinical data provide a basis for testing the dynamics of our feedback model of the GnRH-LH-Te axis (see below).

	III. SPECIFIC METHODS

Top Abstract Introduction Methods Discussion References

A. Constructing GnRH-LH-Te Secretion

To represent explicitly the changing rates of pulsatile and continuous secretion over time, we assume that GnRH of hypothalamic origin signals the gonadotropic cell to increase its release and synthesis of LH molecules. The GnRH stimulus acts over a finite time to increase the rate of de novo LH biosynthesis while concurrently prompting the nearly immediate release of (readily releasable) membrane-associated, storage granule-encapsulated LH. Newly synthesized LH molecules are encapsulated into such secretory granules, which diffuse out toward the gonadotropic cell membrane as facilitated by cytoskeletal elements to direct vectorial drift (30). When storage granules containing LH reach the cell membrane, they can undergo secretory exocytosis (when the GnRH signal is activated) or remain marginated (during secretory quiescence) and thereby accumulate, awaiting the next GnRH signal to trigger release.

For nonprotein hormones (unlike LH), there is no known granule accumulation process but, rather, more nearly continuous release, as presumed here for the steroid hormone Te (7, 10, 30). In contrast, the foregoing granule storage-release process applies to LH and GnRH, for which accumulated intracellular secretory granules become available for rapid initial release at the onset of the next signal, resulting in pulsatile release. The proximate signal for such a pulse of LH release is GnRH, and for GnRH release it is a balance of stimulatory and inhibitory neurotransmitter inputs.

In contrast to pulsatile LH release by gonadotropic cells, some LH molecules are not encapsulated or are encapsulated but not hormonally regulated and diffuse out toward and through the plasma membrane, thus contributing the so-called basal secretion described above (30). Because the pituitary gland contains an array of individual gonadotropic cells, integrating LH secretion over all constituent cells yields the overall LH secretion rate compounded from basal and pulsatile release. Similar integration of secretion over the ensemble of dispersed hypothalamic GnRH neurons is pertinent.

1. Continuous secretion. For continuously released Te (30), the secretion rate at time t will be assumed to depend on X_L( · ), the concentration levels of the input stimulus, LH, over some time-delayed interval, e.g., from t  l₂ to t  l₁, as transduced via some dose-response function H₄( · )

(2)

where _Te is the Te basal release rate. Thus we model Te release via the function H₄( · ) above; the precise assumed form of H₄( · ) is given in sections IVA and VIII.

2. Pulsatile secretion. Pulsatile secretion (of GnRH and LH) is marked by pulse times, when the rate of hormonal secretion rapidly increases, followed by a possibly not so rapid decrease; this secretory episode will be called a pulse (30). Here we allow one principal stochastic pulse generator, i.e., hypothalamic GnRH, from which pituitary LH inherits its episodicity of release (30). Thus, for the GnRH-LH-Te axis, there will be one set of pulse times, T₀, T₁, T₂, ... , where without loss of generality we will assume T₀ = 0.

j  1

(3)

B. Circadian Modulation

We previously formulated a model of the pulsatile secretion (not concentration) of one hormone, without external feedback inputs, as illustrated by episodic LH secretion (16). Pulsing was described by an inhomogeneous Poisson process, for which the deterministic intensity ( · ) function was 24-h periodic, thus describing a circadian rhythm. In the simulations, the ( · ) function was assumed to consist of one harmonic, with amplitudes B₀ and B₁ and phase shift ₁ being chosen so ( · ) would vary between a high of 1 at 4 AM and a low of 0.6 at 4 PM (B₀ = 0.85, B₁ = 0.15, ₁ = 4)

The mass of LH available for release within a pulse, M^j, was assumed to be a linear function of the preceding interpulse length: the longer the interpulse interval, the greater the accumulation of LH-containing granules available for release where M^j represents the pulse masses, ₀ denotes a minimum amount of mass always secreted per pulse, and ₁ is a constant rate of hormone accumulation. Incremental secretion, dS(t) (from the pituitary gland of a horse), was evaluated previously (16). If S(t) is the cumulative secretion up to and including time t, then the incremental secretion between sample times t_i and t_i + 1[S(t_i + 1)  S (t_i)] can be obtained by integrating the differential dS( · ); we previously (16) formulated dS(t) as The Z(t) in our present formulation corresponds to dS(t)/dt in this earlier construction.

Thus the only feedback in the earlier formulation (16) was through the length of the previous interpulse interval, but the instantaneous rate of accumulation, ₁, was constant and, hence, unaffected by feedback; also the point process representing the pulse times was modulated by a deterministic (periodic) intensity function ( · ), but not by feedback. Whereas this preliminary model is very reasonable and appears to fit certain secretory data quite well, occasionally there may be very short interpulse intervals followed by relatively large pulse masses, and vice versa (23). We reasoned that this pulse-to-pulse variability might be explicable if one models the entire system. In addition, by incorporating physiologically pertinent feedback and feedforward connections, one should derive more meaningful estimates of the true variability within the intact system than by modeling release of a single hormone in isolation (6).

	IV. COMPOSITE BIOMATHEMATICAL CONSTRUCT OF THE MALE GONADOTROPIN-RELEASING HORMONE-LUTEINIZING HORMONE-TESTOSTERONE AXIS

A. Dose-Response Linkages [H( · ) Functions]

Experimental evidence suggests that internodal feedback (e.g., GnRH feedforward on LH, LH feedforward on Te, Te feedback on GnRH or LH) is exerted via a time-delayed time average of the effector concentration. Thus, throughout the simulation, we will use the following notation, where 0  t₁  t₂ < , to denote feedback signal intensity (i.e., denotes a time average) Using empirically inferred dose-response (possibly multivariate) logistic functions, which have been evaluated in in vitro and in vivo experiments in humans and animals (1, 30, 35), we then define how such feedback interacts with the various components

(4)

If B_i > 0, the feedback is positive (i.e., feedforward effect); if B_i < 0, the feedback is negative. Accordingly, for the simplified male GnRH-LH-Te axis, there are four major relevant feedback-feedforward dose-response functions: H₁( · ) for the GnRH pulse firing rate as a function of Te concentration, H₂( · ) for the rate of GnRH pulse-mass accumulation as a function of Te concentration, H₄( · ) for the rate of Te secretion as a function of LH concentration, and H_5,6( · ,  · ) for the rate of LH pulse-mass accumulation as a function of Te and GnRH concentrations. The subscripts correspond to feedback-feedforward relationships (1-7 in section VIII). A function H₇( · ), not modeled by a logistic form, describes a refractory condition of the GnRH pulse generator as a consequence of possible short-loop negative autofeedback of GnRH on its own release (see section VIII); also, an optional interaction given by H₃( · ), not part of the basic construct, allows for the Te basal secretion rate to vary with a 24-h periodicity. Feedback connections are thus mediated via relevant dose-response functions H₁( · ), H₂( · ), H₄( · ), and H_5,6( · ,  · ), wherein the maxima, minima, slopes, and midpoints reflect the operating behavior of GnRH-LH, LH-Te, Te-GnRH, and Te-LH linkages. Also, the time delays for the jth feedback interaction will be expressed throughout by l_j,1 and l_j,2. For example, the negative feedback of blood Te concentration (in ng/dl) on the rate of hypothalamic GnRH pulse-mass accumulation (in pg · ml¹ · h¹) would be of the form

(5)

Also, until time t is above the maximum time delay, the feedback will not be over the full time-delay interval, but rather only the amount of time since time 0 (as discussed in section VIII). Figure 1 displays the above feedback connections, and Figure 2 illustrates the corresponding mathematical functions.

View larger version (44K): [in this window] [in a new window]   Fig. 2.   Dose-response feedback and feedforward functions for male GnRH-LH-Te axis. First (leftmost) vertical column displays, for a normal subject, dose-response feedback and feedforward functions of Fig. 1: Te inhibition of GnRH pulse firing rate [H1( · )], Te inhibition of GnRH pulse-mass accumulation rate [H2( · )], LH stimulation of Te secretion rate [H4( · )], and Te and GnRH jointly acting on LH pulse-mass accumulation rate [H5,6( · ,  · )]. Second column depicts H1( · ) and H2( · ) for simulation 1 (consisting of flutamide treatment), H4( · ) for simulation 2 (ketoconazole treatment), and H5,6( · ,  · ) for simulation 1. Third column presents assumed dependency of slope (B1) as a function of Te in H1( · ) for simulation 3 (ketoconazole treatment with Te infusion add back); 2 pulse shapes, G( · ) (solid line) and L( · ) (dashed line); 24-h periodic modulating function ( · ) with its 12-h phase-shifted version superimposed (dashed line); and assumed dependency of lower asymptope (D) as a function of Te in H5,6( · ,  · ) for simulation 3. [Dose-response H( · ) functions are defined in section VIIIA2, and simulations are summarized in sections VB1 and VB4.]

Figure 2 depicts the various dose-response feedback functions for the simulation experiments; the nominal parameter values for the feedback functions and elimination rates are given in section VIIIA along with their motivation. In the sensitivity analysis (section VC), these nominal values are doubled and halved to demonstrate the dependency of the structural dynamics on the parameter values. The pulse shapes _G( · ) and _L( · ) were taken to be generalized  densities as displayed in Fig. 2 (2nd row, 3rd column); the modulating function ( · ) (and its 12-h phase shift) used in the circadian mechanistic simulations is also displayed in Fig. 2 (3rd row, 3rd column). The variances for the combined technical and measurement errors were such that the coefficient of variation was 6%. The parameter  in the construction of the pulse generator (see Eq. 7) was taken to be 2 (see section VIIIB1).

B. Basic SDE Construct

(6)

k  1

(7)

(8)

(9)

(10)

What we then actually observe is a discrete-time sampling of these with attendant technical/measurement error due to sample collection, processing, and assaying In the simulations, the _i values are independent and identically distributed (IID) normal, mean zero and with variances such that the coefficients of variation for Y^G_k, Y^L_k, and Y^Te_k are 6% to fall within an anticipated operating range of 3-15% for typical available assay coefficients of variation (29-36).

C. Mechanisms for Incorporating a Circadian Rhythm

Section VIII gives a further description of the first six possible models for the circadian rhythms in Te and LH. We test (see section V) the predictions of these six putative circadian coupling schemes and show that some, but not all, predict physiological variations in 24-h serum Te (more than LH) concentrations. Mechanism 7 (GnRH autofeedback) is less likely, we speculate, which would entail day-night variations in GnRH feedback inhibition of its own secretion.

	V. COMPUTER-ASSISTED SIMULATIONS OF THE MALE GONADOTROPIN-RELEASING HORMONE-LUTEINIZING HORMONE-TESTOSTERONE AXIS

Simulations from the above discretization of the biomathematical formulation of the GnRH-LH-Te time-delayed feedback axis were implemented (using Matlab). The sampling interval was 30 s, and the simulations were extended for 48 h. The second 24 h were recorded, thus removing any startup effects. The six circadian mechanistic linkages were simulated, as were data from four previously published clinical experiments (see section VB, 1-4). For each simulation, 500 realizations were accomplished. We display the first realization along with the average of the 500 realizations superimposed in Figs. 3-6, which depict pulsatile LH, Te, and GnRH concentrations over a 24-h period for the six circadian mechanisms and simulations of earlier clinical experiments 1-4. For circadian mechanism 4, the secretion rates of LH, Te, and GnRH are also displayed (Fig. 3, bottom).

View larger version (44K): [in this window] [in a new window]   Fig. 3.   Feedback modulation of GnRH pulse firing rates and pulsatile-circadian coupling mechanism 4. Top: 9 panels in the 3 horizontal rows illustrate feedback modulation of GnRH pulse firing rates for 6 circadian mechanisms (1-3 in row 1 and 4-6 in row 2; see sections IVC and VIIIC) and for feedback simulations 1-3 (row 3; see sections VB1 and VB4). Bottom: mechanism 4, i.e., 24-h modulation of feedforward action of LH on rate of Te secretion. Simulations of concentration levels (left) and secretion rates (right) for LH (top), Te (middle), and GnRH (bottom) are shown over 24 h, discretized to every 30 s, sampled every 10 min. Each plot consists of 1 realization with average of 500 realizations superimposed.

View larger version (43K): [in this window] [in a new window]   Fig. 4.   Pulsatile-circadian coupling mechanisms 1, 2, 3, and 5. A: mechanism 1, i.e., periodic 24-h modulation of negative-feedback actions of Te on GnRH pulse generator firing rate. B: mechanism 2, i.e., similar modulation of negative-feedback effects of Te on GnRH secretion rate (mass). C: mechanism 3, i.e., diurnal modulation of basal secretion rate of Te. D: mechanism 5, i.e., nyctohemeral modulation of negative feedback of Te on rate of LH secretion. Simulations of concentration levels for GnRH (top), LH (middle), and Te (bottom) over 24 h, discretized to every 30 s, sampled every 10 min are shown. Each plot consists of 1 realization with average of 500 realizations superimposed.

View larger version (44K): [in this window] [in a new window]   Fig. 5.   Pulsatile-circadian coupling mechanism 6 and feedback simulations of clinical experiments 1-3. A: mechanism 6, i.e., 24-h rhythmic modulation of stimulatory actions of GnRH on rate (mass) of LH secretion. B: simulation 1, i.e., decreased Te negative feedback on LH and GnRH as imposed experimentally via flutamide, an antiandrogen. C: simulation 2, i.e., nearly complete (90%) withdrawal of Te negative feedback achieved via ketoconazole, which blocks testicular steroidogenesis. D: stimulation 3, i.e., nearly complete (90%) withdrawal of Te negative feedback, for 48 h, along with a total daily dose of 8 mg of Te infused continuously over 24 h. Simulations of concentration levels for GnRH (top), LH (middle), and Te (bottom), over 24 h, discretized to every 30 s, sampled every 10 min are shown. Each plot consists of 1 realization with average of 500 realizations superimposed.

View larger version (54K): [in this window] [in a new window]   Fig. 6.   Feedback simulation 4 (45-, 60-, 90-, and 120-min GnRH pulses). In simulation 4, GnRH pulse generator is "clamped" to simulate fixed periodic injections of GnRH every 45 min (A), 60 min (B), 90 min (C), and 120 min (D). Simulations of concentration levels for GnRH (top), LH (middle), and Te (bottom), over 24 h, discretized to every 30 s, sampled every 10 min are shown. Each plot consists of 1 realization with average of 500 realizations superimposed.

A. Implementation of Possible Circadian Rhythm Mechanisms

In the first linkage mechanism considered between pulsatile and circadian rhythms, we test the proposition that Te negative feedback on GnRH secretory burst frequency is modulated by a circadian input, e.g., with greater feedback (inhibition) at night. These simulations elicited relatively little 24-h variation in serum LH and Te concentration profiles (Fig. 4A, top and middle). Simulations of a second mechanism of coupling a circadian oscillator (e.g., speculatively via suprachiasmatic nucleus input) to GnRH neuronal secretory burst mass/rate also yielded relatively small 24-h rhythms for all three primary components of the male axis: GnRH, LH, and Te concentrations (Fig. 4B). In contrast, clinical observations suggest greater Te than LH circadian rhythmicity (34). Third, simulating a construct of 24-h rhythmicity of basal Te secretion (unmodulated by LH) clearly predicts (Fig. 4C, middle) rhythmic serum Te variations over 24 h with lesser changes in GnRH and LH. Fourth, proposed circadian modulation of LH's stimulation of Te secretion, via progressive shifting of the LH-Te dose-response curve sensitivity across 24 h, yields strong (Fig. 3, bottom) day-night variations in LH and Te (but not GnRH) concentrations. This also is consistent with known physiology. Fifth, 24-h modulation of Te's feedback on LH secretory mass yields predictions with predominant rhythms in LH greater than Te concentrations (Fig. 4D). Finally, periodic variation over 24 h of GnRH's feedforward (stimulatory) action on LH release evokes relatively little diurnal rhythmicity in Te concentrations (Fig. 5A). Thus, qualitatively, varying basal Te secretory rates, Te's feedback on LH secretion, or LH's feedforward action on Te via a deterministic periodic (24-h) input will emulate the known (human) physiology of greater Te than LH circadian variation.

B. Simulations of Earlier Clinical Feedback Experiments

1. Simulating decreased Te negative feedback on LH and GnRH. The negative-feedback actions of Te (presumptively on the masses of GnRH and LH secreted and on the GnRH pulse generator frequency) have been antagonized pharmacologically in clinical experiments via the administration of a nonsteroidal antiandrogen, flutamide, which inhibits Te's binding to the androgen receptor (17, 29). This pharmacological treatment would in effect shift the Te negative-feedback dose-response curves in our biomathematical construct to the right. Thus to simulate flutamide's actions, we have introduced an 80% shift (by way of an evident decrease) in feedback sensitivity of GnRH and LH pulse mass and GnRH pulse frequency to Te feedback inhibition (Figs. 1 and 2). As shown in Figs. 5B and 3 (top, 3rd row, 1st column), this simulation disclosed increased serum Te concentrations, with accelerated LH secretory burst frequency and amplified LH secretory burst mass, as reported earlier in young men (17, 29).

2. Simulating nearly complete (90%) withdrawal of Te negative feedback. A low and constant blood concentration of Te has been induced clinically by administration of the drug ketoconazole, which blocks the ability of the testis to synthesize Te (36). To achieve this state in our SDE feedback construct, we fixed Te secretion at 20 ng · dl¹ · h¹ (vs. 1,350 ng · dl¹ · h¹ in nominally normal conditions). The corresponding simulations correctly predicted (Fig. 5C and Fig. 3, 3rd row, 2nd column) the clinical observations reported earlier in young men treated with this inhibitor (36), i.e., an increase in the frequency and mass of LH secreted per burst.

3. Simulating nearly complete (90%) withdrawal of endogenous Te secretion followed by continuous Te (add-back) infusions. In six men treated earlier with the drug ketoconazole for 48 h (see section VB2) to block Te biosynthesis and then intravenously infused with 8 mg of Te continuously over hours 24-48, serum Te concentrations rose on average from 45 to 1,000 ng/dl over the latter interval (36). In this published clinical experiment, virtually all the Te in the system arises from the infusion. Our feedback simulations predicted (Fig. 5D, and Fig. 3, top, 3rd row, 3rd column) the clinically evident pulsatile LH profiles in these young men (36). Here, we assumed that the gradual increase in LH secretion over hours 0-24 of ketoconazole treatment and, hence, Te withdrawal (but before Te add back) tended to reduce the negative-feedback minima for Te on LH secretory burst mass and GnRH pulse frequency, as shown in Fig. 2 (top right and bottom right, respectively).

4. Simulating GnRH pulse generator as "clamped" at various frequencies by fixed periodic injections of GnRH. Simulations consisted of fixed GnRH injections every 45, 60, 90, or 120 min, and the output profiles of serum LH and Te concentrations were recorded (Fig. 6A). Here, GnRH (exogenous) is no longer regulated negatively by itself or by Te negative feedback. This paradigm evinced the key previously reported features of an inverse relationship between LH pulse frequency and burst mass and of a plateauing rise in mean LH concentrations at higher LH pulse frequencies (LH secretion ultimately is inhibited by rising Te negative feedback; Fig. 6), as reported earlier independently in in vivo experiments in the sheep and human (9, 25).

C. Sensitivity (Analysis) of GnRH-LH-Te Output to Changes in Dose-Response Parameter Values

In particular, for feedback-feedforward interactions 1, 2, and 4, each of their corresponding dose-response functions is parameterized by values A, B₁, C, and D; in the case of interactions 5 and 6, the dose-response function was parameterized by A, B₁, B₂, C, and D. The maximum and the minimum were modified by multiplying C and D first by 0.5, then by 2. In the case of interactions 1, 2, and 4, the midpoint was modified by multiplying B₁ first by 0.5, then by 2; in the case of interactions 5 and 6, we modified B₁ by 0.5 and 2 and B₂ by 0.5 and 2. The resulting modified (vs. unmodified) dose-response functions for interactions 1, 2, and 4 are displayed in the first two rows of Fig. 7 and those for interactions 5 and 6 in the third row. Corresponding realizations from our model are shown in Fig. 8, i.e., the resulting concentrations of LH and Te over time.

View larger version (43K): [in this window] [in a new window]   Fig. 7.   Sensitivity analysis: modifications of dose-response functions. Dose-response functions (see Fig. 1) corresponding to feedback-feedforward interactions 1, 2, and 4 (as defined in section VIIIA2) are parameterized by values A, B1, C, and D and interactions 5 and 6 by A, B1, B2, C, and D. For each dose-response function, 2 variations were considered: 1) maximum and minimum of dose-response functions were simultaneously increased, then simultaneously decreased, and 2) midpoint of dose-response function was shifted to left, then to right. Maximum and minimum were modified by multiplying C and D first by 0.5 and then by 2. In case of dose-response (interface) functions enumerated in section VIIIA2 for interactions 1, 2, and 4, midpoint was modified by multiplying B1 first by 0.5 and then by 2; in case of interface function for interactions 5 and 6, B1 was modified by 0.5 and 2 and B2 by 0.5 and 2. Resulting modified H( · ) or dose-response functions for interactions 1, 2, and 4 are displayed in rows 1 and 2 and those for interactions 5 and 6 in row 3. [H( · ) functions are identified by way of relevant interfaces in GnRH-LH-Te axis in Fig. 1.] Within each display row from left to right, dose-response plots modified by 0.5, original plot, and plot modified by 2 are shown. For example, row 1, column 1 shows results of modifications of GnRH firing rate as a function of Te feedback [H1( · ) function]: original function (dotted line), maximum and minimum multiplied by 0.5 (solid line), and maximum and minimum multiplied by 2 (dashed line). (In corresponding location in Fig. 8 are 2 plots, one each for realized LH and Te concentrations; within each are a solid-line plot corresponding to 0.5-fold and a dashed-line plot corresponding to 2-fold modification of maximum and minimum values of Te feedback on GnRH pulse firing rate dose response.)

View larger version (56K): [in this window] [in a new window]   Fig. 8.   Sensitivity analysis. Simulation realizations of SDE model corresponding to each dose-response [H1( · ), ... , H5,6( · )] function modification given in Fig. 7. For each simulation, resulting concentrations of LH and Te are shown. For example, in row 1, column 1 (LH) and column 2 (Te), results of modifications of Te feedback on GnRH firing rate [H1( · ) function] are shown: maximum and minimum multiplied by 0.5 (solid line) and maximum and minimum multiplied by 2 (dashed line). Subplots correspond to those in comparable locations in Fig. 7 (in vertical pairs of LH and Te subplots) from left to right: row 1, H1( · ) and H2( · ); row 2, H2( · ) and H4( · ); row 3, H5,6( · ).

	VI. EXPLORABLE ISSUES, GIVEN A MATHEMATICALLY FORMULATED MULTINODAL FEEDBACK-CONTROL CONSTRUCT

Specific queries in GnRH-LH-Te axis pathophysiology may be addressed, we believe, using a formalized (SDE) construct of neuroendocrine feedback, such as presented here. A foremost aim is to aid in the design of experiments that are likely to help discriminate between alternative plausible explications of pathophysiology; e.g., what circadian coupling mechanisms will link 24-h variations in a deterministic input to ultradian (short-term pulsatile) neurohormone release? In this regard, our six possible models of circadian control of GnRH-LH-Te coupling suggest specific relevant experiments, such as the evaluation of diurnal changes in pituitary responsiveness to GnRH (14), in Leydig cell sensitivity to LH, and/or in GnRH's susceptibility to Te's negative feedback. Further relevant questions are as follows: How responsive are circadian variations in serum Te (and LH) concentrations to small changes in and/or different deterministic (24-h) input schemes (above)? What degree of enhanced circadian regulation is achieved when two or more (rather than a single) coupling mechanisms are implemented in concert? Which deterministic circadian linkage may be susceptible to disruption in a (particular) pathological state? Which linkages are most stable to variations in, e.g., feedback-feedforward signal strength or variability in the biological signals? Which circadian-coupling mechanism(s) emulate alterations in, e.g., puberty, aging, the menstrual cycle, and menopause?

Other physiological questions become relevant given a mathematically formulated GnRH-LH-Te feedback system. For example, how important are feedforward vs. feedback interfaces in defining the quantifiable regularity or orderliness of neurohormone release? Application of approximate entropy measures (20, 21) would be one useful tool in evaluating the following questions: How do changes in the coupling strengths (dose-response functions), intensities, or complexity of feedback and feedforward connections influence the orderliness of neurohormone release? What artifacts in pulsatile secretion or entropy estimates could be introduced by imposing strong circadian variations on the ultradian rhythms?

Another explorable issue is the seemingly more complex and biologically plastic multivalent feedback and feedforward control system of the female reproductive axis. Important questions in the female axis are as follows: How can longer-term cyclical (e.g., monthly) variations in gonadotropin output originate (e.g., by threshold, switch, or trending mechanisms) from short-term diurnal deterministic and ultradian (stochastic) elements? Whereby is overall day-to-day output stability ensured in such a formulation? How is the monthly female preovulatory surgelike increase in gonadotropin secretion generated or disrupted while still allowing preserved ultradian rhythms? What are the earliest predictable changes in system performance with female reproductive aging? How is the feedback altered plausibly in pathophysiology, e.g., the LH hypersecretory states of polycystic ovaries and testicular feminization?

A well-defined and relevantly constructed neuroendocrine feedback model of the GnRH-LH-Te axis should also eventually allow prediction of the activity and regulation of unobserved functional nodes. For example, computer-assisted simulations might help predict when measurements of LH and Te release would allow accurate prediction of (unobserved) GnRH output or GnRH-LH dose-response properties. This is significant, since GnRH cannot be measured in hypothalamic-pituitary blood in the human. What minimal experimental data would be most useful in making valid predictions of unobservable dose-response functions (e.g., GnRH feedforward on LH)? What blood sampling scheme(s) will ensure good statistical power for defining altered GnRH feedback or feedforward properties? Will computer-based simulations help guide the design of in vivo animal experimentation or clinical study? For example, one might ask, What predictions arise from the hypothesis that the GnRH pulse generator is resistant to Te's negative feedback in patients with inborn androgen-receptor mutations? Are these predictions consistent with clinical pathophysiology? What experiments could be done next to further elucidate underlying mechanisms inferred from model simulations?

Generalization of these SDE feedback concepts to other neuroendocrine axes, e.g., the growth hormone-releasing hormone-somatostatin-growth hormone-insulin-like growth factor I axis and the arginine vasopressin-corticotropin-releasing hormone-ACTH-cortisol axis, is plausible. This is practicable by modifying the core SDEs in accordance with, e.g., physiologically relevant internodal connections, particular dose-response functions, pertinent elimination rate constants, and relevant time delays. Indeed, a general SDE feedback construct may also apply on a microscopic scale to autocrine and paracrine feedback regulatory systems at the cellular level within an organ or gland and may possibly be relevant to modeling intracellular biochemical regulatory pathways. Generalization will also likely help clarify other basic issues in physiological control systems: What factors govern normal physiological variability? What is the impact of nonmonotonic dose-response/control functions and/or variably weighted or stochastically varying control functions on short-term and long-term physiological outputs? How sensitive is long-term system behavior to initial conditions, the magnitude and nature of the stochastic contributions, the choice and dispersion of time delays, the introduction of queues, thresholds, switches, or long-term trends, and the complexity of negative- and positive-feedback connections? Is additional stochastic variability in the feedback equations (see section VIII) required to model pathophysiology and, if so, under what conditions? Refinement of explicit SDEs and other formulations of physiological control systems should help address these and other issues.

	VII. DISCUSSION AND SUMMARY

Top Abstract Introduction Methods Discussion References

The salient and novel features of the current biomathematical multinodal SDE feedback control formulation of the male hypothalamic-pituitary-testicular axis include 1) physiologically dictated regulatory nodes (GnRH-LH-Te) that are linked relevantly, 2) internodal interfaces defined by appropriate feedforward/agonistic (GnRH-LH and LH-Te) and feedback/antagonistic (Te-GnRH/-LH) dose-response functions, 3) time delays in feedback and feedforward interactions among input concentrations and GnRH, LH, and Te secretion output, with pertinent individual hormone elimination rates, 4) stochastic elements embodied at potentially three levels, i.e., in the pulse generator waiting times, the system of differential equations itself (see section VIII: biological variability), and the sample observations (technical noise), thus allowing for uncertainties in the model/construct as well as the measurements, 5) continuous functions to allow arbitrary discretization to any coarse or fine time structure, 6) subsidiary linkages to incorporate plausible mechanisms of coupling between a deterministic periodic oscillator (e.g., circadian input) and the pulsatile GnRH-LH-Te feedback system, and 7) comparison of biomathematical simulations with several specific earlier published clinical experiments to assess computer-assisted predictions vs. published LH-Te responses to defined physiological perturbations. These features should offer a more plausible, realistic, and physiologically pertinent formulation of combined feedback and feedforward dynamic regulation of the male hypothalamic-pituitary-testicular axis and, by extension, other physiological neuroendocrine systems.

Our SDE biomathematical construct of the GnRH-LH-Te axis showed stability over multiple realizations, with prolonged realizations, and at high-intensity discretization (e.g., every 30 s). Mechanistic analyses revealed that, in principle, periodic input onto, or in control of, each node and/or dose-response function in this feedback control system could participate in generating or maintaining a circadian pattern of LH and Te release. More realistic were models with a periodic oscillator acting to modulate over 24 h the LH-Te feedforward dose-response curve or Te's negative feedback effects on LH secretory burst mass or basal Te secretion rates. For these linkages, circadian variations emerged in serum Te (>LH > GnRH) concentrations. Diurnal variations in the efficacy or potency of Te's feedback inhibition of GnRH pulse frequency or GnRH pulse mass provide additional possible circadian-pulsatile linkage models. Perhaps most plausible a priori would be circadian variations in GnRH pulse frequency and/or mass and/or in GnRH's sensitivity to feedback inhibition by Te (long-loop negative feedback), since a circadian oscillator mechanism resides in the suprachiasmatic nucleus with connectivity to the (mediobasal) hypothalamus where GnRH neurons originate (30). However, these mechanistic paradigms did not seem to predict strongly 24-h variations in serum Te > LH concentrations. The clinical observation that fixed pulses of GnRH at 90-min intervals result in relatively little day-night variation in serum LH concentrations in prepubertal boys, but significant (30%) variations in serum Te levels in the same individuals (13) would be consistent with a notion of 24-h variations in basal Te release or in LH-Te feedforward coupling. Both of these mechanisms yielded appropriate circadian variations in Te > LH concentrations in our simulations. Thus model simulations suggest corresponding animal experiments to help in ultimately establishing which mechanism(s) is most appropriate. Another circadian mechanism explored here require(s) diurnal differences in LH secretory responses to GnRH stimulation (not generally established in the human to our knowledge; see above). The trivial consideration that neurohormone (GnRH, LH, or Te) elimination rates or distribution spaces vary substantially over 24 h also is not documented to our knowledge. Thus the present work suggests selected relevant experiments to help clarify basic physiological principles of (deterministic) circadian-(stochastic) pulsatile system coupling within the male reproductive axis.

Simulations with the SDE construct embodying time-delayed dose-responsive positive and negative feedback within the male GnRH-LH-Te axis showed appropriate reactivity to selectively altered Te milieus and variable Te feedback potency. The computer-assisted simulations yielded qualitative inferences similar to those published independently in human and animal experiments. For example, a so-called castration response to Te withdrawal is recognized in vivo whether after administration of an androgen-receptor antagonist (17, 29) or a Te biosynthesis inhibitor (36). Our simulations also predicted augmented LH secretory burst frequency and mass in this context. Moreover, clamping GnRH pulse intervals at 120, 90, 60, or 45 min yielded higher but plateau levels of LH and Te concentrations, akin to in vivo responses in endogenous GnRH-deficient patients treated at fixed but escalating GnRH injection frequencies (25). Consequently, the coupled SDE simulator reacts appropriately to selected perturbations in the GnRH-LH-Te feedback and/or feedforward signals.

Application of the present simulator model could forseeably identify new and relevant hypotheses of altered GnRH-LH-Te network regulation underlying the depressed circadian rhythmicity in serum Te concentrations in healthy aging men (4) and the greater serial disorderliness of the LH (and Te) release process also recognized in older men (21). In addition, whereas LH pulse times have traditionally been modeled as essentially a renewal process [typically inferred from studies in normal men (5, 30)], which precludes evident feedback structure in the pulsing mechanism, the existence of strongly negative autocorrelation in successive LH interpulse interval values in midluteal phase women (24) suggests feedback, which could be explored in a suitable mathematical formulation of multivalent time-delayed feedback activity within the female GnRH-LH-progesterone/estradiol axis.

In summary, a multivalent, physiologically structured SDE model embodying relevant time-delayed and dose-dependent agonistic and antagonistic feedback and feedforward connections among principal regulatory nodes exhibits pulsatile neurohormone output that closely emulates that of the normal male GnRH-LH-Te axis. This biomathematical construct with its relevant stochastic components shows appropriate dynamic behavior after defined manipulations of selected feedback (e.g., Te) or feedforward (e.g., GnRH) signals. The SDE model also predicts several possible specific mechanisms by which circadian input may be coupled to a pulsatile axis. Extension of this formulation to the female reproductive axis and to nonreproductive feedback control systems should be possible. Moreover, its embellishment, via complementary features (e.g., in queuing theory) and as additional knowledge emerges regarding other facets of the male axis, will likely help stimulate new insights into the complex nonlinear dynamics underlying highly interactive and integrated neuroendocrine axes.

	VIII. APPENDIX

A. Elimination Rates and Feedback Interactions

1. Elimination rates. The elimination rate  of a biological molecule from a particular sampling space is related to its half-life (t_1/2) as follows: exp (t_1/2) = 1/2. Here, we assume the following 1) GnRH has a nominal t_1/2 of 1  3 min (8, 25, 30); thus _G = 0.23-0.69 min¹; 2) LH has a nominal t_1/2 of 50  80 min (30, 32); thus _L = 0.0087-0.014 min¹; 3) Te has an approximate t_1/2 of 15 min (19, 30); thus _Te = 0.046 min¹.

2. Feedback interactions. For the intact male GnRH-LH-Te axis, we allow for the following possible interactions and nominal time-delayed intervals over which the feedback occurs [these interactions are denoted by corresponding H( · ) functions: H₁( · ), ... , H₇( · )]: 1) the blood Te concentration (ng/dl) exerts a negative time-delayed (25-60 min) feedback on GnRH pulse firing rate (no. of pulses/h) (22, 26, 28); 2) the blood Te concentration (ng/dl) exerts a negative time-delayed (25-60 min) feedback on the rate of GnRH pulse-mass accumulation (pg · ml¹ · h¹); 3) the basal Te secretion rate varies with a periodicity of 24 h; 4) circulating LH concentration (IU/l) exerts a positive time-delayed (20-30 min) feedforward action on Te secretion rate (ng · dl¹ · h¹) (10, 19, 34); 5) the blood Te concentration (ng/dl) exerts a negative time-delayed (25-60 min) feedback on rate of pituitary LH mass accumulation (IU · l¹ · h¹) (11, 13, 27); 6) hypothalamic-pituitary portal blood GnRH concentration (pg/ml) exerts a positive time-delayed (0.5-1.5 min) feedforward effect on the rate of LH mass accumulation (IU · l¹ · h¹) (1, 9, 35); and 7) the hypothalamic GnRH pulse generator because of autonegative feedback exhibits a refractory period of ~1-3 min after each pulse time (30).

1

1

(12)

B. Pulse Generator and Pulse Shape

1. Pulse generator. Experimentally, the pulse times of LH closely mimic those of GnRH but are slightly time delayed (1, 8, 13, 14, 25, 27, 35). Hence, we have defined the male GnRH-LH-Te axis as having one principal pulse generator, i.e., with output (hypothalamic) in the form of GnRH, with pulse times T₀, T₁, T₂, ... . Although in young men there are often ~15-30 GnRH pulses in a 24-h period (30), in the present new model, unlike previous constructs (6, 16), the time structure of the pulses (and their amplitudes) is dynamic, being determined by time-delayed feedback from Te and GnRH itself. The latter feedback introduces a brief refractory state.

(13)

k  1

k  2

(14)

2. Pulse shape. The functions _G( · ) and _L( · ), resulting from an application of Eq. 3 to GnRH and LH, respectively, are the rates of secretion given as mass of hormone released per unit of time and distribution volume. We have used a generalized  family of densities (i.e., normalized to integrate to 1) to define the pulses and accommodate a spectrum of varying asymmetrical shapes

(15)

where ₁ > 1, ₂ > 0, and ₃ > 0 are parameters that model the secretory burst shape. The  and the Weibull families are included in this construction. Appropriate choices of  values allow on the one extreme a nearly Gaussian (symmetrical) secretory burst and, on the other extreme, a host of variably rightward-skewed representations of secretory bursts. Such asymmetry of secretion events is sometimes evident in in vitro and in vivo experiments (1, 8, 18). The maximum (rate of secretion) of ( · ) occurs at ^(r) = ₂(₁  1/₃)^1/₃. The swiftness of the upstroke is controlled by ₁ and ₃, which thus provides some measure of the amount of immediately releasable granule-contained GnRH or LH (30). The inclusion of the parameter ₃ allows for variably peaked events that are not easily accommodated by ₁ and ₂ alone. The rate of decline in secretory rate after the maximum (e.g., when hormone-containing granules are progressively depleted) is controlled by ₂ and ₃. Figure 3 (2nd row, 3rd column) displays the _G( · ) and _L( · ) used in the simulations of section V.

C. Incorporating a Circadian Rhythm

(16)

Within the foregoing structure as developed for the GnRH-LH-Te axis, six theoretically possible mechanisms for the circadian rhythms in Te and LH are 1) circadian modulation of the negative-feedback actions of Te on the GnRH pulse generator firing rate, whereby

(17)

is replaced by

(18)

2) periodic modulation of the negative-feedback effects of Te on the GnRH secretion rate (mass), such that

(19)

3) diurnal modulation of a basal secretion rate of Te itself,

(20)

4) 24-h modulation of the feedforward action of LH on the rate of secretion of Te, where

(21)

5) nyctohemeral modulation of the negative feedback of Te on the rate of secretion of LH, whereby

(22)

is replaced by

(23)

and/or 6) 24-h rhythmic modulation of the stimulatory actions of GnRH on the rate (mass) of secretion of LH, in which

(24)

is replaced by

(25)

Practically, the impact of (t) on loci 1-6 above is accomplished by shifting the sensitivity of the corresponding dose-response curve, i.e., by smoothly varying B in the logistic function (Eq. 4). Not included is a seventh formulation of GnRH autonegative feedback with circadian variation.

D. Time Discretization of the Male GnRH-LH-Te Axis for the Purpose of Simulation

Having constructed the processes up to time t_k with T_j  1 being the last pulse time, we construct T_j by solving the (discretization of the integral) equation, with T_j being the first t_k satisfying because the resulting conditional density of T_j being p[ · |T_j  1, ( · )], where Having constructed T_j, one can then calculate the jth pulse masses In the above, the circadian mechanism is mechanism 4 (24-h modulation of LH feedforward on the Te secretion rate).

Finally, what we observe is the above with "experimental" errors due to sample collecting, processing, and assaying

(26)

where the _i values are IID normal mean zero, and the variances are such that the coefficient of variation for Y^G_k, Y^L_k, and Y^Te_k is 6%.

E. Stochastic Elements

The above additional stochastic elements in the feedback system per se can be incorporated into the core equations (Eqs. 6-10) by replacing Eq. 10 with

(27)

The SDEs above correspond to the core construct (Eq. 1), now having incorporated the feedback-feedforward relationships and stochastic (biological) variability. Thus the true in vivo concentration processes, X_G(t), X_L(t), and X_Te(t) (t  0), are the resulting solutions of the above system of equations. The exact magnitude of this biological feedback system variation is not known but might vary in health and disease.

F. Mathematical Basis for SDE Feedback Time-Delay Construct of the GnRH-LH-Te Axis