|Home | About | Journals | Submit | Contact Us | Français|
Endocrine feedback control networks are typically complex and contain multiple hormones, pools, and compartments. The hormones themselves commonly interact via multiple pathways and targets within the networks, and a complete description of such relationships may involve hundreds of parameters. In addition, it is often difficult, if not impossible, to collect experimental data pertaining to every component within the network. Therefore, the complete simultaneous analysis of such networks is challenging. Nevertheless, an understanding of these networks is critical for furthering our knowledge of hormonal regulation in both physiologic and pathophysiologic conditions.
We propose a novel approach for the analysis of dose-response relationships of subsets of hormonal feedback networks. The algorithm and signal-response quantification (SRQuant) software is based on convolution integrals, and tests whether several discretely measured input signals can be individually delayed, spread in time, transformed, combined, and discretely convolved with an elimination function to predict the time course of the concentration of an output hormone. Signal-response quantification is applied to examples from the endocrine literature to demonstrate its applicability to the analysis of the different endocrine networks.
In one example, SRQuant determines the dose-response relationship by which one hormone regulates another, highlighting its advantages over other traditional methods. In a second example, for the first time (to the best of our knowledge), we show that the secretion of glucagon may be jointly controlled by the β and the δ cells.
We have developed a novel convolution integral-based approach, algorithm, and software (SRQuant) for the analysis of dose-response relationships within subsets of complex endocrine feedback control networks.
The aperiodic pulsatile nature of endocrine hormone secretion is created by an extensive feedback network of interconnected hormones. Figure 1 depicts a network that creates growth hormone (GH) pulsatility.1 The boxes in this figure represent the various compartments (i.e., nodes) of the network. For example, GH exists either in a storage compartment in the anterior pituitary or in the systemic circulation. The lines represent both the stimulatory and inhibitory control pathways (i.e., conduits) within this network.
Ideally, experimental data would be collected to document the concentration time series of each of these nodes and then simultaneously analyzed to obtain information about the properties of each of the information conduits. The complete mathematical description of this network is complex and involves many time-delayed nonlinear differential and integral equations with hundreds of parameters. This is further complicated by the experimental data being heteroscedastic (i.e., there is an unequal variability between measurements). Additionally, the concentration time series for some of the nodes are effectively unobtainable in humans, as catheters cannot be introduced into the brains of normal subjects. The complete analysis of any network of this complexity from incomplete data is thus complex or virtually impossible.
Despite these difficulties, a markedly improved under-standing of endocrine networks described in quantitative terms remains an important goal. To this end, we have developed several approaches for the analysis of small subsets within hormonal networks. The signal-response-quantification (SRQuant) algorithm and associated software presented here models an output response hormone concentration time series by several discretely sampled input signal hormone concentration time series that are delayed, spread, transformed, and combined and by a discrete convolution with an elimination function to produce an output response signal. A key feature of this software is that it was designed to be utilized by clinical investigators who do not have an exhaustive background in either numerical analysis or differential and integral equations. This algorithm is an integral part of the software package Pulse_XP, which can be downloaded for free from http://mljohnson.pharm.virginia.edu?.
This SRQuant algorithm calculates an approximation to the output hormone concentration data set based on the combination of several input signal hormone concentration time series data sets. A model is required for describing the secretion of the output hormone in terms of the input signal data set and for the translation of the combined output data set secretion into a hormone concentration time series. The algorithm adjusts the various values of the model by a weighted nonlinear least squares fitting procedure.
The input and output data sets are each discretely sampled, but they might not be sampled at the same time interval or over the exact range of data as with the data of Moenter and colleagues.2 For those data, the input signal time series is gonadotropin-releasing hormone (GnRH) sampled continuously from the portal hypophyseal blood of a ewe and pooled every 30 s. The corresponding output signal time series is simultaneously measured luteinizing hormone (LH), which was again sampled continuously from the jugular vein and pooled every 10 min.
The mathematical form of the model is a combination of nonlinear integral equations (i.e., convolution integrals) that must be integrated. Since the input function is discretely sampled, the step size for this integration need not be substantially smaller than the interval between the data points. In the present implementation, the numerical integration step size is set to one-fifth the smaller of the input data sets and the output data sampling interval. In addition, it is expected that the delay will be a continuous variable, i.e., not simply restricted to an integer number of data points. Consequently, the input and output data sets are linearly interpolated between the data points.
The next step is to apply optionally a variable width spreading function to the interpolated input series. This is accomplished by a discrete convolution of the input time series with a Gaussian function, a square-wave step function, or a triangular function that is normalized to have a unity area. The convolution with a unity area function will spread the input times series in time but will not alter the amount of mass in the input time series. Each of the input series may be spread independently with different width functions. The widths of these spreading functions are parameters that can be estimated by the fitting procedures described here.
In this context, convolution is a specific mathematical term. Specifically, Equation (1) presents the convolution of a input signal, Input(t), and a Gaussian spreading function:
Since the signals are discretely sampled, a discrete convolution is used where the integration in Equation (1) is replaced by a summation over all the data points. The discrete convolution integral at time t is numerically equal to the cross covariance of the input and the Gaussian with a lag of t.
Each input times series is optionally and independently time delayed. A common problem with explicit time delays is how to evaluate the output at data points that require a value of the input signal function before the first input signal data point because of the time delay. Some software packages will attempt to extrapolate the input signal data set to an earlier time, while others assume that the input signal time series values before the first data point are simply equal to the first data point, equal to zero, or perform an end around approximation. None of these choices is ideal or even statistically valid. In this algorithm, if a delay is expected, then the output experimental observations are truncated by the maximal expected delay so that the model will not be calculated at times earlier than the first input signal data point of each input signal time series. These delays are also parameters that can be estimated by the fitting procedures described here.
These time delays are not simply the delays expected between a hormone secretion and the appearance of change in the hormone levels. These delays are the combination of several potential factors such as a delay time within the circulation while the hormone moves to the sampling site. Another potential delay that is combined within these apparent delays is the time between a hormone reaching its target tissue and that tissue's response, which might be due, for example, to diffusion or protein synthesis.
These interpolated, spread, time-delayed input data sets are then independently transformed by signal-response functions. There are four choices within the SRQuant software for this transfer function. These are defined by Equation (2) a straight line, Equation (3) a quadratic equation, Equation (4) a stimulatory Hill-type function, and Equation (5) an inhibitory Hill-type function. Here the I(t) is the input to the signal-response function, and S(t) is the output. For each input concentration time series, the I(t) corresponds to the concentration at a particular time, t, while the S(t) is the contribution of that particular input time series to the estimated secretion into the serum of the output time series:
If more than one input time series is utilized, the S(t) outputs based on each input series must be combined to yield the estimated secretion into the blood that will be used to predict the estimated concentration of the output time series. This combination is either by multiplication or addition. The parameters of these signal-response functions are also parameters that can be estimated by the fitting procedures described here. The output of this step is the predicted secretion corresponding to the output time series.
The results of the previous step are combined with either a one-component Equation (6) or two component Equation (7) elimination function by a discrete convolution to yield a predicted concentration time series:
The numerical values of the various model parameters can be estimated by fitting predicted concentration time series to the observed output data set. This fitting operation can be performed by either a Nelder–Mead algorithm3 or a Levenberg–Marquardt algorithm.4 The Levenberg–Marquardt algorithm will only perform least squares fits, while the Nelder–Mead algorithm can perform least squares fits, least sum of absolute values fits, or min–max parameter estimations.
Simply selecting all the parameters to be fit simultaneously will likely fail. This fit is to a set of nonlinear integral equations and requires good initial estimates for the parameter values. Even with good initial estimates, the fitting process may require several steps, such as first fitting to some parameters and then fitting to other parameters to refine the initial parameter estimates. Even then, it may be impossible to estimate some of the parameters, because there might not be sufficient information contained within the data set to uniquely determine specific parameters. Because of these possibilities, the Levenberg–Marquardt algorithm was modified to perform two sensitivity analyses of the parameters at every iteration step. This algorithm requires that the derivatives of the fitting function with respect to each of the parameters being estimated be evaluated at every fitted data point. First, the parameter sensitivity is tested by summing squares of these derivatives for each of the parameters over all the data points and testing if these sums vary by more than a factor of 1012 between parameters being estimated. The results will probably not be sensitive to the parameter corresponding to the smaller sum of squared derivatives if the range of squared derivative values exceeds 1012. The cross correlation of the derivatives with respect to every pair of parameters being estimated is also tested. If any of these cross-correlation coefficients has an absolute value greater than approximately 0.9995, then the derivatives of those two particular parameters are almost perfectly correlated, and as a consequence, the parameters cannot be simultaneously estimated. Note, however, that either test may indicate that a problem exists for some particular parameter values while the problem may not exist for other specific parameter values. Since the parameter values are changing with every iteration of the Levenberg–Marquardt algorithm, these sensitivity tests are applied for every iteration. The need for these sensitivity tests is a consequence of the nonlinear nature of the fitting equations.
Some of the parameters are constrained to have physically meaningful values during the fitting operations. These constraints are implemented by a transformation of variables so that the nonmeaningful values do not exist within the number system. Specifically, the half-lives (Equations (6) and (7)), b (Equations (4) and (5)), slopes (Equations (4) and (5)), and thresholds (Equations (4) and (5)) are constrained by fitting to their logarithms instead of fitting to the parameters directly. Thus while the logarithms of these parameters can have either positive or negative values, the actual parameters are restricted to be greater than zero. The fraction of the second half-life, f2 Equation (7), is constrained to be greater than zero and less than one. This constraint is accomplished by creating a new variable, Q, in Equation (8) that can have either positive or negative values. However, after the transformation shown in Equation (8) is performed, the resulting f2 must be greater than zero and less than one:
The full-width at half-height (FWHH) of the input time series spreading functions are constrained to be greater than zero and less than a user-specified maximum value by first using a transformation of variables as given in Equation (9):
The delays for the input time series are constrained in a manner analogous to Equation (9).
To illustrate the potential applicability of the SRQuant software to the analysis of endocrine feedback control networks in general and the pancreatic hormone network in particular, we consider two examples. The first example is of two coupled reproductive hormone concentration time series. The second example studies the possibility that the secretion of glucagon by the α cells is controlled simultaneously by the secretion from the neighboring β and δ cells, insulin, and somatostatin, i.e., somatotropin release-inhibiting factor (SRIF).
An excellent example of two coupled hormone concentration time series is shown in Figure 2.2 This first example is the best illustration we could locate to explain the functioning of the analytical methods and software. While this endocrine example is not directly related to diabetes, it is a system that is well understood, the data were sampled at different sampling rates, and it emphasizes the broad applicability of the approach. The lower curve presents hypothalamic-hypophyseal portal blood of a sheep that was continuously sampled and pooled at 30 s intervals and assayed for GnRH. Simultaneously, jugular blood was continuously sampled and pooled at 10 min intervals and assayed for LH as is shown in the upper curve.
It is intuitively obvious from a cursory examination of these data that these two time series are closely related. It is generally accepted that the pulsatile concentration of GnRH stimulates the secretion rate of LH into the blood.2 However, assigning a high probability to this relationship from these data is somewhat problematic. For example, the data collection is on two different time scales, and the measurement errors are heteroscedastic and large (e.g., the GnRH assay had an optimal coefficient of variation [CV] of approximately 15%). The time intervals between secretion events are generally not constant, i.e., the pulsatile events are aperiodic and of differing amplitudes. In addition, the GnRH data predict the secretion rate of the LH into the blood, while the LH data are for the concentration in the blood, not the secretion rate into the blood.
Cross correlation is commonly utilized in the endocrine literature5 to evaluate temporal relationships between multiple hormone concentration time series. However, in the present case, it is not informative. Cross-correlation analysis requires the two time series to have equal time spacing. Thus performing a cross correlation on the data of Moenter and colleagues2 requires the GnRH data to be binned together so that the sampling interval is the 10 min interval of the LH data. After binning, the maximal cross correlation is 0.63 ± 0.19 at a lag of 20 min. This appears to be a significant relationship. However, these data also show a large autocorrelation of similar magnitude and significance. For GnRH, the maximal autocorrelation is 0.52 ± 0.15, and for LH, the maximal autocorrelation is 0.42 ± 0.15, both with a lag of 70 min. Cross correlations in the presence of autocorrelations are difficult to interpret. Alternatively, with Kendall's τ, a nonparametric measure of cross correlation,6 the 20 min lag cross correlation is 0.23 ± 0.14, which is much less significant.
Another commonly utilized indicator of relationships between coupled hormone time series is normalized cross approximate entropy (XApEn).7 The XApEn is normalized such that a value of one indicates no relationship, while a value of zero indicates a perfect relationship. The normalized XApEn from the binned data of Moenter and colleagues2 is 0.91 ± 0.14 and thus indicates that no relationship exists between the GnRH and the LH time series.
The failure of the XApEn and cross correlation to show a strong relationship between the GnRH and LH time series is a consequence of the waveforms of the events. The LH time series has a substantially longer elimination half-life than the corresponding GnRH data. As a consequence of this difference in elimination half-life, the GnRH concentrations will drop to basal values, while the LH concentrations remain elevated.
Clearly, the relationship of two irregular pulse patterns of hormones that are physiologically paired (e.g., GnRH and LH) cannot be accurately analyzed by conventional methods, and so we have developed a new analytical approach to overcome such limitations. Our SRQuant algorithm and software perform an alternate analysis in that it utilizes multiple experimental measured input data sets. For the present example, a single GnRH concentration time series predicts a single output time series (e.g., a LH concentration time series). It allows for differing sampling time intervals and number of data points, as is required for the data of Moenter and colleagues.2 This algorithm evaluates time-shifted and broadened signal-response functions (i.e., transfer functions) that are combined and then undergo a discrete convolution with an elimination function to predict the output time series. This algorithm estimates the various signal-response function parameters by weighted nonlinear regression procedures. The mathematical form of these signal-response models are consistent with the physiological processes expected for hormonal networks and signaling.
The application of this procedure to the GnRH–LH-coupled data2 is presented in Table 1 and Figure 3. For this calculation, it was assumed that the maximal delay was 8 min and the maximal spreading function was Gaussian, with a standard deviation of 15 min. The upper panel in Figure 3 consists of the LH data points and the calculated LH concentration time series based on the estimated parameters given in Table 1. Note that the first and last two data points are excluded to accommodate the maximal spreading and time delays. The lower panel in Figure 3 is the input signal (i.e., GnRH) time series. The middle panel is the predicted LH secretion that undergoes a discrete convolution with a single elimination half-life to generate the calculated LH concentration time series. The analysis described 96.5% of the variance of the LH concentration data.
The insert in the middle panel of Figure 3 is the signal-response relationship found by this analysis. In this case, it was found to be a quadratic equation with an upward curvature. This curvature is significant based on the Akaike information criterion (AIC).8 The AIC is a model selection criterion where a lower value indicates a better or more appropriate model description of data. For an analysis similar to the one shown in Figure 3 but with the quadratic term assumed to be zero, the AIC was found to be 42.48, while it was 37.58 in Figure 3 and Table 1.
Similarly, the spread FWHH and delay time are also required by the AIC. The corresponding analysis with these set to zero yielded 77.83 for the AIC. It is also of interest to consider one possible origin of these terms. Both the GnRH and LH data were sampled continually and pooled at 30 s and 10 m intervals, respectively. This means that the bloods were sampled over a different span of time with the average time for a particular LH sample occurring 4.75 min before the corresponding GnRH sample. The 4.75 min expected delay is very close to the fitted 5.5 min given in Table 1 and Figure 3. The same sampling paradigm will also require the introduction of a spreading function onto the GnRH, because the corresponding LHs were pooled and sampled over a wider interval. The expected width of this spreading is dependent upon the spreading function waveform, but it is clearly in the 10 min range. When the delay is fixed at 4.75 min and the FWHH of the spreading function is set to 10 min, the AIC was found to be 37.24 with a Gaussian spreading function, 37.30 with a triangular spreading function, and 38.31 with a square spreading function. Clearly, for these data, either Gaussian or triangular spreading functions appear equivalent and consistent with these experimental observations. Consequently, most, if not all, of the required spreading and time delay can be explained by the sample collection protocol and thus need not be the consequence of the hormonal network interactions.
In a second example, we exemplify the use of this algorithm and software to investigating the possibility that within the endocrine pancreas, the secretion of glucagon is regulated jointly by the insulin from β cells and SRIF from δ cells. The analysis uses published data in which portal vein insulin, SRIF, and glucagon were simultaneously measured in a female Sprague–Dawley rat.9 The secretion of glucagon is modeled by the product of two inhibitory Hill functions, and the formulation assumes that the release of glucagon is negatively regulated by insulin and SRIF, specifically,
The results of the application of SRQuant are shown in Table 2 and Figure 4. This secretion pattern is shown in the center panel of Figure 4. When convolved with a single pharmacokinetic half-life of 0.4 min, this secretion gives the solid curve in the upper panel of Figure 4. The vertical lines in the upper panel correspond to the experimental data with an assumed ±5% CV. The experimental insulin and SRIF data are shown in the lower two panels of Figure 4. The inhibitory Hill functions for insulin and SRIF are shown in the corresponding inserts in these two lower panels.
The second panel from the top of Figure 4 presents the standard deviation weighted residuals, i.e., the differences between the experimental glucagon data and the predicted glucagon data. It is expected that the residuals should be Gaussian distributed and random. There are many statistical tests that can be applied to assess if these residuals are actually Gaussian distributed and random. In the present case, it appears by inspection that systematic differences exist between the experimental and the predicted data between approximately 12 and 14 min. However, only one of the first fifteen autocorrelations is significantly nonzero (p < .025), and the probability of too few runs is 86%. Thus based on the autocorrelations and runs test, these residuals are not significantly nonrandom.
The cross correlations between the fitted parameters provide information pertaining to the difficulty of the numerical analysis problem. These cross correlations are a function of the range of values within the experimental data and are dependent upon the exact form of the equations that are being fit to the data. These cross correlations are between ±1, with zero being no correlation and both ±1 being intractable. Usually, cross correlations between approximately ±0.96 are reasonable. For the present example, the highest cross correlation of −0.734 is between the scaling factor and the elimination half-life.
The modeling results shown in Figure 4 and Table 2 describe 89.3% of the variance within this data, with 141.016 for the AIC with a multiplicative model. If the contributions of insulin and SRIF are added instead of being multiplied, then only 87.5% of the variance is described, with 149.921 for the AIC. The insulin without the SRIF described 84.71% of the variance, with 152.057 for the AIC. The SRIF without the insulin described 79.27% of the variance, with 166.828 for the AIC. The combination of these indicates that the best description of the data include the multiplied, not added, effects of both insulin and SRIF. This preference for the multiplied rather than the added terms implies specific mechanisms for the nature of the control of glucagon secretion.
It is also interesting to note that, for this second example, there were no explicit delays between insulin, glucagon, and SRIF. The original publication9 indicated that SRIF secretion was slightly delayed compared to that of insulin and was antisynchronous to that of glucagon. The difference in interpretation reflects the expected delay introduced by the integration of the secretion time series required to obtain the corresponding concentration time series.10
Hormonal secretion feedback control networks are typically very complex. They generally contain multiple hormones, pools, and compartments, and the hormones themselves commonly interact with multiple pathways and targets within the networks. A complete description of such networks will typically include hundreds of model-specific parameters and dozens of simultaneous differential equations. In addition, it is difficult, if not impossible, to collect specific experimental data pertaining to every hormone, pool, and compartment within the network. Thus the complete analysis of hormonal secretion feedback control networks is difficult. Nevertheless, an understanding of these networks is vital for furthering knowledge concerning hormonal regulation in both health and disease.
Typically, these, and most other feedback networks, are analyzed by formulating systems of differential equations that simultaneously describe the time course of the concentrations for each of the hormones, pools, and compartments within the network. Commonly, these models include time-delayed differential equations. These differential equation models and software packages can be daunting for the uninitiated. Within this context, we have developed an approach, algorithm, and software for the analysis of subsets of hormonal secretion feedback control networks that are not based on systems of differential equations. Our goal is to present a method that reconstructs subsets of these endocrine networks. After all subsets are reconstructed, the information can be further unified to appraise the behavior of the entire network. One of the major points of this article is that hormonal networks can be studied without formulating a dynamic feedback model and that SRQuant does not use differential equations or compartmental models. 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. The analysis of these subsets of networks is clearly more tractable than the simultaneous analysis of the complete networks.
It is also important to have the option of utilizing input signals other than actual experimental data. For example, an investigator might wish to analyze a data set where the control mechanism is thought to be a pulse generator that is amplitude modulated by a second signal. Thus one input could be a synthetic pulse generator with specific frequency and phase relationships. A second signal could be the sum of several sine and cosine waves to produce an arbitrary modulating signal. When multiplied, these would produce a modulated pulse generator. The numerical analysis will then provide information about both the pulse generator and the independent modulation signal.
To illustrate the potential applicability of the SRQuant software to the analysis of endocrine control networks, we provide two examples. In the first example, the relationship between the GnRH concentration and the secretion of LH is modeled as a quadratic polynomial. Clearly the biochemical pathways involved within this step are numerous and complex. However, given the limited amount of available data, it is clearly impossible to model these reactions completely with specific differential equations for every step of hormone synthesis, transport, and binding. Our approach is to combine these into a generic dose-response relationship, i.e., a transfer function. The most commonly used dose-response relationships are based on Hill functions, as in Equations (4) and (5). However, these are problematic when only limited amounts of data are available. Consequently, for the first example, a quadratic polynomial was employed as a limited range surrogate for the more complex Hill binding equations. Even though this system is not related to diabetes and blood glucose regulation, the example highlights the potential for the application of the software to the analysis of tightly coupled hormones for which appropriate data are available. It illustrates the advantage of our new method in comparison to other traditional approaches and the physiological information that can be obtained by reconstructing the existing dose-response interaction.
The second example is more complicated. It studies the possibility that the secretion of glucagon by the α cells is simultaneously controlled by the secretion from the insulin from neighboring β cells and SRIF from neighboring δ cells (see Figure 5). This possibility is part of the putative network mechanism suggested to regulate the pancreatic glucagon response to hypoglycemia.11 For the first time, we present evidence that the secretion of glucagon may be jointly controlled by both the β and the δ cells, which supports the hypothesized connectivity11 and demonstrates the applicability of SRQuant to problems related to blood glucose regulation and diabetes.
This algorithm and the supporting software are suitable for use by investigators who have limited fluency in differential equations. The algorithm does not contain any explicit differential equations and is based on discrete convolution integrals. The user simply defines the controlling input signals, the output signal, and the general types of interactions between the signals and the eliminations. Then SRQuant provides a rigorous statistical analysis of the particular subset of the network being studied.
The SRQuant software is an integral part of the Pulse_XP hormone pulsatility analysis software.
The authors thank Dr. Sue Moenter and Dr. Fred Karsch for supplying their 30 s sampled GnRH data used in Figure 2.