System for quantifying asynchrony between signals5846189Abstract A data processing environment processes a plurality of signals to detect a system state. The specific system can be, for example, a living human, where the signals are biological parameters. A processor operates the data points from two signals by defining respective classes of contiguous runs of a prescribed length of the data. The processor then assigns quantitative values to measure the regularity and stability of similar patterns between the first and second sets of classes of data points from the defined classes. These assigned quantitative values are aggregated to quantify the degree of asynchrony or conditional irregularity between the signals. Claims The invention claimed is: Description BACKGROUND OF THE INVENTION
TABLE 1a
______________________________________
TABLE OF RESULTS (HEALTHY)
HEALTH ApEn VAR WT(GMS.) AGE(WKS)
SEX
______________________________________
1. HEALTHY 0.94 5.57 2050 36 M
2. HEALTHY 1.08 6.99 1750 33 F
3. HEALTHY 1.11 6.69 2010 31 M
4. HEALTHY 1.12 10.29
1890 33 F
5. HEALTHY 1.16 8.13 1800 34 F
6. HEALTHY 1.20 9.42 550 24 F
7. HEALTHY 1.24 8.53 1820 37 F
8. HEALTHY 1.25 17.65
2020 41 M
9. HEALTHY 1.27 8.56 3650 40 M
10. HEALTHY 1.27 11.08
1300 34 F
11. HEALTHY 1.29 11.95
1600 36 F
12. HEALTHY 1.30 10.31
1730 33 F
13. HEALTHY 1.30 9.54 3490 40 F
14. HEALTHY 1.38 14.31
3100 40 M
15. HEALTHY 1.40 15.10
4360 42 M
______________________________________
SUMMARY STATISTICS: (MEAN +/- SD)
ApEn: 1.22 +/- 0.12
VAR: 10.27 +/- 3.33
WT: 2210 +/- 1000
AGE: 35.6 +/- 4.7
Hormone Secretion Study A study was performed to examine the potential applicability of ApEn to clinical endocrinology, and to quantify pulsatility in hormone secretion data. The study evaluated the role of ApEn as a complementary statistic to widely employed pulse detection algorithms, represented herein by ULTRA (Van Cauter, E., "Quantitative Methods for the Analysis of Circadian and Episodic Hormone Fluctuations," In Human Pituitary Hormones: Circadian and Episodic Variations, edited by E. Van Cauter and G. Copinschi, The Hague:Martinus Nyhoff, 1981:1-25), via the analysis of two different classes of models that generate episodic data. ApEn is able to discern subtle system changes and to provide insights separate from those given by ULTRA. ApEn evaluates subordinate as well as peak behavior, and often provides a direct measure of feedback between subsystems. ApEn generally can distinguish systems given 180 data points and an intraassay coefficient of variation of 8%. Additionally, the models and the extant clinical data are both consistent with episodic, not periodic, normative physiology. Thus, approximate entropy (ApEn), as a statistic, is applicable to hormone secretion data. Given the presence of a non-trivial amount of noise, there are two steps in performing hormone secretion pulse analysis. The first is separating the "true" secretion time-series from the noise. The second step is in evaluating the resulting "true" time-series. While these two steps are typically commingled in each algorithm, this is a complementarity between ApEn and the pulse-identification algorithm, due to their different approaches to the second step. ApEn summarizes the time-series by a single number, whereas the pulse-identification algorithms identify peak occurrences and amplitudes. ApEn discerns changes in underlying episodic behavior that do not reflect in changes in peak occurrences or amplitudes, while the pulse-identification algorithms ignore such information. Implicit to current models of hormone release is a periodicity assumption, with deviations attributed to noise. Two models which are capable of generating by themselves episodic, but not periodic, data are presented herein. In each model, there are several parameters that are varied, to generate a variety of data sets. For each model, ability of ApEn and a widely-used pulse-identification algorithm, ULTRA, to distinguish among the data sets generated by these models is evaluated. It is not suggested that these models represent known physiological systems, but rather these are offered as representative of alternative hypotheses to be considered when explaining observed episodic hormonal secretion. The present focus is not to propose a model that best mimics physiological reality, but rather to propose a new use of a statistic that gives different insights than are given by pulse-counting algorithms. Episodic Hormone Secretion Episodic, or pulsatile, secretion of hormones is an increasingly general finding in endocrinology. With the availability of sensitive radioimmunoassays (RIAs), which require only small sample volumes, protocols employing frequent sampling became possible. Furthermore, methods which help distinguish assay noise from biological variability make pulse detection a more rigorous endeavor. Studies employing such techniques in humans and diverse animal species have characterized pulsatile secretion of a large number of hormones, including luteinizing hormone (LH), insulin, progesterone, glucagon, growth hormone, ACTH, cortisol, prolactin, aldosterone, and HCG. Elucidating the secretory patterns of hormone release has not only shed light on endocrine physiology, but also clarified the pathophysiology and improved the treatment of some diseases. For example, derangement in the episodic secretion of LH underlies some common disorders in humans, such as polycystic ovary syndrome, and hypogonadotropic hypogonadism. Administration of LHRH in a periodic fashion, designed to produce a normal LH secretory pattern, improved the pharmacologic therapy of these disorders. Similarly, elucidation of pulsatile insulin secretion in normal subjects laid the groundwork for the discovery of abnormal insulin secretory patterns on diabetic, and improved the efficacy of insulin replacement therapy by administration of the hormone in a periodic fashion. Current Pulse-Identification Algorithms The tools currently employed by endocrinologists to analyze the pulsatility of hormone secretion data fall under the aegis of peak-identification algorithms. The philosophy of these methods is to identify the "true" peaks in the data, distinct from apparent peaks generated by the random variations due to assay imprecision. Once these true peaks are identified, one may be able to determine normal and abnormal ranges of pulse frequency, amplitude, and duration, and hence potentially identify abnormal secretion. There are considerable differences among the algorithms, due to a variety of approaches in handling the intraassay noise. This intraassay variation typically has a coefficient of variation (CV) of between 6% and 14% (e.g., Fuchs, A. R., K. Goeschen, and P. Husslein, "Oxytocin and the Initiation of Human Parturition III: Plasma Concentration of Oxytocin and 13, 14-dihydro-15 Keto-prosaglandin F2-alpha in Spontaneous and Oxytocin-Induced Labor at Term," Am. J. Obstet. Gynecol. 147 (1983):497-502), an amount of noise that can in some instances make true peak detection very difficult. Nonetheless, for all of these algorithms, in the absence of noise, (i) one achieves identical peak detection, and (ii) changes in subordinate patterns that do not result in new or altered peaks are ignored. The following eight pulse-detection programs are among those most widely available and extensively employed: Santen and Bardin (Santen, R. J., and C. W. Bardin, "Episodic Luteinizing Hormone Secretion in Man: Pulse Analysis, Clinical Interpretation, Physiologic Mechanisms," J. Clin. Invest. 52 (1973):2617-2628); modified Santen and Bardin; ULTRA; PULSAR (Merriam, G. R., and K. W. Wachter, "Algorithms for the Study of Episodic Hormone Secretion," Am. J. Physiol. 243 (1982):E310-318); Cycle Detector (Clifton, D. K., and R. A. Steiner, "Cycle Detection: A Technique for Estimating the Frequency and Amplitude of Episodic Fluctuations in Blood Hormone and Substrate Concentrations," Endocrinology 112 (1983): 1057-1064), Regional Dual Threshold (Veldhuis, J. D., J. Weiss, N. Mauras, A. D. Rogol, W. S. Evans, and M. L. Johnson, "Appraising Endocrine Pulse Signals at Low Circulating Hormone Concentrations: Use of Regional Coefficients of Variation in the Experimental Series to Analyze Pulsatile Luteinizing Hormone Release. Pediatr. Res. 20 (1986):632-637), Cluster (Veldhuis, J. D., and M. L. Johnson, "Cluster Analysis: A Simple, Versatile, and Robust Algorithm for Endocrine Pulse Detection," Am. J. Physiol. 250 (1986):E486-493); and Detect (Oerter, K. E., V. Guardabasso, and D. Rodbard, "Detection and Characterization of Peaks and Estimation of Instantaneous Secretory Rate for Episodic Pulsatile Hormone Secretion," Comput. Biomed. Res. 19 (1986):170-191). The similarity of the pulse-identification algorithms in the presence of negligible noise, the apparent relative robustness to non-trivial CVs, the usefulness with 50 to 200 data points, and the philosophy of peak analysis as the means to evaluate pulsatility bond this class of algorithms together. ULTRA has been chosen as representative of these algorithms in performing the comparisons with ApEn below. It is expected that another choice of pulse-detection algorithm, for the purpose of comparison with ApEn, would give quite similar results. Based on published time-series of hormonal concentration levels, there is the need for an added dimension in the analysis of episodic hormone release, beyond monitoring the pulse count and related statistics. Lang et al., (Lang, D. A., D. R. Matthews, and R. C. Turner, "Brief, Irregular Oscillations of Basal Plasma Insulin and Glucose Concentrations in Diabetic Men," Diabetes 30 (1981):435-439) conclude that brief, irregular oscillations in plasma insulin levels, in maturity-onset diabetics, are superimposed on longer term oscillatory fluctuations commonly observed in the non-diabetic. ApEn provides a quantification of the regularity of these data, which is useful for distinguishing a diabetic's insulin secretion patterns from those of a non-diabetic. Furthermore, episodic variation in hormones often has revealed complex patterns, challenging existing programs to characterize, and then differentiate, a "diseased" pattern from a healthy one. Finally, frequency distributions of discrete LH pulse properties, given by Urban (Urban, R. J., W. S. Evans, A. D. Rogol, D. L. Kaiser, M. L. Johnson, and J. D. Veldhuis, "Contemporary Aspects of Discrete Peak-Detection Algorithms: I. The paradigm of the Luteinizing Hormone Pulse Signal in Men," Endocrine Revs. 9 (1988):3-37) and based on nearly 200 pulses, show significantly non-Gaussian distributions for both pulse frequencies and amplitudes. The asymmetry of these distribution is not consonant with the typical assumption of periodic pulses in the presence of symmetrically distributed noise. One thus either concludes a lack of periodicity in these LH pulses, or at least must entertain the possibility of such a periodicity in constructing algorithms to analyze such series. The crucial difficulty in applying conventional entropy measurements to hormone secretion data is that hormone secretion data are relatively few in number (at most, several hundred data points), whereas an accurate conventional entropy calculation for an underlying system of dimension d typically requires from 10.sup.d to 30.sup.d data points (Wolf, A., J. B. Swift, H. L. Swinney, and J. A. Vastano, "Determining Lyapunov Exponents from a Time-Series," Physica 16D (1985):285-317). The number of data points is key, because there is no reason to anticipate, and no evidence to show, that data typically encountered from such complex interacting systems of glands and hormones that form endocrine systems be low-dimensional. Furthermore, one cannot assume that hormonal secretion is correctly modelled by deterministic chaos, as opposed to a stochastic model. ApEn has many of the characteristics deemed important for effective characterization of episodic hormone release as described by Urban et al. ApEn is objective, simple to use, via existing FORTRAN and C-language computer programs, and ApEn produces a single number. ApEn has minimal dependence on the specific type of signal or noise present in the underlying data. ApEn is versatile; it can be used for any time-series data analysis, to compute a measure of regularity. For hormone pulse detection, ApEn is readily adaptable to differences in sampling frequency and duration, assay performance, and signal-noise ratios. ApEn is very stable to small changes in noise characteristics, infrequent and significant data artifacts, and changes in sampling frequency. ApEn is concordant with visual inspection. ApEn accounts for a variety of dominant and subordinate patterns in data; notably, ApEn is affected by changes in underlying episodic behavior that do not reflect in changes in peak occurrences or amplitudes. Additionally, ApEn provides a direct barometer of feedback system change in some coupled systems. Thus, ApEn is useful in shedding insight into interactions among hormones, indicating a source of underlying physiologic deviations, such as a breakdown in the normal system feedback process. Model Comparison Framework Results from ApEn and ULTRA calculations for test data from two models are discussed below. To calculate ApEn and ULTRA for these data, certain inputs in each algorithm must be specified. For ApEn, m is set to 2 throughout the models, and r is chosen, fixed for each model, to equal about 20% of the standard deviation of a typical data set. This results in choices of r=0.4 for the Ueda model and r=0.1 for the Rossler model, consistent with guidelines given by Pincus (Pincus, S. M., I. M. Gladstone, and R. A. Ehrenkranz, "A Regularity Statistic for Medical Data Analysis," J. Clin. Monit. 7(4) (October 1991):335-345) For ULTRA, 3 CVs are chosen as the threshold for the Ueda model, and 2 CVs are chosen as the threshold for the Rossler model. This is consistent with Van Cauter's guidelines (applied to the "predominant" pattern in each instance). To determine concentration ranges and CV values for each range, one works backwards from the noise standard deviation data given in each version of each model. In each version, there is a model in which Gaussian noise of a fixed standard deviation (sdev) is superimposed on all the data, to model the inaccuracy of assay measurements. The output concentration ranges, for each time series, is divided into 8 pieces, each with a mean m. For each range, the CV is set to be sdev/m. The output of ApEn is a number, while the output of ULTRA is an identified set of peaks of pulses in the data. From each ULTRA output, the number of pulses, the average and standard deviation of pulse lengths, and the average and standard deviation of pulse amplitudes are calculated. For each model, a table is used to summarize the runs. Each line in the tables lists the run number, number of data points in the time series, and input model characteristics: parameter choices, and standard deviation of superimposed Gaussian noise. The output data includes the mean and standard deviation of the time series, ApEn value, number of pulses, and mean and standard deviation for both pulse frequency and pulse amplitude. Ueda Differential Equation Model The equation x+Ax+x.sup.3 =Bcost (5) is a differential equation that has received considerable attention in recent years, due in great part to studies by Ueda (Ueda, Y., "Steady Motions Exhibited by Duffing's Equation: A Picture Book of Regular and Chaotic Motions," In New Approaches to Nonlinear Problems in Dynamics, edited by P. J. Holmes, Philadelphia, PA.:SIAM, 1980:311-322) showing that the long-term dynamics of the solution represent steady-state chaotic behavior, for parameter values A=0.05, B=7.5. This equation, where the dots denote differentiation with respect to time t, describes the behavior of the variable x over time; for each time, the corresponding value of x can be calculated (via numerical methods), to deduce a time-plot of x as illustrated in FIGS. 7A-C. Equation 5 may be used in mechanical engineering, e.g., to model the motion of a sinusoidally forced structure undergoing large elastic deflections. The solution is bounded, episodic, yet nonperiodic. Here, Equation 5 is analyzed for five (A,B) pairs: (0.05, 7.5), (0.05, 8.5), (0.05, 12.0), (0.09, 7.5), and (0.21, 7.5). For each pair, equation 5 is solved as a function in time by an explicit time step method, .DELTA.t=0.002. A time series is extracted from the solution by sampling every 0.5 t-units. This sampling rate was chosen to yield about 12 data points per episode, and generates the baseline series. This is consistent with Yates, (Yates, F. E., "Analysis of Endocrine Signals: The Engineering and Physics of Biochemical Communication Systems", Bio. Reprod. 24 (1981):73-94), where samples of at least 6 times the expected frequency are seen as necessary to deduce periodicities, and with Veldhuis (Veldhuis, J. D., W. S. Evans, A. D. Rogol, C. R. Drake, M. O. Thorner, G. R. Merriam, and M. L. Johnson, "Intensified Rates of Venous Sampling Unmask the Presence of Spontaneous, High Frequency Pulsation of LH in Men," J. Clin. Endocrinol. Metab. 59 (1984):96-102), which notes the clinical need for intensified sampling rates. The solution time-series is post-processed by converting x to x+6.0 for all data values. This is done to ensure positive values, in the range 3.0 to 9.0, to mimic endocrine data. Uniform white noise is added to each baseline value to deduce the final series. For each pair, two different lengths of series are analyzed, 180 points and 900 points. For (0.05, 7.5) and (0.05, 8.5), the series is analyzed with 2000 points. This model is analyzed for two primary reasons. First, the model exemplifies a simple system which gives rise to highly nontrivial, putatively pulsatile behavior. Second, the model forces a careful examination of the meaning of pulsatility, to ensure that the quantitative tools used reasonably correspond to intuitive expectations. The crucial property of the solution to the Ueda equation is that it is episodic, but truly non-periodic. The equation's recurrent nature is evidenced by the fact that certain patterns in the waveform repeat themselves at irregular intervals, but there is never exact repetition. There is an apparent baseline frequency per episode (pulses), though there is temporal variation of a non-periodic nature. Furthermore, there are second-order, irregularly varying wiggles in the episodes that are generated by the model itself. This system is an appropriate model for hormone secretion, with normal secretion given by a model with A and B as stated above fixed, with A=0.5, and B between 8.5 and 12.0. On the basis of time-series data alone, the system can detect that certain data came from an "abnormal" system for which A=0.5 and B=15.0. As reported by Table 2, runs 1-10, the pulse count for noiseless systems is given as half the number of sign changes. This property is common to many of the current pulse identification algorithms, in which a pulse is flagged as a measured rise and fall, with both the rise and fall indicated by some percentage rise and fall times the noise level.
TABLE 2
__________________________________________________________________________
ULTRA Statistics
Ueda
Run Parameters
Input Noise No. of Sign
No. of
Avg. Avg.
No No. of Points
K B SD Mean
SD ApEn
Changes
Pulses
Freq.
SD Freq.
Amp
SD Amp
__________________________________________________________________________
1 180 0.05
7.5
0.0 5.808
1.593
0.677
56 28 6.185
2.450
7.368
1.834
2 180 0.05
8.5
0.0 6.401
1.612
0.543
50 25 6.792
3.400
7.771
1.276
3 180 0.09
7.5
0.0 6.074
1.650
0.574
66 33 5.406
1.292
7.363
1.596
4 180 0.05
12.0
0.0 5.992
1.781
0.762
81 40 4.410
1.044
7.406
1.692
5 180 0.21
7.5
0.0 6.174
1.549
0.676
46 23 7.818
3.390
7.624
1.338
6 900 0.05
7.5
0.0 5.973
1.597
0.894
275 137 6.463
2.839
7.476
1.576
7 900 0.05
8.5
0.0 6.533
1.567
0.466
213 105 8.452
4.031
8.028
0.769
8 900 0.09
7.5
0.0 6.034
1.637
0.590
333 166 5.388
1.447
7.335
1.570
9 900 0.05
12.0
0.0 6.068
1.758
1.153
401 200 4.462
1.131
7.531
1.540
10 900 0.21
7.5
0.0 5.908
1.550
0.666
227 112 7.991
3.361
7.635
1.550
11 900 0.05
7.5
0.05 5.971
1.597
0.904
281 117 7.578
3.492
7.731
1.486
12 900 0.05
7.5
0.1 5.968
1.599
0.953
287 104 8.534
3.694
7.925
1.393
13 900 0.05
7.5
0.2 5.963
1.609
1.091
299 97 9.156
3.675
8.098
1.250
14 900 0.05
7.5
0.4 5.953
1.647
1.336
367 84 10.59
3.425
8.417
0.975
15 900 0.05
8.5
0.05 6.530
1.567
0.473
213 77 11.57
2.568
8.311
0.611
16 900 0.05
8.5
0.1 6.528
1.569
0.510
253 76 11.72
2.408
8.320
0.603
17 900 0.05
8.5
0.2 6.523
1.577
0.742
289 78 11.41
2.769
8.305
0.613
18 900 0.05
8.5
0.4 6.513
1.614
1.196
379 77 11.57
2.806
8.479
0.400
19 900 0.09
7.5
0.05 6.031
1.638
0.602
333 165 5.421
1.486
7.332
1.575
20 900 0.09
7.5
0.1 6.029
1.640
0.634
335 155 5.773
1.881
7.371
1.582
21 900 0.09
7.5
0.2 6.024
1.650
0.909
339 121 7.408
3.245
7.789
1.390
22 900 0.09
7.5
0.4 6.014
1.687
1.292
365 92 9.769
3.715
8.340
1.029
23 2,000 0.05
7.5
0.0 6.075
1.597
0.871
588 294 6.782
3.116
7.544
1.469
24 2,000 0.05
8.5
0.0 6.559
1.556
0.443
460 229 8.715
4.077
8.069
0.639
__________________________________________________________________________
In the absence of noise, any rise is considered a pulse ascent, and any fall considered a pulse descent. Therefore, to evaluate ULTRA as a pulse-counting algorithm, it suffices to examine the statistical properties of the algorithm that counts the number of sign changes. This statistic has been extensively examined (Sen, P. K., "Signed-Rank Statistics," In Encyclopedia of Statistical Sciences 8, edited by S. Kotz and N. L. Johnson. New York:John Wiley, 1988:461-466) and provides useful information. It does not, however, utilize any information contained in the magnitudes associated with the sign changes, so that a tiny wiggle counts as much as a large wave. An instance of an improved measure is given by the Wilcoxon signed-rank statistic, a standard non-parametric statistical test. In this context, ranks would be given to the sign changes, with the largest rank to the greatest sign change. Hence, big pulses "count" more than little pulses, possibly a desired characteristic in the goal to distinguish normal from abnormal behavior. A central issue for this model is apparent upon examination of FIGS. 7A-C. Time-series output is shown in FIGS. 7A-C for three pairs of parameter values, (a) K=0.05, B=7.5, (b) K=0.05, B=8.5, and (c) K=0.09, B=7.5, respectively. These series are apparently different, but quantitative tools to distinguish them are not a priori apparent. These series have mean approximately equal to 6, standard deviation approximately equal to 1.6. Each series has a "period" of 2.pi., but no two periods are identical; there are different peak amplitudes, shapes, and subordinate "wiggles" throughout. Both ApEn and ULTRA distinguish versions of this model, but the results require scrutiny, because they appear to be in disagreement. First, runs 1-10, are summarized in Table 2. They represent runs for the five K-B pairs specified above, for two series lengths, 180 points and 900 points. According to ApEn, these versions rank (from most random to least random, in descending order) as (0.05, 12.0), (0.05, 7.5), (0.21, 7.5), (0.09, 7.5), (0.05, 8.5). This order is maintained for both 180 and 900 points, although several distinctions are sharper for 900 points than for 180 points. For this model, 900 points yields good convergence for ApEn; comparing run 6 to run 23 (900 vs. 2000 points, K=0.05, B=7.5), ApEn changes from 0.894 to 0.871. Similarly, comparing run 7 to run 24, (900 vs. 2000 points, K=0.05, B=8.5), ApEn changes from 0.466 to 0.443. According to ULTRA, these versions rank (from most random to least random) as (0.05, 12.0), (0.09, 7.5), (0.05, 7.5), (0.21, 7.5), (0.05, 8.5). This order is nearly maintained for both 180 and 900 points, although the last two versions reverse order in the 180 and 900 point cases. Furthermore, with the exception of the (0.05, 8.5) case, a five-fold increase in point count corresponds to virtually a five-fold increase in pulse number. This ratio of pulses to points is maintained in the two 2000-point runs, hence the 900-point runs are sufficiently long to extract the salient pulse information here. However, there is an apparent conflict over which of (0.05, 7.5) or (0.09, 7.5) is more random (unpatterned). The Poincare section is a tool to resolve this impasse. First, a phase space plot is generated (for each series), plotting the trajectory of x versus its time derivative, dx/dt. To insure a sequence of strictly comparable points, the trajectory is marked stroboscopically at times that are an integer multiple of the forcing period 2.pi.. The resulting plot, in the x-dx/dt plane, shows only the strobed points as the Poincare section. If the motion of the system were strictly periodic with the frequency of the forcing, the strobe point would all be the same point, repeating indefinitely. If the true motion were multiply periodic, then a sequence of n dots would appear, repeated indefinitely. More complicated dynamics are represented by more filled out Poincare section portraits, which correspond to greater ApEn. It can be shown (data not shown) that FIG. 7A has the most complicated dynamics, FIG. 7C has the next most complicated dynamics, and FIG. 7B has the least complicated dynamics. This corresponds to a greatest randomness for (0.05, 7.5), then (0.09, 7.5), followed by (0.05, 8.5), the order given by ApEn. Furthermore, the respective ApEn values, 0.894, 0.590, and 0.466, seem to correspond to the intuition that the (0.09, 7.5) case is closer to the (0.05, 8.5) case in randomness than to the (0.05, 7.5) case. The apparent inconsistency in ULTRA is explained by its equal weighting of each of many tiny wiggles and the larger sign changes. The (0.09, 7.5) case has the greatest number of sign changes of the three cases examined in FIGS. 7A-C, but these sign changes, particularly the "small wiggles," tend to occur near similar locations in each major "pulse." This can be virtually expressed by areas of darker clustering in phase portraits. Greater randomness would be marked by a greater spread of these dark clusters. The last point reemphasizes the foibles of the sign change algorithm, as opposed to a weighted sign change algorithm. Returning to Table 2, runs 11-22 further illustrate the difficulties that these small wiggles pose for ULTRA. For each of the versions (0.05, 7.5), (0.05, 8.5), and (0.9, 7.5), four different noise levels, standard deviations of 0.05, 0.1, 0.2, and 0.4, corresponding to CVs of approximately 1%, 2%, 4%, and 8%. In the (0.05, 7.5) case, ULTRA noted 137 pulses at 0 noise, compared to 117 pulses at 0.05 noise, and 104 pulses at 0.1 noise. This represents a computational loss of about 15% of pulses at 1% CV. In the (0.05, 8.5) case, ULTRA noted 105 pulses at 0 noise, compared to roughly 77 pulses in the presence of at least 0.05 noise. These 77 pulses represent, almost solely, the large pulses of approximate duration 2.pi.. Virtually all the small wiggles were effectively ignored in the presence of the noise levels noted above. This represents a computational loss of about 27% of pulses at 1% CV. In the (0.09, 7.5) case, ULTRA behaved more robustly at low noise levels, with 166 pulses at 0 noise, 165 pulses at 0.05 noise, and 155 pulses at 0.1 noise. ApEn performs more robustly at low noise levels. In the (0.05, 7.5) case, ApEn is 0.894 at 0 noise, 0.904 at 0.05 noise, and 0.953 at 0.1 noise. In the (0.05, 8.5) case, ApEn is 0.466 at 0 noise, 0.473 at 0.05 noise, and 0.510 at 0.1 noise. In the (0.09, 7.5) case, ApEn is 0.590 at 0 noise, 0.602 at 0.05 noise, and 0.634 at 0.1 noise. These all represent about a 1% to 2% change at 1% CV, and a 7% to 9% change at 2% CV. At each noise level ApEn maintains the order of randomness of these versions, although system distinction is much less marked at 0.4 noise level, as shown in FIGS. 8A-C, at which ApEn values are 1.336 for the (0.05, 7.5) case, 1.196 for the (0.05, 8.5) case, and 1.292 for the (0.09, 7.5) case. ULTRA also maintains its order of ranking these versions, with pulse counts of 84, 77, and 92 in the same three cases at 0.4 noise standard deviation. It is not surprising that the distinctions among the versions are muddied at this noise level; some of the small wiggles in the base physiological cases are accentuated, some are eliminated, and some new small wiggles emerge with 0.4 level noise. From analysis of this model, in the presence of noise, ULTRA tends to smooth out the time-series data, in effect eliminating some small wiggles in the process. In some contexts, that may be desirable, but in instances such as this model, in which numerous small, subordinate pulses are present, ULTRA is discarding physiological information. Rossler Feedback Model The Rossler Feedback Model is a coupled system of three variables, represented by three ordinary differential equations. This is considered as a putative model for the male reproductive endocrine system, with variables the pituitary portal concentration of LHRH, and the serum concentrations of luteinizing hormone (LH) and testosterone (T). These concentrations are modelled by a coupled feedback system: the LHRH secretion rate is given as a function of the local concentrations of LH and serum testosterone. The LH secretion rate is given as a function of the concentration of LHRH, plus a rate proportional to its own concentration. The testosterone secretion rate is given as a rate proportional to its own concentration, plus a term proportional to the product of the LHRH and testosterone levels. This feedback system is represented as follows, with K to be specified: ##EQU4## For each time, and each value of K, the corresponding concentration levels are calculated by an explicit time step method, .DELTA.t=0.005. A time series is extracted from the solution by sampling every 0.5 t-units. For suitable choices of K, the solutions have many of the qualitative features seen in clinical endocrine data. Here each version is defined by a choice for K. Changes in K can be thought to mirror the intensity of interaction between testosterone and LHRH levels. This system is analyzed for coupling levels K=0.4, 0.7, 0.8, 0.9, and 1.0. All this is done in a post-transient setting, in which the first 90 t-units are omitted from consideration. The solution time-series is then "post-processed" as follows, to ensure positive values: convert LHRH to 0.1(LHRH)+3.0, convert LH to 0.1(LH)+3.0, and T to 0.1(T)+3.0. Add white noise to each baseline value, for each of LHRH, LH, and T to deduce the time-series solution to the coupled system. For each noiseless version, two different lengths of series, 180 points and 900 points are analyzed. For the 900 point series, analyze three different versions of the model, K=0.7, 0.8, and 0.9, each under four different noise levels, noise standard deviations of 0.02, 0.05, 0.1, and 0.2. For K=0.8 and K=1.0, also analyze the series with 2000 points. This model is similar to one examined by Rossler (Rossler, O. E., "An Equation for Continuous Chaos," Phys. Lett. 57A (1976):397-398) as an example of a system that produced chaotic behavior for certain parameter values of K. It is also thematically similar to models by Smith, which are meant to plausibly model the male endocrine system, and are shown to capture some of the essential physiological dynamics of the true reproductive system. The Rossler model was analyzed, rather than the Smith model, for pedagogic reasons: distinctions among versions are sharper for the Rossler model than for the Smith model, though qualitatively quite similar. In any case, this model was analyzed for some of the reasons given by Smith. Relatively simple versions of this system can explain a number of possible qualitative modes of hormonal dynamics: serum concentrations that are constant in time, periodic in time, or chaotic in time. Most importantly, different behavioral modes can result solely from changes in defining system parameters, or internal interactions among the system subcomponents, and need not be produced by and external, driving force. For example, the onset of puberty, in one version of Smith's model, is seen to be generated simply by an appropriate change in certain system parameters, without an external switch or component entering into the fray. Furthermore, the Rossler model is substantially different from the Ueda model. In particular, the Rossler model is a function of several variables, and is an explicit feedback system. Thus, it is possible that neither ApEn nor ULTRA may detect changes in the feedback (coupling) rate, as seen by varying K. The model analyzed here was chosen to give either periodic, multiply periodic, or chaotic output for the behavior of LH with time, depending on K. In general, with increasing K, there is increasing system complexity: the LH behavior evolves from periodic to multiply periodic to chaotic. FIGS. 9A-E illustrate LH time-series output for the five coupling parameters of the Rossler coupled differential equation model, K=0.4, 0.7, 0.8, 0.9, and 1.0 in a noiseless environment. Virtually an identical pulse count is apparent in each of these systems. For K=0.4, the system is strictly periodic, while for K=0.7, the system is "twice-periodic," with a higher pulse always followed by a smaller pulse. The system is "four-times periodic" for K=0.8 (high-low-highest-lowest), and chaotic for K=0.9, and K=1.0. In these last two instances, no pattern of multiple pulses forms a fundamental period of its own. The increase in system complexity with increasing K can be further confirmed by phase-space plots. Phase- space plots serve a similar purpose to Poincare sections, to geometrically capture complexity via an appropriate perspective on the data. In a phase-space plot, the trajectory of LHRH versus LH is plotted so each point represents a single "LHRH-LH" pair of values at a fixed instant. Increased complexity manifests itself in more complicated phase-space portraits, which here is with increasing K. If the motion of the LH-system were singly periodic, the portrait would be a simple closed curve. Multiple periodicity is shown by multiple loops in a closed curve. Chaotic behavior is not represented by closed curves. Fine system structure in these versions is apparent with phase- space portraits produced from much longer time-series input than considered here. ULTRA's evaluation of the respective noiseless model versions is considered in runs 1-10 of Table 3. Runs 1-5 and 6-10 are 180 and 900 points long, respectively, with each set of 5 runs arranged in increasing K. Runs 1-5 give either 15 or 16 pulses for each series, and runs 6-10 give between 77 and 82 pulses for each series, indicating little version distinction based on pulse count. For the other statistics there is a distinct difference between the K=0.4 case and the other four versions, all of which produce quite similar values. In runs 7-10, the average frequency ranges from 11.28 to 11.63, the standard deviation of the frequency ranges from 0.636 to 0.650, the average amplitude ranges from 3.50 to 3.56, and the standard deviation of the amplitude ranges from 0.130 to 0.175. Only the last of these statistics, the amplitude standard deviation, shows any spread among the four versions, and even for this statistic, the lowest value is achieved for K=0.9, with both K=0.8 and K=1.0 versions slightly higher. This conflicts with intuition, which suggests that a lowest value for each of these statistics should be for either the least or the most complex system.
TABLE 3
__________________________________________________________________________
ULTRA Statistics
Coupling
Run Parameters
Input Noise No. of Sign
No. of
Avg. Avg.
No No. of Points
K SD Mean
SD ApEn
Changes
Pulses
Freq.
SD Freq.
Amp
SD Amp
__________________________________________________________________________
1 180 0.4 0.0 2.774
0.692
0.165
32 16 10.93
0.258
3.764
0.015
2 180 0.7 0.0 2.871
0.529
0.266
32 15 11.26
0.775
3.540
0.227
3 180 0.8 0.0 2.882
0.499
0.442
31 15 11.29
0.727
3.512
0.187
4 180 0.9 0.0 2.902
0.476
0.472
30 15 11.50
0.650
3.513
0.142
5 180 1.0 0.0 2.918
0.453
0.489
30 15 11.64
0.497
3.498
0.138
6 900 0.4 0.0 2.772
0.689
0.165
164 82 10.91
0.283
3.765
0.014
7 900 0.7 0.0 2.872
0.527
0.262
159 79 11.28
0.643
3.560
0.175
8 900 0.8 0.0 2.887
0.495
0.431
157 78 11.42
0.636
3.535
0.160
9 900 0.9 0.0 2.897
0.475
0.495
155 77 11.55
0.641
3.528
0.130
10 900 1.0 0.0 2.910
0.454
0.510
154 77 11.63
0.650
3.500
0.139
11 900 0.7 0.02 2.872
0.529
0.323
159 79 11.28
0.662
3.560
0.178
12 900 0.7 0.05 2.872
0.532
0.633
161 79 11.28
0.804
3.563
0.184
13 900 0.7 0.2 2.874
0.571
1.453
354 83 10.72
2.251
3.650
0.242
14 900 0.7 0.1 2.873
0.540
1.112
227 79 11.27
1.101
3.587
0.196
15 900 0.8 0.02 2.887
0.495
0.503
157 78 11.42
0.695
3.537
0.162
16 900 0.8 0.05 2.888
0.497
0.761
161 78 11.42
0.950
3.548
0.168
17 900 0.8 0.2 2.890
0.535
1.544
365 83 10.72
2.229
3.636
0.269
18 900 0.8 0.1 2.888
0.505
1.167
213 79 11.27
1.429
3.583
0.178
19 900 0.9 0.02 2.897
0.474
0.565
155 77 11.57
0.680
3.527
0.132
20 900 0.9 0.05 2.898
0.475
0.822
167 77 11.57
0.854
3.533
0.141
21 900 0.9 0.2 2.900
0.505
1.503
395 84 10.59
2.833
3.606
0.273
22 900 0.9 0.1 2.898
0.478
1.187
247 80 11.13
1.957
3.544
0.191
23 2,000 0.8 0.0 2.888
0.495
0.430
349 174 11.43
0.639
3.537
0.156
24 2,000 1.0 0.0 2.908
0.453
0.505
343 171 11.62
0.652
3.498
0.143
__________________________________________________________________________
For runs 1-5, in increasing K, ApEn values are 0.165, 0.266, 0.442, 0.472, and 0.489, monotonically increasing with K. For the 900 point runs (6-10), corresponding ApEn values are 0.165, 0.262, 0.431, 0.495, and 0.510, again steadily increasing with K. With these longer runs, distinction is sharper between the K=0.8 case and the K=0.9 and K=1.0 cases. Furthermore, ApEn is (slightly) larger for the K=1.0 case than for the K=0.9 case, establishing system distinction despite the presence of chaos in both instances. In addition, the ApEn values remain nearly constant for run lengths greater than 900 points, as indicated by runs 23 and 24. For K=0.8, no noise, ApEn=0.431 with 900 points, while ApEn=0.430 with 2000 points; for K=1.0, no noise, ApEn=0.510 with 900 points, while ApEn=0.505 with 2000 points. Hence for this model, (i) ApEn distinguishes all the versions from each other; (ii) ApEn, via monotonic increase, directly verifies the growing complexity and increased feedback with increased K; and (iii) establishes points (i) and (ii) with no more than 180 points necessary. Runs 11-22 indicate the effects of noise on the ULTRA and ApEn computations. For each of the three versions K=0.7, K=0.8, and K=0.9, four different noise levels were examined, standard deviations of 0.02, 0.05, 0.1 and 0.2, corresponding to CVs of approximately 1%, 2%, 4%, and 8%, for 900 point time-series. ULTRA maintained its pulse count of roughly 80 total pulses throughout these runs, increasingly slightly with noise level of 0.2 to 83, 83, and 84 pulses for the K=0.7, K=0.8, and K=0.9 cases, respectively. As above, this provided little distinction among these three systems. At 0.02 noise, ApEn maintained increasing order with complexity (0.323 vs. 0.503 vs. 0.565); similarly for the 0.05 and 0.1 noise levels (0.633 vs. 0.761 vs. 0.822, 1.112 vs. 1.167 vs. 1.187). With the 0.1 level noise, the system distinctions were becoming blurred, and with 0.2 noise, the system distinctions were obliterated, especially in the K=0.8 vs. K=0.9 cases (1.453 vs. 1.544 vs. 1.503), in which complexity is slightly reversed (due to "realization" and finite sample size issues). This blurring is evident in phase portraits comparing the K=0.8 and K=0.9 cases at 0.2 noise (data not shown). FIGS. 10A-C compares the LHRH, LH, and T time-series from the K=1.0 version of this model, and raises an important issue. The LHRH and LH time-series are visually similar; both have 16 pulses, similar amplitudes and general pulse characteristics. It could be expected that these hormones belong to a single autonomous system. The behavior of T, however, is visually discordant with the behavior of LHRH and LH; there are 12 pulses, long stretches of flat tracings, spiked pulses, and three pulses that are much greater in amplitude than the others. Thus, dissimilar pulsatile characteristics of hormonal plasma concentrations do not eliminate the possibility that the hormones may be derived from a single system, with no external influences. Discussion and General Conclusions Some general conclusions from the above runs can be inferred. ApEn and ULTRA provide different and complementary information from the data. ULTRA gives a first-order measure of the pulsatility of the system, via the pulse count and related statistics. ULTRA can be applied to data with 10% CV with 180 data points, typical values for current studies. For some systems, such as those defined by the Ueda and Rossler models above, ULTRA is relatively ineffective at distinguishing distinct versions of the systems, and may possibly give counterintuitive results. Subordinate pulses create difficulties for ULTRA, as do models in which pulse timing is reasonably constant, where the variation is in the patterned versus random behavior of the respective pulse amplitudes. The first-order, as opposed to finely-tuned behavior of ULTRA is further evidenced by the observation that in noiseless systems, ULTRA is statistically equivalent to a sign-change identifying algorithm. This algorithm was noted earlier to be useful, but to lack the greater versatility that appropriately weighted versions maintain. In contrast, at low intra-assay noise levels, with the stated input parameters ApEn can effectively distinguish all the distinct versions of each model from one another. In directly assessing the regularity of the data, ApEn can distinguish between versions of episodic behavior, as well as between episodic versus more random behavior. By considering all the time-series data, not just the data that make-up the pulse acmes, ApEn evaluates subordinate behavior. There is a significant increase in ApEn with increasing CV, though it is still possible to compare systems with identical intraassay CVs, even as high as 8%, via ApEn to discern system distinction. Such analyses produce ApEn values that are much larger than the corresponding values in noiseless systems; in a few cases, systems that are distinguished by ApEn at low CV are no longer distinguished at 8% CV. From the Ueda model, it is noted that there may be important regularity information in time-series data that can be effectively extracted only in the presence of a small intra-assay CV. For such purposes, ApEn is well suited, with a finer focus than that of the pulse-detection algorithms currently employed. The required decrease in intraassay CV from current levels is consistent with the direction in which endocrinologists are actively moving. To validate the above claim of effective distinction of model versions by ApEn, an estimate of ApEn standard deviation is determined. The estimate (Monte Carlo estimates, 100 replications per computed standard deviation) is provided for two quite different processes: the MIX(p) model introduced in Pincus (Pincus, S. M., "Approximate Entropy as a Measure of System Complexity," Proc. Natl. Acad. Sci. 88 (1991):2297-2301) and a paradigm for chaos, the parametrized logistic map, f.sub.a (x)=ax(1-x), 3.5<a<4.0. First define MIX(p):fix 0.ltoreq.p.ltoreq.1. Define X.sub.j =.sqroot.2 sin(2.pi.j/12) for all j, Y.sub.j =IID (Independent, Identically Distributed) uniform random variables on ›-.sqroot.3, .sqroot.3!, and Z.sub.j =IID random variables, Z.sub.j =1 with probability p, Z.sub.j =0 with probability 1-p. Then define MIX(p).sub.j =(1-Z.sub.j)X.sub.j +Z.sub.j Y.sub.j. This is a family of stochastic processes that samples a sine wave for p=0, is IID uniform for p=1, and intuitively becomes more "random" as p increases. For m=2, r=20% of the process standard deviation, and 900 points, the standard deviation of ApEn (MIX(p)), calculated for each 40 values of p equally spaced between 0 and 1, is less than 0.055 for all p. For 180 points, ApEn (same m and r) standard deviation is less than 0.07 for all p. For the logistic map, the "randomization" needed to make this deterministic map fit a Monte Carlo scenario is given by different choices for the initial condition. For m=2, r=20% of the process standard deviation, and 900 points, the standard deviation of ApEn (f.sub.a (x)), calculated for each of 50 values of "a" equally spaced between 3.5 and 4.0, is less than 0.015 for all a. For 180 points, ApEn (same m and r) standard deviation is less than 0.035 for all a. Thus ApEn values of a=1.1 and b=0.9 would have very high probability of coming from different processes, for either of these two model classes. The MIX process computation is appealing, in that the process is nearly IID (uncorrelated iterates) for p near 1. Because larger ApEn standard deviation generally corresponds to more uncorrelated processes, it is expected that the standard deviation bounds for ApEn for MIX(p) will provide bounds for a large class of deterministic and stochastic processes. Given the ApEn sensitivity to intra-assay CV, several caveats must be noted to ensure appropriate application of this method. If the same process is analyzed in two different laboratories, one with CV 2%, the other with CV 8%, the ApEn values can be significantly different. Also, if the same process is analyzed under two very different sampling regimens (e.g., samplings every 5 minutes, versus every 20 minutes), ApEn values can be quite different; in effect, the relative noise levels can be dissimilar. Thus, until CVs and other "noise levels" that vary from system to system are markedly reduced from present values, comparison of ApEn values should be restricted to data sets produced from similar settings (e.g., same laboratory and sampling frequency), which would ensure a relatively constant CV across samples. The comparisons done for the two models above, at a fixed CV level, model such a "homogeneously noise" environment, and as already noted, show valid ApEn distinction, given CVs at presently observed levels. Along the same lines, it is critical to distinguish between the comparison of ApEn (with fixed m and r) values for two data sets, given N data points, from the questions of convergence of ApEn for a specific system. The results from the two models analyzed above indicate that ApEn typically needed on the order of 900 points for convergence. In comparing systems with 180 data samples, ApEn distinguished most systems that were distinguished with 900 points, occasionally less sharply. Thus, a fixed sample length should be used for all data sets under study. The models analyzed above were chosen to illustrate different types of physiologically plausible behavior, and while there was no substantial effort to model a particular endocrine system, it would seem likely that a true endocrine system would be at least as mathematically complex as either of these models. Thus it is imperative that statistics, meant to evaluate pulses generated by true endocrine system hormones, be capable of effective discrimination of versions of the above models. A key observation from these models is that nonlinear systems can produce highly nontrivial, episodic, yet non-periodic output behavior from equations that are simple in appearance. Output that appears as a sequence of identical, sine wave-like pulses is usually associated with uncoupled, linear systems. Such linear systems have been extensively studied because they readily yield exact, analytic mathematical solutions. There is no a priori reason to anticipate that true endocrine systems be either linear or devoid of feedback. Hence, the likelihood that episodicity (no exactly repeating patterns) is physiologically normative must be considered. In addition to those considered above, stochastic models, such as Markov processes and networks of queues, could have been analyzed. Similar qualitative conclusions to those realized herein are anticipated. In complex systems of glands and hormones, a direct barometer of feedback, or interaction between systems would likely be insightful. Either a breakdown in or an excessive amount of feedback may mark the onset of disease, and a method that could directly mark such a change in feedback has added value. For the Rossler model, as the coupling parameter K was increased, ApEn steadily increased, thus providing a direct measure of increasing system complexity. In general, ApEn appears to increase with greater system coupling and greater attendant complexity. While coupled systems currently must be individually analyzed to ensure this increase of ApEn with feedback parameter, this property holds significant potential utility in practical applications. Above, potential near-term applicability was indicated, by observing that with 180 points, or with 8% CV, ApEn still was useful in drawing distinctions between most model versions. In a preferred embodiment of the invention, a randomized version of ApEn is applied to hormone level data. This randomized version of ApEn has the advantage that it can be coupled with bootstrapping methods (Efron, B., The Jacknife, the Bootstrap, and Other Resampling Plans, Philadelphia:SIAM, 1982:27-36) to yield a statistic that distinguishes data sets of 100 points with high probability (via a small variance), in the presence of nontrivial noise. Hence, greater applicability of ApEn to hormone level data can be achieved both by more accurate and numerous clinical data, and by statistical advances outside the clinical setting. In summary, the potential use of approximate entropoxy (ApEn) to quantify regularity in endocrine hormone data has been described. ApEn offers new insights in the detection of abnormal behavior, especially given modest increases in the number of data samples and in the accuracy of the serum concentration level at each sampling. Turbulence Measurement and Flow Control When a fluid impinges on an object, the undisturbed fluid pressure and the velocity of the fluid changes. Depending on the shape of the object, a wake may be formed, which sheds eddies. The eddies may be aperiodic or periodic. The formation of wakes is dependent on the Reynolds number, which is a dimensionless ratio of the inertial force to the viscous force of the fluid. An object in a fluid stream may be subject to the downstream shedding of vortices from alternating sides of the upstream object. As the wake frequency approaches the natural frequency of the structure, the periodic lift force increases asymptotically in magnitude. When resonance occurs, the structure fails. The neglection of this phenomenon has accounted for failures of numerous structures, including electric transmission lines, smokestacks, and bridges. Turbulence also affects the amount of fraction or drag between the object and the medium. As the fluid flow transitions from laminar to turbulent, the coefficient of drag increases. Increased drag results in inefficient flow of the medium past the object. The inefficiency caused by turbulence requires that additional energy be exerted to maintain the flow of the medium. For example, a vehicle (or vessel) in motion consumes more fuel when the air flow (or water flow) in the wake is turbulent instead of laminar. Hence, it is desirable to maintain laminar flow as long as possible. A preferred embodiment of the invention uses a negative feedback system to create maximum or minimum turbulence of a fluid flowing around a primary solid. Classical control or optimum control techniques are used. A secondary solid, smaller than the primary solid is placed in the fluid in such a fashion as to either encourage or discourage turbulence. The turbulence is controlled by critically pulsing, shaping, slowing or otherwise metering the fluid. The system can be adjusted for a small amount of turbulence that minimizes stress on the surface of the primary solid while maximizing flow, or any other complex combination of variables with a desired result. FIG. 11 is a schematic block diagram of a preferred embodiment of a turbulence measurement and flow control device according to the present invention. A medium 110, such as a fluid or a gas, is shown flowing across or through a region constrained by a primary solid 120. The medium 110 includes but is not limited to air, water, and blood. The primary solid 120 may be an airfoil or a hydrofoil (e.g. a wing, a propeller blade, or a rudder), a valve, a tube, a pipe, a channel, or any other structure that partially interferes with the flow of the medium 110. In particular, the primary solid 120 may be an artificial heart valve. The flow parameters of the medium 110 are measured by at least one sensor 130. The sensors 130 detect and quantify parameters such as speed, pressure, and direction of flow at specific locations in proximity to the primary solid 120. The measured parameters from sensors 130 are provided to a computational unit 140, which employs digital computations of approximate entropy to determine a time-varying parameter ApEn(t) for the medium 110 in proximity to the primary solid 120. The ApEn(t) parameter is continuously provided to a compensated negative feedback control 150. The system uses ApEn as a time-varying measure of turbulence, rather than the classical Reynolds Number, because ApEn is easier to measure, more immune to measurement noise and error, scale length independent, and completely shape independent. The feedback control 150 generates an optimum time-varying signal to provide to an actuator 160. The actuator 160 moves a secondary solid 170. The secondary solid 170 may be a constrictor, a flap, a vibrating plate, or any other structure that affects the flow of the medium 110, so as to change and optimize the flow characteristics of the medium 110 in proximity to the primary solid 120. In a stagnate medium, the approximate entropy wil | ||||||