## Abstract

We develop, implement, and test a feedback and feedforward biomathematical construct of the male hypothalamic [gonadotropin-releasing hormone (GnRH)]-pituitary [luteinizing hormone (LH)]-gonadal [testosterone (Te)] axis. This stochastic differential equation formulation consists of a nonstationary stochastic point process responsible for generating episodic release of GnRH, which is modulated negatively by short-loop (GnRH) and long-loop (Te) feedback. Pulsatile GnRH release in turn drives bursts of LH secretion via an agonistic dose-response curve that is partially damped by Te negative feedback. Circulating LH stimulates (feedforward) Te synthesis and release by a second dose response. Te acts via negative dose-responsive feedback on GnRH and LH output, thus fulfilling conditions of a closed-loop control system. Four computer simulations document expected feedback performance, as published earlier for the human male GnRH-LH-Te axis. Six other simulations test distinct within-model coupling mechanisms to link a circadian modulatory input to a pulsatile control node so as to explicate the known 24-h variations in Te and, to a lesser extent, LH. We conclude that relevant dynamic function, internodal dose-dependent regulatory connections, and within-system time-delayed coupling together provide a biomathematical basis for a nonlinear feedback-feedforward control model with combined pulsatile and circadian features that closely emulate the measurable output activities of the male hypothalamic-pituitary-Leydig cell axis.

- neuroendocrine
- biomathematics
- stochastic differential equations
- reproductive axis

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*+ d*t*) (from Ref. 16) - d
*W*_{i}(*t*) - Stochastic noise term (via differential of Brownian motion) (see section VIII)
- d
*X*_{i}(*t*) - Incremental change in concentration of
*hormone i*at*time t* *Z*_{i}(*t*)d*t*- 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
- (
*l*_{j,1},*l*_{j,2}) - Time-delayed interval for
*j*th feedback interaction - λ( ⋅ )
- Cosine function specifying periodic (circadian) input
- Λ( ⋅ )
- (Stochastic) pulse generator intensity function
- LH
- Luteinizing hormone
- Pulse mass
*j*for*hormone i* *p*(*t*)- Pulse generator (stochastic) intensity
- ψ
_{i}( ⋅ ) - Pulse shape for
*hormone i* - SDE
- Stochastic differential equation
- Pulse
*time j*for*hormone i* - Te
- Testosterone
*W*_{i}( ⋅ )*i*th Brownian motion process (one of “noise” processes) (see section VIII)*X*_{i}(*t*)- Concentration of
*hormone i*at*time t* - Observed concentration of
*hormone i*at*time t*_{k} *Z*_{i}(*t*)- Secretory rate for
*hormone i*at*time t*

## I. GENERAL METHODS

### A. Background Physiology

An important initial question in capturing physiological behavior of the male hypothalamic-pituitary-gonadal axis in a biomathematical construct is: What is an appropriate level at which to formulate the model? We have chosen to focus on the rate of change in blood hormone concentration, since this is a principal measurement variable in health and disease. Here, we do so in continuous time. By evaluating resultant hormone concentrations one can test predictions of the biomathematical formulation experimentally via available measurements, and by structuring in continuous time one can compare data under different sampling schemes.

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.

#### 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
Equation 1where *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).

Next we incorporate feedback and feedforward relationships within the hypothalamic-pituitary-Leydig cell axis by defining mathematically, via “*H*” functions, how each hormone secretion rate (at any instant in *time t*) of *Z*
_{G}(*t*),*Z*
_{L}(*t*), and*Z*
_{Te}(*t*) depends, in a nonlinear manner, on (all or some of) the prior ( ⋅ ) hormone concentrations [*X*
_{G}( ⋅ ),*X*
_{L}( ⋅ ), and*X*
_{Te}( ⋅ )] over appropriately time-delayed intervals. In addition, we will allow for the stochastic variability, which occurs in several respects within the axis.*Equation 1
* will first be replaced by stochastic differential equations (SDEs), *Eqs. 6-10,* which will allow for a stochastic (feedback-driven) pulse generator. In section VIII*E*, a more extended stochastic formulation is presented, which, in addition to the foregoing, attempts to describe further in vivo biological variability. Experimental uncertainty also arises from sampling and measurement errors due to sample collection and processing and subsequent assay of the hormone concentration. For the LH axis, such sampling-related technical variability has been estimated as 3–15% (29-36). Biological variability would influence the time evolution of the actual unobserved “true concentration” levels, whereas technical variability does not. Both factors affect any individual “realized” continuous concentration series or its discrete-time sampling.

Our feedback construct of the in vivo male GnRH-LH-Te axis thus consists of a system of coupled SDEs, one for each of GnRH, LH, and Te. The SDEs are of the most basic type, i.e., random ordinary differential equations, in that the forcing function is now a random process; in section VIII a more extended formulation is considered. The elimination function acting on hormone concentrations will be assumed initially to be exponential (31, 32). We will impose time-delayed, nonlinear (dose-responsive) feedback by relevant hormone concentrations on future secretion rates and on the intensity of the point process, which describes the GnRH pulse generator’s bursting activity (below). We posit that these principal biological features produce the observed oscillatory nature of such dynamic physiological systems.

## 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

### A. Constructing GnRH-LH-Te Secretion

The pituitary gland releases a small, approximately time-invariant amount of LH into the blood, which is often referred to as constitutive or basal secretion (1, 8, 25, 30), here designated as β. In addition, regulated (nonbasal) LH secretion is driven by feedforward (GnRH) and inhibited by feedback (Te) interactions within the GnRH-LH-Te axis. We assume that this potentially regulated secretion of LH can occur in two forms: pulsatile release and continuous release. To accommodate pulsatile LH release, we assume a strict dependence of each LH pulse on a GnRH pulse signal (8, 13, 25, 35). We designate activity of the so-called GnRH pulse generator via a stochastic point process, which has a probability of activating a GnRH (LH) pulse in the next time increment (*t* + δ*t*) that varies on the basis of time-delayed negative feedback by Te. Because the feedback is time delayed, our point process will not be a Markov process (see section VIII), and the distribution of waiting times will not be simply exponential.

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*
_{2}to *t* − *l*
_{1}, as transduced via some dose-response function *H*
_{4}( ⋅ )
Equation 2where β_{Te} is the Te basal release rate. Thus we model Te release via the function*H*
_{4}( ⋅ ) above; the precise assumed form of *H*
_{4}( ⋅ ) is given in sections IV*A* 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*
_{0}, *T*
_{1},*T*
_{2}, … , where without loss of generality we will assume *T*
_{0} = 0.

To accommodate GnRH and LH synthesis and accumulation in storage granules, let *M*
^{j} denote the amount (mass) of hormone accumulated from the last pulse time (*T*
_{j − 1}) to the present pulse time (*T*
_{j}) this accumulation will be available for release at *time T*
_{j}. Then, to define the time course of hormone release, let ψ( ⋅ ), where ψ(*s*) = 0, *s* ≤ 0, called the pulse shape, be a function normalized to integrate to 1, which represents the instantaneous rate of secretion in units of mass of hormone released per unit time and distribution volume.

Let us represent a pulse at *time t*, having started at pulse*time T*
_{j}, by a function*M*
^{j}ψ(*t* − *T*
_{j}), where *M*
^{j} is the mass (or amplitude) of the pulse. The overall rate of secretion at *time t*,*Zt*), is then assumed to be given by the sum of basal (β) and (summated) pulsatile secretion
Equation 3where *H*( ⋅ ) designates feedback control functions that regulate the mass of GnRH or LH accumulating. [In section IV*A* these will be*H*
_{2}( ⋅ ) in the case of GnRH and*H*
_{5,6}( ⋅ ) in the case of LH.] This is a plausible approximation of LH release in many physiological contexts (1, 8, 9, 11-14, 19, 22, 24-27, 30-36). In sections IV and VIII we explicitly expand on the form of the above*H*( ⋅ ) control functions and the feedback interactions.

### B. Circadian Modulation

In addition to basal and pulsatile GnRH and LH release and continuous Te secretion, there is a prominent circadian (24-h) rhythm in serum Te concentrations. In men, serum Te concentrations exhibit an ∼24-h cycle, with maximal values at around 4–6 AM, dropping later in the day and evening by 15–70% and then rising again in the night (34). An important issue is how to appropriately incorporate this “rhythm” in the GnRH-LH-Te axis in a manner consistent with experimental data.

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*
_{0} and *B*
_{1} and phase shift φ_{1} being chosen so λ( ⋅ ) would vary between a high of 1 at 4 AM and a low of 0.6 at 4 PM (*B*
_{0} = 0.85, *B*
_{1} = 0.15, φ_{1} = −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, η_{0} denotes a minimum amount of mass always secreted per pulse, and η_{1} 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*)/d*t* 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, η_{1}, 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]

The detailed mathematical motivations of the present formulation, i.e., elimination rates, feedback interactions, pulse shape, and pulse generator, and time discretization of the output are given in section VIII. Here we present the (time-delayed) feedback-feedforward dose-response functions that embody the interactions within and confer nonlinear dynamics on the (male) pulsatile GnRH-LH-Te axis. In section IV*B* we present six possible adaptations of this basic structure, each incorporating a circadian rhythm in a different mechanistic manner.

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*
_{1} ≤ *t*
_{2} < ∞, 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
Equation 4If *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*
_{1}( ⋅ ) for the GnRH pulse firing rate as a function of Te concentration,*H*
_{2}( ⋅ ) for the rate of GnRH pulse-mass accumulation as a function of Te concentration,*H*
_{4}( ⋅ ) 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*
_{7}( ⋅ ), 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*
_{3}( ⋅ ), 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*
_{1}( ⋅ ),*H*
_{2}( ⋅ ),*H*
_{4}( ⋅ ), 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 *j*th 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^{−1} ⋅ h^{−1}) would be of the form
Equation 5Also, 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 Figure2 illustrates the corresponding mathematical functions.

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 VIII*A*along with their motivation. In the sensitivity analysis (section V*C*), 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 VIII*B1*).

### B. Basic SDE Construct

Our formulation of the male GnRH-LH-Te feedback system is thus governed by the following coupled SDEs. The pulse times of GnRH (*T*
_{1}, *T*
_{2}, …) are the result of a point process with a so-called stochastic pulse intensity
Equation 6where *H*
_{7}(*t*) is a refractory condition whereby GnRH release may be inhibited by autofeedback (see section VIII). The conditional density for*T*
_{k} given*T*
_{k − 1} and Λ( ⋅ ) will be required to satisfy
Equation 7The secretory pulse masses for GnRH and LH, respectively, are given by
Equation 8
Equation 9The basic components of the male axis (GnRH, LH, and Te) are thus
Equation 10The SDEs in *Eqs. 6-10
* correspond to the core construct (*Eq. 1
*), now having incorporated the feedback-feedforward relationships and a stochastic pulse generator (feedback-driven) 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. In section VIII*E*, a generalization of *Eq. 10
* is discussed with additional stochastic biological variability in the feedback equations per se.

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

We can include a circadian rhythm acting deterministically via several possible mechanisms, in which there is deterministic circadian input to the GnRH-LH-Te feedback axis at any of one or more nodes or control functions: *1*) Te (negative) feedback on GnRH burst frequency,*2*) Te (negative) feedback on GnRH burst mass, *3*) basal Te secretion rate, *4*) LH (positive) feedforward on Leydig cell Te secretion, *5*) Te (negative) feedback on LH secretion mass,*6*) GnRH (positive) feedforward on LH secretion, and *7*) GnRH autonegative feedback.

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 V*B, 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*).

### A. Implementation of Possible Circadian Rhythm Mechanisms

As illustrated in Figs. 3, 4
*A*, and 5
*A*, individual and mean (of 500) realizations of pulsatile LH, Te, and GnRH output show a range of 24-h periodic rhythms that emulate circadian variations in these neurohormone concentrations. Figure 3, *top* displays individual and mean (of 500) realizations of the (instantaneous) GnRH pulse firing rate (no./day) for each of the six possible circadian rhythm models and for simulations of clinical *experiments 1–3* described below (see sections V*B* and VIII*C*); the GnRH pulse generator was “clamped” in feedback *simulation 4,* and hence there is no entry for it in Fig. 3, *top*.

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. 4
*A*, *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. 4
*B*). 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. 4
*C*, *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. 4
*D*). Finally, periodic variation over 24 h of GnRH’s feedforward (stimulatory) action on LH release evokes relatively little diurnal rhythmicity in Te concentrations (Fig. 5
*A*). 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

To explore further the performance of this multinodal SDE feedback model of the GnRH-LH-Te axis, the following four computer-assisted simulations, corresponding to known prior clinical experiments (as referenced below), were performed. The incorporation of the circadian rhythm for these experiments was via *mechanism 4,* i.e., the 24-h periodic modulation of LH feedforward on Te secretion rate (discussed in section IV*C*).

#### 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. 5
*B* 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^{−1} ⋅ h^{−1}(vs. 1,350 ng ⋅ dl^{−1} ⋅ h^{−1} in nominally normal conditions). The corresponding simulations correctly predicted (Fig. 5
*C* 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 V*B2*) 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. 5
*D*, 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. 6
*A*). 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

To ascertain the dependency of the structural dynamics of our SDE feedback model on the parameters of the dose-response functions, we varied their respective values in two respects (cf. Fig. 2): *1*) the maximum and minimum of the dose-response functions were simultaneously increased, then simultaneously decreased, and *2*) the midpoint of the dose-response function was shifted to the left, then to the right.

In particular, for feedback-feedforward *interactions 1*,*2*, and *4*, each of their corresponding dose-response functions is parameterized by values *A*, *B*
_{1},*C*, and *D*; in the case of *interactions 5* and*6*, the dose-response function was parameterized by *A*,*B*
_{1}, *B*
_{2}, *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*
_{1} first by 0.5, then by 2; in the case of *interactions 5* and *6*, we modified*B*
_{1} by 0.5 and 2 and *B*
_{2} 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.

## 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

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}) = ½. 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^{−1}; *2*) LH has a nominal *t*
_{1/2} of 50 − 80 min (30,32); thus α_{L} = 0.0087–0.014 min^{−1};*3*) Te has an approximate *t*
_{1/2} of 15 min (19, 30); thus α_{Te} = 0.046 min^{−1}.

#### 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*
_{1}( ⋅ ), … , *H*
_{7}( ⋅ )]:*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^{−1} ⋅ h^{−1});*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^{−1} ⋅ h^{−1}) (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^{−1} ⋅ h^{−1}) (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^{−1} ⋅ h^{−1}) (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).

The formulation of *interaction 7,* a refractory condition, denoted by *H*
_{7}( ⋅ ), will not be via a dose-response function; this is discussed in section VIII*B1*; also, *interaction 3,* allowing Te basal secretion rate to vary with a 24-h periodicity, will not be part of the basic GnRH-LH-Te axis but represents one possible mechanism for incorporating a circadian rhythm (see *Eqs. 16
* and *
20
*). Accordingly, for the simplified male axis, there are four principal feedback dose-response functions: *H*
_{1}( ⋅ ) for Te feedback on the GnRH pulse firing rate (*interaction 1*),*H*
_{2}( ⋅ ) for Te feedback on the rate of GnRH pulse-mass accumulation (*interaction 2*),*H*
_{4}( ⋅ ) for the rate of LH-driven Te secretion (*interaction 4*), and*H*
_{5,6}( ⋅ , ⋅ ) for Te feedback and GnRH feedforward on the rate of LH pulse-mass accumulation (*interactions 5* and *6*). We then use empirical dose-response logistic functions to define how such feedback interacts with the various axis components
The values of *A*, *B*
_{1}, *C*, and *D*and *A*, *B*
_{1}, *B*
_{2},*C*, and *D* used in the computer experiments were empirically determined on the basis of presumed normative physiology in the healthy man (2, 4, 5, 13, 14, 17, 19, 21, 25, 29-36); the values were 2.06, −0.005, 28, and 3, for*H*
_{1}( ⋅ ), 3.57, −0.008, 60, and 1 for *H*
_{2}( ⋅ ), −2.76, 0.8, 900, and 10 for *H*
_{4}( ⋅ ), and 2, −0.0074, 0.35, 5, and 1 for*H*
_{5,6}( ⋅ , ⋅ ). Sensitivity analysis (Figs. 7 and 8) was used to explore the physiologically relevant limits and impact of varying absolute parameter values on key measures of GnRH, LH, and Te output. The time delays for the above *j*th feedback interaction will be expressed throughout by *l*
_{j,1} and*l*
_{j,2}. So far, in vivo lag times are not well established, but they have been estimated by cross-correlation analysis or by catheterization (1, 2, 7-11, 14, 19, 21, 25, 30,34).

We integrate various processes over time-delayed intervals. Consequently, until *time t* is above the maximum time delay, the feedback will not originate over the full time-delay interval but, rather, only from the amount of time since *time 0*; to accommodate this we could use the notation *s* ∨ 0 to define the maximum of *s* and zero. For example, the feedback of Te concentration (in ng/dl) on the rate of GnRH pulse-mass accumulation (in pg ⋅ ml^{−1} ⋅ h^{−1}) would be of the form
Equation 12However, for simplicity we have not used this notation, but such limits (*s* ∨ 0) have been analytically and computationally implemented.

### 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*
_{0}, *T*
_{1},*T*
_{2}, … . 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.

To formulate the refractory condition of *interaction 7*mathematically, let *N*(*t*) denote the number of pulses up to and including *time t,* for then*T*
_{N(t)} is the time of the last pulse. Let ε be a small positive value. The time delays for*interaction 7* are *l*
_{7,1} = 0 and*l*
_{7,2} = 1 min. We can define the refractory condition as *H*
_{7}(*t*)
Thus, at a *time t* during the interval (*T*
_{j}, *T*
_{j} + *l*
_{7,2}) the value of *H*
_{7}(*t*) will be the small value ε; at other times the value of*H*
_{7}(*t*) will be 1. (In the simulations of section V, ε was taken to be 0.)

In the present model the evolving concentrations of Te and GnRH exert feedback on the GnRH pulse generator through a pulse intensity function. For a stationary or nonstationary Poisson process, the intensity function λ is a fixed deterministic function; it is not using any feedback. In addition to including feedback, one would like to be able to control (in a probabilistic sense) the pattern of interpulse lengths. The Weibull distribution has an additional parameter for that purpose, call it γ, and as γ increases it forces a more regular pattern (the Poisson is the special case γ = 1). In our formulation the function λ is now a stochastic process (because of the feedback), not a deterministic function. The present formulation of the stochastic pulse generator allows for a modulation of the instantaneous rate of pulsing by time-delayed feedback, in addition to control over the regularity of interpulse length. Just as the deterministic intensity function in the nonstationary Poisson can be viewed as a deterministic time transformation of the stationary Poisson, our model can be viewed as a stochastic time transformation of a Weibull renewal process. In the simulations (section V) the value of γ was taken to be 2, chosen to produce more regularity than the Poisson but still allowing a fair amount of flexibility.

More precisely, there is a “pulse generator” intensity Λ(*t*), for which Λ(*t*)d*t* describes the probability of firing in the infinitesimal time increment (*t*, *t* + d*t*)
Equation 13In Fig. 3 (1st row, 1st column), we show the feedback function*H*
_{1}( ⋅ ) of Te on frequency output of the GnRH pulse generator. The conditional density for*T*
_{k} given*T*
_{k − 1},*T*
_{k − 2}, … ,*T*
_{0} and Λ( ⋅ ) will then be required to satisfy
Equation 14We allow the probabilistic structure of the pulse times to depend on time-delayed feedback and flexibility in the variability of the interpulse lengths when the mean is (roughly) known. As a consequence of the time delays, the *T*
_{k} values cannot satisfy a Markov property (at least not in a finite-dimensional sense). However, conditional on the pulse intensity Λ( ⋅ ), the *T*
_{k} values will, in fact, be Markov; hence, our pulse generator construct incorporates time-delayed feedback while at the same time allowing for relative simplicity. Also, the parameter γ ≥ 1 allows for flexibility in the variability of the interpulse lengths when the mean is (roughly) known. It was the inflexibility of the stationary and nonstationary Poisson processes that suggested the above generalization. If Λ( ⋅ ) were a deterministic function and, hence, feedback is absent, then the pulse generator reduces to a Markov process, including the stationary and nonstationary Poisson processes as special cases.

#### 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
Equation 15where β_{1} > 1, β_{2} > 0, and β_{3}> 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)} = β_{2}(β_{1} − 1/β_{3})^{1/β3}. The swiftness of the upstroke is controlled by β_{1} and β_{3}, which thus provides some measure of the amount of immediately releasable granule-contained GnRH or LH (30). The inclusion of the parameter β_{3} allows for variably peaked events that are not easily accommodated by β_{1} and β_{2} alone. The rate of decline in secretory rate after the maximum (e.g., when hormone-containing granules are progressively depleted) is controlled by β_{2} and β_{3}. Figure 3 (2nd row, 3rd column) displays the ψ_{G}( ⋅ ) and ψ_{L}( ⋅ ) used in the simulations of section V.

### C. Incorporating a Circadian Rhythm

Given the above basic construct, we can consider the inclusion of a circadian rhythm acting deterministically via several possible mechanisms. The effects of the 24-h rhythm are observed primarily in the cyclical pattern of serum Te concentrations over the day and, to a lesser extent, in LH. An unresolved physiological issue is how the circadian rhythm might interface with various pulsatile and feedback loci within the GnRH-LH-Te axis. A general form of any (24-h) rhythmic function, composed of *m* harmonics, is simply
Equation 16where appropriate values of the amplitudes*B*
_{0}, … , *B*
_{m} and the (possible) phases φ_{1}, … , φ_{m} will ensure the positivity of λ ( ⋅ ); in the simulations of section V, one harmonic (*m* = 1) was used, with *B*
_{0} = 0.85,*B*
_{1} = 0.15, and φ_{1} chosen so that the maximum of λ( ⋅ ) occurred at an appropriate time (e.g., 4 AM).

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
Equation 17is replaced by
Equation 18
*2*) periodic modulation of the negative-feedback effects of Te on the GnRH secretion rate (mass), such that
Equation 19
*3*) diurnal modulation of a basal secretion rate of Te itself,
Equation 20
*4*) 24-h modulation of the feedforward action of LH on the rate of secretion of Te, where
Equation 21
*5*) nyctohemeral modulation of the negative feedback of Te on the rate of secretion of LH, whereby
Equation 22is replaced by
Equation 23and/or*6*) 24-h rhythmic modulation of the stimulatory actions of GnRH on the rate (mass) of secretion of LH, in which
Equation 24is replaced by
Equation 25Practically, 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

The discretization of the above system of SDEs, by use of an Euler scheme, results in the SDEs that follow. Let Δ*t* =*t*
_{i + 1} −*t*
_{i} be the scale of the discretization (e.g., Δ*t* = 30 s); this is ordinarily smaller than the actual sampling increment (e.g., sampling every 10 min). Here, we take the sampling increment to be the same as the scale of discretization, Δ*t* = 30 s. A sequence of IID uniform (0, 1) random variables is created (*U*
_{k}, *k* ≥ 1), which will be used in the construction of the pulse times.

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*j*th 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
Equation 26
where the ε_{i} values are IID normal mean zero, and the variances are such that the coefficient of variation for
,
, and
is 6%.

### E. Stochastic Elements

The hypothalamic-pituitary-Leydig cell axis is a highly coupled system with variability within its dynamics. Such variability arises not only via internodal interactions but also at the levels of *1*) the stochastic GnRH pulse generator, which is further perturbed by feedback, and *2*) and multicellular hormone (GnRH, LH, or Te) secretion and subsequent mixing within the blood. In addition, random experimental variations arise in the course of sample withdrawal, processing, and analytic assay. The motivation of our extension from the deterministic differential equations (*Eq. 1
*) to the SDEs (*Eqs. 6-10
*) is to allow for precisely such aggregate within-system variability in a continuous time formulation. Using the probabilistic concept of Brownian motion, we can formulate the component of biological variability resulting from instantaneous secretion arising from nonuniformly secreting cells that are variously arrayed about capillaries and the subsequent mixing of secreted neurohormone molecules within the turbulent bloodstream (30); the use of Brownian motion to describe such variation in diffusive behavior of fluids (e.g., mixing and turbulence) can be mathematically justified because of the cumulative effects of a large number of interactions (e.g., cellular and molecular), each occurring on a very small time scale. The biological variables, such as molecular admixture in the bloodstream, intra- and intercellular nonuniformities in hormone synthesis and secretion rate, local intraglandular variations in microvascular perfusion nonuniform cellular exposure to paracrine and autocrine regulatory molecules (distinct from GnRH, LH, and Te), and unequal secretory cell energy stores (30), can be recognized in aggregate by inclusion of an additional stochastic (Brownian) contribution in the SDEs (see below).

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
Equation 27The 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

The mathematical basis for the above formulation (*Eqs. 6
*,*
9
*, and *
27
*), as well as that of *Eqs. 6-10,*is currently being studied. There is a proper probabilistic framework in which such a system, with all the imposed time-delayed, dose-responsive feedback interactions, stochastic elements, etc., is achievable, with the realizations from the resulting processes,*X*
_{G}(*t*),*X*
_{L}(*t*), and*X*
_{Te}(*t*) (*t* ≥ 0), being continuous real functions. Such mathematical verification shows that the various structures (e.g., feedback and pulsing), which are “local” in nature, produce the proper long-term stable behavior. The present model, in which time-delayed feedback from the concentration processes, and not just the history of the prior pulse times themselves, modulates the pulse generator, is novel to stochastic modeling in general. For general references on point processes and SDEs, see Refs. 3 and 23.

## Acknowledgments

D. M. Keenan was supported by Office of Naval Research Contract 00014-90-J-1007. J. D. Veldhuis was supported in part by the National Science Foundation Center for Biological Timing, the National Institutes of Health General Clinical Research Center, NIH Research Career Development Award 1K04-HD-00634 and P-30 Reproduction Research Center Grant HD-28934 from the National Institute for Child Health and Human Development.

## Footnotes

Address for reprint requests: J. D. Veldhuis, Div. of Endocrinology, Dept. of Internal Medicine, Box 202, University of Virginia Health Sciences Center, Charlottesville, VA 22908.

- Copyright © 1998 the American Physiological Society