Biomathematical modeling of pulsatile hormone secretion: a historical perspective

A logical extension of this improved understanding of both the fundamental physiological underpinnings of pulsatile hormone release noted earlier and the potential importance of pulsatile secretion in terms of clinical applications was the recognition of the need to define and subsequently characterize in an objective manner hormone pulses within concentration-time series. Examples of some initial efforts to achieve this goal included defining pulses as an increase followed by a decrease of an arbitrary number of units (e.g., 5 mIU in the case of LH; Yen et al ., 1972 ) or an increase in the serum concentration of a hormone of at least 20% over the proceeding nadir (the so-called “Santen and Bardin” algorithm; Santen and Bardin, 1973 ). For whatever the reason, however, nearly a decade passed before the power of computer-assisted algorithms for the identification and characterization of pulsatile hormone secretion—and importantly the impact of how blood sampling protocols may affect the ability to define such pulses—began to receive investigative attention.

An alternative method proposed by Keenan and Veldhuis (1997) to appraise hormone secretion uses a combination of a model of random, but structured, variations in pulse amplitudes superimposed on basal hormone release with a nonstationary Poisson process generating the timing of the hormone pulses. This approach utilizes stochastic differential equations to evaluate certain properties of irregular hormone time series, including pulse onset times, basal secretion, and pulse shape, and has the advantage of using an automated procedure for initial estimate of the pulse positions. However, in this method the amplitude of each pulse is assumed to depend on the previous interpulse length, which is a restrictive condition that may not be valid in all hormone release control systems.

In the same general time frame that attention was increasingly being focused on how the regulation of individual hormone-secreting components might be appraised in quantitative terms, methods with which to define the activity of coupled hormonal systems were being proposed. One of the earliest such approaches utilized cross-correlation, a standard and widely applied method for estimating the degree to which two data series are related. It was argued that cross-correlation could be used to support the existence of a relationship between two nodes in the hormone network and indeed this method has been applied successfully in various circumstances ( Booth et al ., 1996 ; Veldhuis et al ., 1988 ). However, the straightforward use of cross-correlation suffers from several disadvantages, including spurious values due to autocorrelation and the linear nature of the assumed relationship. These limitations prompted the development of other methods, specifically adapted to endocrine applications. Within this context, and shortly after the development of ApEn to appraise regulation of a single hormone-secreting node, Pincus and Singer (1996) suggested that the signal-response relationship between coupled nodes within a hormonal network could be described by estimating their cross-approximate entropy. Indeed, the program Cross-ApEn assesses the degree of asynchrony between two hormonal concentration-time series known to be coupled and determines whether the coupling is more “ordered” or more “random.” Taking the LH–testosterone system as an example, application of Cross-ApEn demonstrated more irregularity or “asynchrony” in data obtained from older vs younger men, even though the mean serum concentrations of LH and testosterone were not found to differ as a function of age ( Pincus et al ., 1996 ). Thus, Cross-ApEn reinforced the notion that information of potentially significant physiologic interest could be obtained by quantitating the interactions of coupled systems within the overall hormonal network. Quite recently, approaches using distinctly different mathematical methodologies from those subserving cross-correlation or Cross-ApEn have been developed. Farhy and colleagues have proposed a method that seeks to define and quantitate relationships whereby one hormone modulates the release of a second hormone ( Nass et al., 2008 ). This method (referred to as FeedStrength ) overcomes some of the problems associated with the application of the more traditional statistical methods (e.g., cross-correlation) for estimation of the relationships between two simultaneously sampled hormones. To test the hypothesis that hormone A modulates the secretion of hormone B, FeedStrength first determines the individual secretion bursts of hormone B using deconvolution procedures. The amplitude of each secretory burst of hormone B is compared to the average concentration of hormone A over a fixed time interval during and before each hormone B burst. Rather than comparing the entire hormone A and B concentration profiles (as would be the case if cross-correlation is used), the method relates the amplitude of an individual hormone B burst to the concentration of hormone A accompanying the development of this burst. This approach has the advantage that it disregards those parts of the release profile of hormone A during which there is no pulsatile secretory activity of hormone B. Consequently, the method is suitable for analyzing endocrine relationships in which one of the hormones is presumed to be an amplifier/modulator and not the sole regulator of the secretory activity of the second hormone. The application of FeedStrength to estimate the amplifying effect of ghrelin on growth hormone secretion resulted in the first clear demonstration of a relationship between these two hormones in normal subjects ( Nass et al. , 2008 ). That this approach is also suitable for estimating the variability in feedback relationships was proposed when the method was applied to synthetic data ( Farhy et al. , 2006 ). As shown in Fig. 14.6 , network feedback modeling was used to create a system mimicking the feedback effect of progesterone on LH secretion. Synthetic data were modeled to mimic serum LH (hormone B) stimulating progesterone (hormone A) secretion in normal midluteal phase women. The network also assumes a self-entrained GnRH pulse generator (G) stimulating the secretion of B, an ovarian-hypothalamic delayed feedback component, and a system-independent external regulation (F). Estimates of the level of feedback control (LFC) were computed using FeedStrength as described ( Nass et al. , 2008 ). Also shown in Fig. 14.6 are results of the application of this method to synthetic data sets in which the strength of the feedback signal (A-to-G) is progressively degraded and the external regulation (F-to-G) becomes dominant. As can be seen, the LFC for the A-to-B coupling identifies correctly gradual changes of the impact of A on the system relative to the effect of F across 11 different levels (L0–L10): (i) minimal LFC when only F, but not A, affects G (L0); (ii) a gradual LFC increase when the feedback becomes progressively dominant over the impact of F (L1–L9); and (iii) maximal LFC when only A, but not F, controls G (L10). Thereby, the LFC as calculated by FeedStrength provides a robust estimate of the feedback control exerted by one hormone on another when the main pulse-generating hormone is unknown and can estimate feedback alterations associated with such factors as age, gender, and pathophysiologic conditions. An additional approach with which to assess the possibility that one hormone regulates the secretion of another hormone in a quantitative manner is to use dynamic modeling. To this end, Johnson and colleagues (2008) developed a computer-based approach for Signal-Response Quantification ( SRQuant ), which reconstructs the dose–response interactions between hormones by fitting parameters of a dynamic model to experimental data. Within this framework, several input signals can be individually delayed, spread in time, transformed by dose–response relationships, combined, and then convolved with a pharmacokinetic elimination function to predict the time course of the concentration of an output hormone. As an example, SRQuant has been used to estimate the effect of the endogenous opioid peptide antagonist naloxone on GnRH (and consequently LH) secretion using an animal model (the ovariectomized ewe) in which GnRH and LH were measured simultaneously in blood samples obtained from the hypothalamic–hypophyseal portal circulation (GnRH) and the jugular vein (LH)( Goodman et al ., 1995 ). A minimal model was constructed based on the assumption that GnRH stimulates the secretion of LH, which is, in turn, eliminated after its release into the circulation. The model has several physiological parameters (rate of elimination and secretion of LH, ED 50 , and slope of the GnRH–LH dose–response interaction)( Farhy, 2004 ), which were determined simultaneously by SRQuant via an iterative procedure to achieve the best possible fit of the model output to the LH data if GnRH was used as input. The result of this process is illustrated in Fig. 14.7 and depicts a very good coincidence between the model predicted and observed LH concentration: 80% of the variance was explained, confirming the existence of a strong feed-forward dose–response relationship between GnRH concentration and LH secretion in this animal model. Further analysis showed that before naloxone administration, the model provides a much better fit of data by explaining approximately 88% of the variation in LH, whereas during the treatment only 68% of the variation was explained. In addition, a clear change in the “shape” of the dose–response interaction was detected after the treatment. This outcome suggests that at least part of the effect of naloxone treatment on the GnRH–LH axis reflects a direct effect on GnRH–LH relationships. It is interesting to note that the analysis of such GnRH–LH relationships with other methods, including cross-correlation and Cross-ApEn , does not provide satisfactory results: In this particular case the cross-correlation was only 0.64, which would account for about 36% of the variance in LH by the linear correlation model. The (ratio) Cross-ApEn values were 0.94 ± 0.08 ( m = 1) and 1.07 ± 0.2 ( m = 2), which are typical for a random rather than a strong dose-dependent relationship between the two hormones. Thus, SRQuant provides a viable alternative and a powerful addition to the more traditional methods for the appraisal of coupled endocrine systems.

Whereas the methods described earlier have allowed for a detailed appraisal of systems in which hormone-secreting component A stimulates or inhibits hormone-secreting component B or vice versa, it is obvious that in many cases endocrine systems function in a dynamic fashion in which there are multiple feedback interactions among the system components. Given that such systems fulfill commonly accepted criteria to be considered as networks , an increasing amount of investigative interest has focused on defining such networks in quantitative terms. A separate approach with which to address this issue and as adopted by Farhy and colleagues ( Farhy, 2004 ; Farhy and Veldhuis, 2005 ; Farhy et al ., 2007 , 2008 ; Kovatchev et al ., 1999 ) and others ( Bergman et al ., 1979 ; Dalla Man et al ., 2005 ; Keenan and Veldhuis, 1997 , 1998 , 2004 ; Wagner et al ., 1998 ) suggests that endocrine networks can be appraised by deterministic differential equations that describe the rate of change of individual hormones and their interaction. Utilization of this method has the potential to test whether a physiological regulatory hypothesis is consistent with the experimental observations and typically results in formulation of physiological assumptions (additional to the main hypothesis) that would explain the observed data specifics ( Farhy, 2004 ). In effect, this is in silico testing of physiological hypotheses, which addresses specific features of the in vivo behavior of endocrine systems. This approach permits estimation of the integrative action of individual inferred or postulated feedback and feed-forward interactions within an endocrine axis and combines multiple steps ( Farhy, 2004 ), including (i) preliminary analysis of available experimental data, (ii) design of a formal network describing the axis linkages and its approximation with coupled differential equations, and (iii) determination of the model parameters either functionally or by fitting the model output to available data. As an example, application of this approach may allow insights into the physiological parameters of a specific endocrine system, which may subserve its oscillatory behavior. Evaluation of simple two- or three-hormone feedback networks, such as the one shown in Fig. 14.8 , illustrates the capacity of delayed nonlinear feedbacks to generate hormone pulsatility ( Farhy, 2004 ). In these examples (that may appear in different real endocrine axes), hormones (A, B, and C) are secreted continuously and are subject to elimination. However, in contrast to the previous examples, two or more (delayed; D 1 , D 2 ) feedback control loops exist between the system components. It is these loops that generate the oscillatory secretory patterns. For example, the two delayed control loops in the network shown in the middle left panel of Fig. 14.8 can drive a particular hormone oscillatory pattern if certain physiological conditions are met. As is depicted in Fig. 14.9 , this complex multiple-volley secretory pattern is typical for the growth hormone dynamics in the adult male rat, and results of the simulations show that the underlying construct ( Fig. 14.8 , middle left panel) wherein B may be thought to be growth hormone (GH), A GH-releasing hormone, and C somatostatin, represents one minimal set of interactions capable of reproducing the typical in vivo GH pattern in the male rat. However, as evident from additional model-based analysis, this network is too simplistic and needs to be extended and modified to address other key experimental observations ( Farhy and Veldhuis, 2005 ; Farhy et al ., 2007 ). Another model-associated outcome could be a detailed reconstruction of the endocrine network dose–response interactions (similarly to the use of the aforementioned procedure SRQuant ). This would be the case when suitable data exist to permit the precise iterative determination of the model parameters ( Bergman et al ., 1979 ; Dalla Man et al ., 2005 ). In this case the physiological hypothesis that the postulated network unifies the primary endocrine axis control mechanisms would be strongly supported if the model can provide a good fit of data based on objective statistical criteria.

Given the complexity of endocrine axes, a logical extension of previous attempts to assess isolated signal–response endocrine interactions would include the development of automated methods with which to reconstruct the dose–response interactions within an entire hormone network. To this end, a dynamic modeling method, such as that described earlier, could be utilized within a computer-based model similar to SRQuant. This approach could make possible the approximation of endocrine networks exhibiting multiple feedback interactions and the subsequent reconstruction of their parameters based on existing experimental data. The result would be a powerful tool capable of testing regulatory hypotheses and estimating the alterations in hormone release control mechanisms related to various factors including, but not limited to, inter- or intraspecies differences, natural processes such as aging, pathophysiologic situations, and interventional procedures.

The second half of the 20th century witnessed an extraordinary number of scientific advances that coalesced to provide a much improved understanding of how endocrine systems in general—and the reproductive axis in particular—function as “networks” in which individual components communicate with one another in a highly coordinated fashion. Central to the reproductive system are cells within the anterior pituitary gland, which secrete two so-called “gonadotropins”—luteinizing hormone (LH) and follicle-stimulating hormone (FSH)—which control both gametogenesis and hormone production in the testes and ovary. Although our recognition of the role played by LH and FSH in relation to gonadal function dates back to the 1920s ( Fluhmann, 1929 ), it was midcentury before the possibility that the hypothalamus might serve to regulate pituitary function, including that of the gonadotropin-secreting cells, was put forward by Green and Harris (1947) . As depicted in Fig. 14.1 , these investigators suggested that neuronal systems within the hypothalamus could both synthesize and secrete chemical (hormonal) messengers, which could then travel to the pituitary gland via the hypothalamic–hypophyseal portal circulation where they would influence the secretion of pituitary hormones. However, another two decades would elapse before the first of these hypothalamic hormones, including gonadotropin-releasing hormone (GnRH), which stimulates the synthesis and release of LH and FSH, and thyrotropin-releasing hormone (TRH), which is responsible for thyrotropin stimulating hormone (TSH) production, were isolated, characterized, and synthesized by Guilleman ( Burgus et al ., 1972 ) and Schally ( Schally et al ., 1971 ; Matsuo et al ., 1971 ). The significance of this work was recognized by the Nobel Prize for Physiology or Medicine being awarded to these investigators in 1977. Shortly after the discovery of GnRH, and capitalizing on the availability of a relatively new measurement technique known as radioimmunoassay, studies in the rat ( Gay and Midgley, 1969 ), Rhesus monkey ( Atkinson et al ., 1970 ), and the human ( Midgley and Jaffee, 1971 ; Santen and Bardin, 1973 ; Yen et al ., 1972 ) suggested that LH and FSH were secreted in a pulsatile rather than a tonic fashion. The profound physiologic significance of this observation was demonstrated by Knobil and colleagues using the Rhesus monkey as a model ( Knobil, 1980 ). Knobil’s seminal studies demonstrated that whereas the administration of GnRH in a pulsatile fashion to monkeys in which the GnRH-secreting cells had been ablated did indeed effect the secretion of LH and FSH, the administration of the peptide as a nonvarying infusion failed to stimulate secretion of the gonadotropins ( Fig. 14.2 ). Moreover, GnRH infused into monkeys with an intact hypothalamic–pituitary system actually diminished the secretion of the gonadotropins. These studies strongly suggested that GnRH itself might well be secreted in a pulsatile manner, thus stimulating the pulsatile secretion of LH and FSH. Approximately a decade later, studies in which GnRH and LH were measured simultaneously in the rat ( Levine and Ramirez, 1982 ) and sheep ( Clarke and Cummins, 1982 ; Levine et al ., 1982 ) within the hypothalamic–hypophyseal portal circulation and in the vasculature draining the anterior pituitary gland confirmed this direct relationship. Subsequent basic science and clinical studies have served to underscore the extraordinary significance of the pulsatile nature of hormone secretion within the reproductive axis. Using the rat as a model, Marshall and colleagues have shown that the gene expression of LH and FSH is controlled by the frequency and amplitude characteristics of the pulsatile GnRH signal ( Haisenleder et al ., 1994 ). Within the clinical setting, the significance of the characteristics of the GnRH signal has been exploited by clinicians who wish either to restore normal gonadotrope function by means of the exogenous pulsatile administration of GnRH (e.g., in cases of primary hypothalamic failure such as Kallmann’s syndrome or hypothalamic failure-associated infertility) or to inhibit gonadotropin secretion by means of a long-acting GnRH agonist (e.g., in the setting of prostrate cancer, precocious puberty, or endometriosis), which provides a nonpulsatile signal to the pituatary and serves to desensitize the pituitary and essentially abolish LH and FSH secretion.

While deconvolution procedures represented a major advance in our ability to model hormonal data, it became readily apparent that the earlier versions had significant limitations, including the subjective nature of the choice of candidate secretory bursts, the lack of robust statistical verification of resolved secretory bursts, and the user-unfriendly interface of the programs. To address these concerns, a novel deconvolution procedure known as AutoDecon has been developed and validated for LH ( Johnson et al ., 2008 ) and growth hormone. This program finds the optimal number of secretory bursts and initial parameter estimates while simultaneously performing the deconvolution. Moreover, the fact that AutoDecon is fully automated and statistically based addresses the issues noted earlier, including the user-unfriendly nature of prior deconvolution methods and the previous lack of statistical verification of resolved secretory events. Analysis with Auto-Decon results in a substantially higher true-positive rate of identification of hormone secretory events than analysis with Cluster . Moreover, both false-positive and false-negative rates are lower with AutoDecon compared to Cluster and sensitivity is substantially higher with AutoDecon .

Following the development of methods for hormone pulse and secretory burst identification as described previously, attention focused on whether quantitative techniques could be utilized to appraise the regulation of individual “nodes” (i.e., hormone-secreting components) within overall hormone networks. Such components often respond to a primary input signal but are, at the same time, subject to modulation by a separate signal in the form of “feedback” from the target gland. Again using the reproductive network as an example, LH- and FSH-secreting gonadotropes within the anterior pituitary gland respond to GnRH from the hypothalamus (acting as the primary signal), but the pituitary response to the GnRH signal may be influenced significantly by ambient concentrations of hormones such as estrogen and progesterone secreted by the ovary (i.e., the target gland for LH and FSH). An early attempt to characterize regulation of a single hormone secretory component sought to determine, using quantitative methods, the degree of “orderliness” within a hormone concentration-time series, which was considered to reflect the degree to which the gland is regulated by both primary and feedback signals. To assess such orderliness, Pincus (1991) developed a mathematical approach with which to estimate “approximate entropy” within the data series. This approach (known as ApEn ) is complementary to hormone pulse and secretory burst activity in that it evaluates both dominant and subordinate patterns in concentration-time series. Given its ability to detect changes in underlying episodic behavior, which will not be reflected in pulse occurrences or their characteristics, ApEn acts as an explicit “barometer” of feedback system change in many coupled systems. In addition to its application to components of the reproductive system ( Pincus et al ., 1997 ), ApEn has been utilized to probe the regulation of other hormones, including growth hormone ( Friend et al ., 1996 ; Hartman et al ., 1994 ), aldosterone ( Siragy et al ., 1995 ), cortisol ( van den Berg et al ., 1997 ), and insulin ( Meneilly et al ., 1997 ).

The early to mid-1980s saw a flurry of activity focused on the development of a range of sampling protocols (e.g., protocols in which blood samples were obtained anywhere from every 1 to every 60 min over intervals of 2–24 h) and of computerized pulse analysis techniques, each subserved by quite different mathematical constructs. The latter included Ultra developed by Van Cauter and colleagues (1981) , Pulsar by Merriam and Wachter (1982) , and Cycle Detector by Clifton and Steiner (1983) . The use of sampling protocols of differing intensity and interval—and analysis of the resultant hormonal concentration-time series with computer-assisted methods based on nonuniform mathematical and statistical approaches—perhaps not surprisingly resulted in quite diverse estimates of pulsatile hormone secretion, in many cases under theoretically similar physiological conditions. Thus, application of these earlier pulse detection programs to LH concentration-time series obtained in normal women at several phases of the menstrual cycle resulted in the identification of anywhere from 14 to 24 LH pulses/24 h in early follicular phase women, from 17 to 29 in late follicular phase women, and from 4 to 16 in midluteal phase women ( Backstrom et al ., 1982 ; Burger et al ., 1985 ; Crowley et al ., 1985 ; Filicori et al ., 1984 , 1986 ; Murdoch et al. , 1989 ; Reame et al ., 1984 ; Rossmanith et al ., 1990 ; Schweiger et al ., 1989 ; Soules et al ., 1984 ). To define the impact of sampling protocol and of pulse identification method on the ability to detect LH pulses in men and women, we applied several of the available pulse detection algorithms to serum LH concentration-time series obtained from normal men ( Urban et al ., 1988 ) and normal women at different phases of the menstrual cycle ( Evans et al ., 1992 ). Both the sampling protocol and the pulse detection algorithm were found to have a profound effect on the number of pulses detected. With regard to sampling intensity, up to an eightfold difference in the number of pulses detected within the same data set could be found when the complete data set (sampled at 4-min intervals) was progressively “deintensified” by omitting successive samples to simulate sampling at, for example, 12-, 20-, and 40-min intervals. Similar concerns were apparent when data from normal women during the midluteal phase of the menstrual cycle were analyzed with four of the independently developed computer-based techniques mentioned previously. As is shown in Fig. 14.3 , when data obtained with a standard protocol (samples obtained at 10-min intervals over a 24-h interval) were analyzed, estimates for the number of pulses varied from 3 to 17 per 24 h. Such divergent results prompted the development of yet more pulse detection algorithms. Nearly simultaneously, two novel programs became available to the scientific community: Detect , written by Oerter, Guardarraso, and Rodbard ( Oerter et al ., 1986 ), and Cluster , provided by Veldhuis and Johnson ( Veldhuis and Johnson, 1986 ). Although using distinctly different mathematical approaches, both programs shared the properties of being rigorously statistically based. Whereas Cluster identifies statistically significant increases and decreases in the concentration-time series with a sliding T test, Detect defines pulses as statistically significant increases and decreases in the numerical first derivatives of the concentration-time series. Of interest, when these programs were applied to the LH concentration-time series mentioned earlier, the number of LH pulses identified was statistically indistinguishable.