Abstract
Negativefeedback (inhibitory) and positivefeedforward (stimulatory) processes regulate physiological systems. Whether such processes are themselves rhythmic is not known. Here, we apply crossapproximate entropy (crossApEn), a noninvasive measurement of joint (pairwise) signal synchrony, to inferentially assess hypothesized circadian and ultradian variations in feedback coupling. The data comprised simultaneous measurements of three pituitary and one peripheral hormone (LH, FSH, prolactin, and testosterone) in 12 healthy men each sampled every 10 min for 4 days (5,760 min). Ergodicity, due to the time series stationarity of the measurements over the 4 days, allows for effective estimation of parameters based upon the 12 subjects. CrossApEn changes were quantified via movingwindow estimates applied to 4day time series pairs. The resultant ordered windowed crossApEn series (in time) were subjected to power spectrum analysis. Rhythmicity was assessed against the null hypothesis of randomness using 1,000 simulated periodograms derived by shuffling the interpulseinterval hormoneconcentration segments and redoing crossApEn windows and spectral analysis. By forward crossApEn analysis, paired LHtestosterone, LHprolactin, and LHFSH synchrony maintained dominant rhythms with periodicities of 18–22.5, 18, and 22.5 h, respectively (each P < 0.001). By reverse (feedback) crossApEn analysis, testosteroneLH, testosteroneprolactin, and testosteroneFSH synchrony cycles were 30, 18, and 30–45 h, respectively (each P ≤ 0.001). Significant 8 or 24h rhythms were also detected in most linkages, and maximal bihormonal synchrony occurred consistently at ∼0400–0500. Collectively, these analyses demonstrate significant ultradian (<24 h), circadian (∼24 h), and infradian (>24 h) oscillations in pituitarytestis synchrony, wherein maximal biglandular coordination is strongly constrained to the early morning hours.
 approximate entropy
 orderliness
 feedback
 patterns
physiological systems operate by coupling diverse signals through specific pathways (processes) (8, 12, 13). Pathophysiology typically results in pathway disruption, which is reflected in altered forward (amplification) and/or reverse (inhibition) control of overall output. Examples include autonomous endocrine tumors with reduced feedback responses, resulting in sustained increases in hormone secretion (9), transpubertal adaptations in pituitarygonadal coordination (39), and agingassociated adaptations of pathwayspecific orderliness (29, 31). Experimental methods have often entailed agonist and inhibitor experiments, including studies after in vitro separation of cell types (5, 35, 40). However, physical separation of otherwise integrated components definitionally disrupts regulatory linkages.
A wellestablished, modelfree metric for examining signal disruption on the basis of irregularity or complexity is the approximate entropy (ApEn) measurement (19). ApEn quantifies and discriminates subtle erosion of pattern reproducibility in system output captured as time series. In numerous biological applications, and for diverse classes of theoretical mathematical models, ApEn identifies small changes in sequential regularity of time series that correspond to changes in underlying system state (19, 27, 30, 41). A bivariate analog of ApEn, crossApEn, quantifies joint synchrony of paired signals and is thus applicable to coupled data sets from (possibly very complicated) networks (30). CrossApEn is able to distinguish insidious disruption of pathway fidelity (e.g., pituitaryadrenal and pituitarygonadal coupling) in aging individuals (29, 42). Moreover, the application of directional (A→B vs. B→A) crossApEn permits appraisal of primarily forward visàvis reverse pathway changes in health, disease, aging, and development without disturbing the underlying system (16, 17).
Although circadian, ultradian, and other rhythms are well established for meanlevel variations in hormonal dynamics, virtually nothing is known about the temporal variation of the dynamics per se. Since the study of the dynamics of hormonal secretory patterns has elucidated much toward our understanding of both physiology and pathophysiology of endocrine systems (5, 8, 35, 36), it would seem most appropriate and important to attempt to evaluate the nature of temporal variations in signal dynamics themselves. The present data set, consisting of pituitarytestis data from 4 days of intensive sampling of four hormones at 10min intervals, provides a nearly unprecedented opportunity to test the postulate that forward/reverse pathway integration varies in both a basic circadian (∼24h rhythmic) and an ultradian (rhythm with <24h period) fashion.
METHODS
Subjects.
The present analyses were applied to archival LH, FSH, and prolactin data reported earlier (14), where only univariate ApEn was assessed. Matching testosterone time series have not been described previously. In brief, 12 healthy men provided written informed consent approved by the Institutional Review Board at the Mayo Clinic, Rochester, MN. Blood was sampled from normal men every 10 min for 4 days without interruption (total blood loss: 510 ml). Three regular meals were served daily at 0800, 1200, and 1700, and room lights were turned off at 2230. Room lights were turned on at 0700, and sampling began at 0800. Subjects had lavatory access and were allowed to ambulate in the room, but they were not allowed to exercise vigorously or smoke. LH, FSH, testosterone (T) and prolactin were assayed in each of the 576 samples in each of the 12 volunteers (except for FSH, which was assayed in 11 volunteers). There were no adverse events. Figure 1A illustrates the time series data of one man.
Overview of analysis.
The data consisted of 12 subjects, with each subject having measurements taken every 10 min for 4 days (576 observations). In some contexts 12 subjects is a small sample, whereas in others it isn't necessarily so. For instance, in time series one typically observes a single timeevolving sequence of dependent observations, the values being from one subject in a population. Based upon the property of ergodicity, one is able to estimate the parameters for the population from the single subject's time series. In that setting, although there would be but one subject, the degrees of freedom are on the order of the number of observations. This is the basis for time series and spectral analysis. Ergodicity occurs in time series under the conditions of strict stationarity for which the time dependencies die out sufficiently quick. One perspective for a time series is that it consists of a large number of oscillations with random amplitudes and phases. It is this perspective that is the basis for the estimation of the ultradian/circadian/infradian rhythms of the present study. For many of the calculations of the present study, the estimates are based upon the time series for each of the 12 subjects, e.g., the power spectrum for each subject and the crossApEn calculations. We then pool these estimates over the 12 subjects, with each estimate having a large number of degrees of freedom. By so pooling the 12 estimates, large statistical power is achieved by the population estimates based upon the 12 subjects. Hence, from the time series perspective, although we only have 12 subjects, there are a large number of degrees of freedom that result in highly accurate estimates for each subject. To identify circadian (and other) rhythms, one needs data that consist of multiple days. At the same time, one must sample sufficiently fast so as to be able to detect pulse times and to estimate secretion rates. Thus, the present data, sampled every 10 min for 4 days, were ideal, and the large number of observations (576) per subject more than compensates for the relatively small sample size of 12 subjects.
Proper application of crossApEn, as for virtually all statistical measurements, is to relatively stationary or steadystate epochs. However, here the primary point of our study is to examine the potential temporal variation in feedback and feedforward signal exchange, e.g., to assess diurnal or ultradian adaptations in network and pathway function. Accordingly, to ensure proper interpretability in statistical comparisons, we require that crossApEn be applied to relatively stationary epochs and thus to relatively short windowed blocks of data. The twostep procedure of first evaluating crossApEn relatively locally in time, and then secondly sliding the underlying local windows along the entirety of the 4day records, produces a statistically valid means to evaluate temporal changes in bihormonal synchrony, which often correspond to changes in associated feedback and feedforward interactions.
Analysis of joint synchrony (crossApEn) variations over 4 days.
To test the hypothesis that LH feeds forward onto T (LHT) in a rhythmic manner, we calculated the periodogram of a movingwindow crossApEn (MVcrossApEn; details are shown below). In the appendix, definitions of ApEn, movingwindow ApEn (MVApEn), and MVcrossApEn are provided for completeness. To this end, we first standardized (performed zscore transformation of) the series of LH and T concentrations so that LH and T patterns could be compared with respect to relative reproducibility. Next, we calculated MVcrossApEn(T LH), defined as LH→T feedforward synchrony, with parameters m = 1, r = 0.2, and w = 36, where m is the pattern vector length, r is the pattern similarity threshold [expressed as a percentage of the whole (n = 576 data) series SD], and w is the window size in MVcrossApEn. Specifically, we normalized each series so that SD = 1 and then computed MVcrossApEn using r = 0.2. Note that MVcrossApEn is based on 36 points (details are given in MVcrossApEn), but r is defined on the basis of all 576 points. MVcrossApEn (T LH) yields a series of ordered crossApEn values, which can be viewed as a time series, and its periodogram is calculated. The postulate is that MVcrossApEn exhibits rhythmic changes that can be assessed from the periodogram, thus denoting systematic variations in twohormone synchrony over the 4day interval. The same concept was used for LH→prolactin and LH→FSH feedforward pairs. This analysis was applied also to the reverse direction for T→LH, T→prolactin, and T→FSH, which are termed here as feedback pathways (periphery to pituitary).
To evaluate which frequencies in the periodogram showed a power that could not be explained by chance alone, we performed a resampling procedure whereby we shuffled (resampled) interpulseinterval segments (details are given below) in the original concentration data, creating simulated series (Fig. 1B; note that in this figure the interpulse interval at 400 is shuffled to 155 in the simulated series). Pulse onsets, and thereby the intervening concentration segments, were identified initially using a previously described pulse detection method (15). It is important that interpulseinterval segments are shuffled, and not at all possible time points; in the former, the null hypothesis will be stationary processes but without the crossApEn rhythm. If all points are shuffled, the null hypothesis is an independent and identically distributed random (i.e., “white noise”) process; this is inappropriate because the observed data is obviously not iid. To obtain accurate empirical P values, the resampling process was repeated to generate a null set of 1,000 simulated 4day time series for each LHcontaining pair and each subject. For each such simulated series, we calculated MVcrossApEn and the periodogram of MVcrossApEn. The 1,000 simulated (randomly derived) periodograms were used to construct a cumulative frequency distribution from which to compute an empirical P value for the strongest rhythmic cycle in the authentic data. Given a group of k subjects, we obtained such a P value for each subject and each hormone pair also for the mean of each group for each hormone pair. Then we used Fisher's method to combine the k P values to get an overall chisquare measurement of group significance by −2 × ∑_{i}[ln(p_{i})], I = 1, . . . ,k. The overall (group effect) P value was obtained from the chisquare distribution, with degree of freedom 2 × k.
To assess the difference (in terms of the significance of the dominant period) between two pathways, for example, LH→T vs. LH→FSH, we transformed the corresponding P values in each subject to zscores and applied the twosample KolmogorovSmirnov (KS) test (zscores of LH→T vs. zscores of LH→FSH) or a matchedpair ttest (matched by subject). One can also plot the empirical cumulative frequency distribution of the zscores to visualize the difference (figures not included).
Periods of special interest (24 and 8 h) were examined by Fourier analysis to determine the amplitude and clock time of peak bihormonal synchrony. We also assessed the difference (over all periods) between two pathways by comparing the shape of their group mean spectrum. The nonparametric method of comparing cumulative periodograms by Diggle and Fisher (4) was applied (details are given below).
CrossApEn.
ApEn was developed as a quantification of irregularity in sequences and time series data motivated by both application needs (19) and fundamental questions within mathematics (27, 30). Two input parameters, a block or run length m and a tolerance window r, must be specified to compute ApEn. Briefly, ApEn measures the logarithmic probability that blocks or templates of length m that are close (within r) for m contiguous observations remain close (within the same tolerance width r) on the next incremental comparisons. The ApEn calculation provides a single nonnegative number, which is an ensemble estimate of relative process randomness wherein larger ApEn values denote greater irregularity, and smaller values to more instances of recognizable patterns or features in data. It has been widely applied to hormonal data such as these (7, 24, 29).
Analogously, crossApEn quantifies joint pattern synchrony between two time series (20, 30); the formal mathematical definition of crossApEn is given in Ref. 29. As for ApEn, it is a twoparameter family of statistics, with m and r taking the same meaning as in the ApEn setting, CrossApEn measurements, within tolerance r, the (conditional) regularity or frequency of vpatterns similar to a given upattern of window length m. It is typically applied to standardized u and v time series. Greater asynchrony indicates fewer instances of (sub)pattern matches, quantified by larger crossApEn values. For the studies below, we calculate crossApEn values for all data sets applying widely utilized and validated parameter values m = 1 and r = 20% of the SD of the specified time series. Previous evaluations, including both theoretical modelbased analysis (19, 25, 26, 29) and numerous diverse applications (16, 17, 29, 33, 37), demonstrate that the input parameters indicated above produce good reproducibility for crossApEn for time series of the lengths considered herein. In particular, the SD of crossApEn is less than or equal to 0.06 for virtually all processes analyzed for the present data lengths. This establishes the model independence for crossApEn; i.e., it assures robust qualitative inferences across diverse model configurations. In effect, the choices for m and r are made to insure that the conditional frequencies that underlie the crossApEn specification are reasonably estimated from the input data time series. For much smaller r values or larger m values than those utilized herein, one usually achieves poor conditional frequency estimates (i.e., statistical replicability is compromised), whereas for larger r values, too much detailed system probabilistic information is lost (collapsed into coarser agglomerations).
ApEn and crossApEn evaluate both dominant and subordinant patterns in data; notably, they can detect changes in underlying episodic behavior not reflected in peak occurrences or amplitudes or even subtler, more insidious changes in instances in which no apparent features, e.g., pulses or spikes, are readily identifiable (28). Additionally, ApEn and crossApEn are robust to noise and artifacts and are scale and translation independent. Importantly, changes in ApEn and crossApEn have been shown mathematically to correspond to mechanistic inferences concerning subsystem autonomy, feedback, and coupling in very distinct model settings (22, 28, 38).
Elsewhere (9), we provide a detailed stepwise description of the ApEn calculation, accompanied by several figures, that should enhance the visual understanding of the ApEn and crossApEn calculations. A recent study of paired ACTHcortisol dynamics in Cushing's disease (33) further illustrates crossApEn quantification (specifically, Fig. 1) in Ref. 33, with greater ACTHcortisol secretory asynchrony in the diseased subject compared with the control.
In Ref. 7, we provide a theoretical basis for understanding why ApEn and crossApEn provide a more general and robust measurement of feature persistence than do linear correlation and spectral measurements, either univariate or bivariate. Descriptively, correlation and spectral measurements assess the degree of matching or recurring features (characteristic subblocks) at fixed spectral frequencies, whereas the ApEn formulation implicitly relaxes the fixed frequency mandate in evaluating recurrent feature matching. Thus both ApEn and crossApEn provide a sharper measurement of equidistribution and asynchrony and can identify subtle yet persistent pattern recurrences in both data and models that the aforementioned alternatives fail to do (18, 29).
MVcrossApEn.
A movingwindow version of crossApEn was implemented to capture changes in local synchrony over time, i.e., a local time inhomogeneity of the underlying network interactions. The calculation is dynamic in providing a series of crossApEn values, since the calculation window is shifted en bloc in 10min (1sample) increments across the 4day paired hormone data. For example, consider two time series, X = (x_{i}) and Y = (y_{i}), I = 1, . . . ,N, and window size 36 samples. MVcrossApEn values are calculated on the basis of every 36 consecutive pairs (series A and series B) of data points. The first crossApEn value depends on ([x_{1},y_{1}], . . . , [x_{36},y_{36}]), and the second crossApEn value depends on ([x_{2},y_{2}], . . . , [x_{37},y_{37}]). In general, the jth crossApEn, 1 ≤ j ≤ n − 36 + 1, is calculated from data ([x_{j},y_{j}], . . . , [x_{36 + j − 1},y_{36 + j − 1}]). Thus, in MVcrossApEn, each crossApEn_{j} is calculated using only a local region of contiguous data in the 4day paired time series data.
Window size 36 (6 h) was chosen after estimates were compared for 18 (3 h), 36 (6 h), 48 (8 h), and 72 (12 h). This reflects a balance between statistical and probabilistic concerns (19). Smaller windows lose power, and larger series lose discrimination. A window size of n = 18 was too small for the current data set because there were some unreliably small crossApEn values due to lack of matches (nonreplicability), whereas a window size of n = 72, although statistical very stable, was sufficiently coarse to most effectively detect more rapid ultradian rhythms. The patterns of the MVcrossApEn series with window sizes of 36 and 48 were qualitatively quite similar. Figure 2 presents the LHT crossApEn time series with the four different choices of window size in a subject.
Periodogram.
A power spectrum, or periodogram, was used to determine the amplitudes of frequencies of periodic variations in the crossApEn time series. The goal was to quantify periodic components in the presence of noise. A peak in the periodogram indicates an important contribution to variance of the frequency (or frequencies) near the value that corresponds to the peak.
Power vs. frequency (no. of cycles/5,410 min) is plotted in the periodogram. Power is calculated as scaled square of the fast Fourier transformation of the time series. The length (period) of a crossApEn cycle (in min) for the present 4day (5,760 min) data [created from 576 separate 10min concentrations using a moving 36point (360 min) window] is 5,410 min/no. of cycles. Power was construed as significant when the value of the true observations exceeded 95% of the values of the 1,000 random samples of the interpulseinterval sequence of the time series. To summarize the results from all men, we also calculated the mean periodogram by averaging the power over all men at each frequency and analogously 1,000 mean periodograms of the random resamples to get the P value.
Resampling interpulse concentration segments.
To obtain empirical P values of the significance of a power, we broke the original time series (hormone concentrations) into pieces according to a vector of pulse onset times (15) and resampled the pieces to get a new process comprising randomly reordered pulse times as follows: given process X and its pulse onset times I, decompose X into pieces P = (p_{i}) according to I, e.g., I = (9, 53, 120, 300), then P = [(x_{1}, . . . ,x_{9}), (x_{10}, . . . ,x_{53}), (x_{54}, . . . ,x_{120}), (x_{121}, . . . ,x_{300})]; resample I according to process X [e.g., I_{X} = (53, 300, 120, 9)], and use I_{x} to resample pieces of P to get P*, e.g., P* = [(x_{10}, . . . ,x_{53}), (x_{121}, . . . ,x_{300}), (x_{54}, . . . ,x_{120}), (x_{10}, . . . ,x_{53})].
Given processes X and Y and pulse times I_{X} of X, define P_{X} and P_{Y} by I_{X}. Treat P_{x} and P_{y} as a twodim process, and resample paired pieces of P_{X} and P_{Y} to get P*_{X} and P*_{Y} simultaneously.
Again, as stated above, it is important that interpulse interval segments are shuffled and pieced back together to produce a shuffled series and not a shuffling of all time points. Roughly speaking, the end of one pulse interval and the starting of another is very much like a renewal process. The segments of a hormonal concentration profile that are closest to being independent are the interpulse intervals.
Comparison of periodograms.
We compare the shape of two periodograms by comparing their normalized cumulative periodograms (4). The details of the analytical procedure are illustrated by the following example.
We compare the two periodograms of moving block CrossApEn series LT (LH→T) and LF (LH→FSH). Let I_{LT}(ω_{j}) and I_{LF}(ω_{j}), j = 1, . . . , [(541 − 1)/2] be the power of corresponding crossApEn series at frequency ω_{j} (j cycles/5,410 min), and 541 is the length of the crossApEn series.
Calculate the normalized cumulative periodograms, denoted as F_{LT}(ω_{j}) and F_{LF}(ω_{j}), j = 1 . . . 270 as F_{LT}(ω_{j}) = ∑_{i = 1 to j} I_{LT}(ω_{i})/∑_{i = 1 to 270} I_{LT}(ω_{i}) and F_{LF}(ω_{j}) = ∑_{i = 1 to j} I_{LF}(ω_{i})/∑_{i = 1 to 270} I_{LF}(ω_{i}).
Calculate the KS difference as d = sup_{j}∣F_{x} (ω) − F_{y} (ω)∣, j = 1, . . . ,270.
Randomly interchange I_{LT}(ω_{j}) and I_{LF}(ω_{j}) for j = 1, . . . ,270, to get I*_{LT}(ω) and I*_{LF}(ω) and then calculate F*_{LT}(ω), F*_{LF}(ω) and their KS difference d*.
Repeat step 3 1,000 times, and calculate the total number of those d* at least as large as d and divide by 1,000 to get the significance probability.
RESULTS
Twelve men participated and completed all 4 days of 10min blood sampling without interruption. One subject did not have sufficient serum to permit FSH measurements. The median (range) age and body mass index were 33 (18–75) yr and 26 (21–29) kg/m^{2}, respectively. Deconvolution analysis of the uninterrupted 4day LH time series in the 12 subjects yielded (per 4 days) 50 (42–67) LH secretory bursts, rapid and slowphase LH halflives (min) of 22 (6.5–25) and 90 (71–119), an LH secretory burst mode (min) of 15 (12–20), basal and pulsatile LH secretion (IU·l^{−1}·4 days^{−1}) of 143 (64–289) and 198 (125–315), mass per burst (IU/l) of 3.6 (2.3–7.2), and γ (a regularity term): 2.7 (2.1–3.2), for which values > unity denote greater regularity of the Weibull process than a Poisson process.
CrossApEn values without windowing were determined first. These are single values calculated over all 4 days (576 data points) for six physiologically relevant pairs. Since crossApEn values are asymptotically normal under broadly valid assumptions, Table 1 gives the mean ± SE (n = 11, 12) crossApEn (and the median) for the pairwise (LH, FSH, prolactin, and T) combinations studied here. By paired parametric testing, feedback crossApEn values were higher (denoting lower joint pattern synchrony) than feedforward crossApEn values: TLH crossApEn > LHT crossApEn (P < 0.01), TPRL > LHPRL (P < 0.01), and TFSH > LHFSH (P < 0.01). The descending rank order of joint synchrony was LHPRL (maximal synchrony and lowest crossApEn) > LHT = LHFSH > TFSH > TLH = TPRL.
To test the hypothesis that joint synchrony varies over time, windowed crossApEn (1, 0.2, 36) was evaluated repeatedly in 10min increments across the 4 days in each subject for all 6 hormone pairs. Figure 3 shows crossApEn time series of LHT and all other hormone pairs for one subject. The total number of successive crossApEn values is 541 (spanning 5,410 min or ∼90 h) in each series, reflecting the window size and original series length.
Simulated random crossApEn time series were obtained by shuffling the multivariate process by LH (for feedforward) and the T (for feedback) interpulse segments of the original concentration data 1,000 times each. This produced 1,000 (random) crossApEn time series each containing 541 estimates for each subject (n = 12 or 11) and each hormone pair (n = 6). Fourier analysis (power spectrum estimation) was then performed on the authentic crossApEn series derived from the unshuffled hormone pair and on the 1,000 randomly generated crossApEn time series. This allowed direct comparison of each authentic crossApEn power spectrum (plot of squared amplitude against frequency) with the distribution of 1,000 random crossApEn power spectra (periodograms) in the 12 subjects.
At each frequency, a mean power spectrum was calculated by averaging individual spectra and can be compared with its randomly shuffled counterparts to get a significance measurement. Figure 4 depicts the results of six crossApEn power spectrum analyses. Each of the six graphs in this figure presents the mean power spectra (mean periodograms) of one hormone pair. The solid (irregular) curve gives the mean spectrum for the authentic (unshuffled) crossApEn series and the dashed (smoother) curve the 95% limit for the random series. The xaxis defines the number of cycles (frequency) per 5,410 min (or ∼90 h) and the yaxis the strength of the periodic signal. Period length (h) of the dominant cycle detected in the crossApEn time series is listed below the rectangles. The P value within each graph designates the probability of observing the detected dominant cycle in 1,000 shuffled pairs. The periodicities with the greatest power (each P ≤ 0.001) were 18–22.5 (LHT), 22.5 (LHFSH), and 18 h (LHPRL) for feedforward and 30 (TLH), 30–45 (TFSH), and 18 h (TPRL) for feedback, suggesting mainly circadian (nearly 24 h) and ultradian (more than 1 cycle/24 h) rhythms.
The periods with the greatest power (primary period) for each subject and each pair of hormones are summarized in Table 2. Those periods are not necessarily the same across subjects for a given hormone pair synchrony. For the null hypothesis of no significant periodicities of crossApEn in the cohort as a whole, chisquare values were 89 (22 df), 104 (24 df), and 115 (24 df) for LHFSH, LHprolactin, and LHT (each P < 10^{−6}), and pooled zscores were 6.8, 7.4, and 8.1, respectively (each P < 10^{−6}). Chisquare values for T feedback crossApEn cycles were comparable at 80 (24 df) for TLH, 110 (22 df) for TFSH, and 112 (24 df) for Tprolactin (each P < 10^{−6} vs. chance periodicities). Corresponding cohortdefined zscores were comparably significant, namely 5.9, 7.9, and 7.8 (each P < 10^{−5}).
Although most of the primary periods are detectable (P < 0.05) for each subject and each synchrony, their level of significance may differ across the distinct synchronies. We compared synchronies by their levels of the significance of the primary periods (Table 3). P values were first transformed to zscores, with a higher score corresponding to a higher level of significance. By paired ttest (P = 0.008) and by the KS test (P = 0.002) of the zscores, the primary period of LHT feedforward variation is more significant than that of the TLH feedback variation in joint synchrony. In addition, primary TFSH crossApEn rhythms are generally more significant than primary TLH rhythms (P = 0.018, paired ttest; P = 0.016, KS test). Primary Tprolactin periods are more significant than primary TLH feedback synchrony by ttest (P = 0.05) or by post hoc nonparametric testing (P = 0.01).
Significant (P < 0.05) ranges of periodicities, which include 24 h (∼4 cycles/541 observations), were disclosed for all linkages except T→LH (where 3 and 5 cycles/541 observations are significant). Significant 8h periodicities (∼11–12 cycles/541 observations) were obtained for all linkages except (LH→FSH). Cosine regression is illustrated in Fig. 5, and phase results are tabulated in Table 4. These show maximum synchrony at 0400–0500 for many 24h rhythms and for most 8h rhythms.
To visually compare the shapes of the mean spectra of hormone pairs, we created Fig. 6 by plotting the normalized cumulative mean periodogram of the first hormone pair vs. the normalized cumulative mean periodogram of the second hormone pair. The cumulative periodogram was normalized so that the range was 0 to 1. If the two periodograms are identical, the black solid curve should fall along the equal line of identity (black dashed line). The gray dashed curves create an envelope of the black solid curve by plotting its analog, using the randomly interchanged data. A black solid curve falling outside of the envelope corresponds to a significant difference between the two hormone pairs. If the xaxis denotes the measurement of the first series and the yaxis denotes that of the second one, a black solid curve significantly below the equal line indicating the power spectrum of the first is shifted to lower frequency compared with that of the second one. Roughly speaking, the relative power (compared with higher frequencies) of the low frequencies is greater in the first series.
Four comparisons are presented in Fig. 6, and the P values of all nine comparisons are listed in Table 5. Four were found to be significantly different: LH→T vs. LH→FSH (P = 0.033), LH→T vs. LH→PRL (P = 0.028), LH→FSH vs. LH→PRL (P = 0.002), and T→FSH vs. TPRL (P < 0.001). The power spectrum of LH→FSH was shifted to lower frequencies compared LH→T and LH→PRL. LH→T was shifted to lower frequencies compared with LH→PRL. The shape of T→FSH power spectrum was shifted to lower frequencies compared with T→PRL. Note that this method allowed us to compare the shape of the periodograms, and thus the comparison is over a combination of all periods.
DISCUSSION
The present analyses unveil for the first time that quantifiable synchrony of reproductive hormonal linkages varies considerably across the day and night. This is true for both 24 and 8h rhythms and for almost all bihormonal pathways studied here. We repeated the analysis using partial data (2.24 days with 323 data points, starting at 2 PM on the 1st day); the 24h rhythm was still detected in all pathways, and the 8h rhythm was detected in LHPRL and TPRL pathways. The consistency of these findings is remarkable and suggests entrainment of reproductive hormonal networks by the suprachiasmatic nucleus and other Zeitgebers. The phase analysis using full data showed that maximal synchrony occurs at ∼0400–0500 in most linkages. Maximal regularity of hormonal linkages just prior to awakening may reflect a physiological process that partly regulates or is influenced by the sleepwake transition. The exact relationship between paired hormone coordination and sleep or circadian pacemakers could potentially be assessed if sleepwake cycles were desynchronized from circadian time. Although substantial knowledge of circadian clock function is available in other physiological contexts (3), whether the rhythms that we have observed are energy responsive and temperature compensated, as observed in animal models (1, 2), will require further elucidation.
A strength of the crossApEn approach is that complete system specification or a full network model is not necessary to certify the statistical validity of changes in ApEn (or crossApEn) with concomitant changes in joint synchrony and thereby the network function (22). This work also introduces an empirical estimator of random feedforward relationships, e.g., LHT crossApEn, calculated repeatedly after (LH concentration) pulse segments were resampled (performing 1,000 random permutations on the original pulse onsetdefined intervals). The empirically constructed cumulative frequency distribution was used to estimate the probability of the authentic crossApEn rhythm arising from chance LHT pattern similarities in each of the 12 subjects. By this construction, the cohortdefined zscore for the probability of falsely rejecting the null hypothesis of no crossApEn rhythms was estimated as z = 8.1 (P < 10^{−8}). Likewise, for LHFSH and LHprolactin, group zscores were z = 6.8 (P < 10^{−6}) and z = 7.4 (P < 10^{−7}). Thus, according to this estimation strategy, joint hormone feedforward (pituitarygonad and pituitarypituitary) synchrony exhibits remarkable nonrandom rhythmicity over a 4day interval. The individual factors supervising these changes will require further study, including such entities as food intake, exercise, sleep, and stress.
The present data show that bivariate hormone synchrony oscillates in an infradian (periodicity >24 h) fashion in normal individuals for TLH and TFSH feedback. The latter rhythm was foreshadowed in a rat model of Tinduced nycthemeral FSH cycles (Fig. 2 in Ref. 10). Moreover, androgens are capable of modifying the master circadian pacemaker (11). Whether infradian rhythms in feedforward synchrony exist with periodicities exceeding 4 days (5,760 min) cannot be ascertained from our data (which is only 4 days in length). Important infradian rhythms operate in a variety of physiological systems (43, 44), but precisely how they relate to the basic feedback and feedback synchrony cycles discovered here is not yet known.
Forward crossApEn oscillated at a median ultradian (>1 cycle/24 h) periodicity of 18 h in the case of LHprolactin synchrony. Whether such joint synchrony arises from prolactin's feedback effects on the hypothalamus (34) or from intrapituitary gonadotropelactotrope interactions is not yet known (37).
Caveats include the relatively small cohort size (12 subjects), which is more than sufficiently compensated by the degrees of freedom resulting from the 2,308 separate hormone measurements in each person, restriction of analyses initially to men, the limitation of sampling only every 10 min for 4 days, and the ultimate need to appraise whether feedback rhythms are free running, Zeitgeber entrainable, and temperature compensated. The methods presented here should offer an analytical platform for pursuing such studies, including in other feedback systems such as cellular and molecular networks. The requirement is that relevantly paired signals be monitored frequently without biological coupling being disrupted and that sufficient data be obtained to capture several cycle lengths.
In conclusion, reproductive hormone homeostasis, evaluated at the level of bivariate signaling pathways, proceeds by way of infradian, circadian, and ultradian dynamics. If verified further, these outcomes introduce a new basis for experimental and modelbased investigations of feedback/feedforward regulation in health, premorbidity, disease, and aging.
Perspectives and significance.
Endocrine glands communicate with each other and with remote target tissues via intermittent signal exchange (6, 12, 13, 40). Whereas parametric models of single hormone signals are widely utilized in the endocrine literature, modeling the joint dynamics of more than one hormone (e.g., the pituitary and target gland hormones together) is not well developed. This is due mainly to the fact that biological interactions are typically nonlinear, and the modeling of nonlinear dynamics imposes severe technical limitations. Linear methods such as auto and crosscorrelations will not necessarily detect nonlinear interactions (29). Hence, methods such as crossApEn, which can detect nonlinear relationships without requiring explicit modeling of the nonlinear structure, are important (21, 23, 32). This concept was used to show, for the first time, that the dynamics of hormonal secretory patterns follow circadian and ultradian patterns with maximal bihormonal synchrony at about 0400–0500 just before waking. The magnitudes of the observed changes across the 24h day are similar to those observed in many other nycthemeral rhythms. Feedforward linkages were uniformly more synchronous than feedback linkages among reproductive hormones. The consistency of these findings suggests broader relevance to other classical (hypothalamopituitaryend organ) hormonal networks and introduces the need to conduct experiments at a strictly consistent time of day, to evaluate bihormonal regulatory complexity in other physiological contexts, and to explore pathophysiological disruption and molecular bases of underlying mechanisms.
GRANTS
This work was supported in part via Center for Translational Science Activities Grant No. 1UL1RR024150 from the National Center for Research Resources (Rockville, MD) and by AG031763 and DK050456 (Metabolic Studies Core of the Minnesota Obesity Center) from the National Institutes of Health (Bethesda, MD).
DISCLOSURES
The authors have nothing to disclose.^{ . . . }
ACKNOWLEDGMENTS
We thank Jill Smith and Ashley Bryant for their support in manuscript preparation and Ashley Bryant for data analysis and graphics, the Mayo Immunochemical Laboratory for assay assistance, and the Mayo research nursing staff for implementing the protocol.
APPENDIX: DEFINITIONS OF ApEn AND MOVINGWINDOW VERSIONS
ApEn can be defined for a deterministic or stochastic realvalued sequences of finite length and for an infinite realvalued sequence, assuming the limit exists. Let m be a positive integer and r a positive real value. Let x_{t} = {u_{t}}, t = 1, . . . ,N, be a sequence of length n. Define x_{t} = [u, . . . ,u_{t + m}] for t = 1, . . . ,N − m, and the distance between x_{t} and x_{s} as the maximum difference in their respective components:
If d_{m + 1}(x_{t},x_{s}) ≤ r, the two segments of length m + 1 are said to be close. The basis of ApEn and crossApEn is the quantification of the predictability of future values based upon the past. For a given x_{t}, the fraction that is close is
Let
In C_{t}^{m + 1}(r) above, x_{t} is being compared with all possible x_{s}, s = 1, . . . , N − m.
In the present study, we are interested in movingwindow versions, where we will define the window length via a parameter w; the window at time t, 1 ≤ t ≤ N − m − w, is the set of times s such that 0 ≤ s − t ≤ w. To define MVApEn, let
Once the C_{t}^{(w,m + 1)}(r)(vu) and C_{t}^{(w,m + 1)}(r)(uv) are defined, the construction of the MVcross − ApEn(v u) and the MVcross − ApEn(uv) proceeds exactly as in that of the MVApEn.
Footnotes

↵* S. M. Pincus is an independent mathematician in Guilford, CT.
 Copyright © 2011 the American Physiological Society