INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

PRENATAL COCAINE EXPOSURE is believed to affect infant heart rate control and sleep-wake state organization. However, the number of studies that have investigated these effects are few, and the conclusions that have been drawn from the various data pools remain ambiguous (20). On one hand, some studies (21, 33) have found mean heart rate to be elevated and heart rate variability (HRV) to be reduced in infants with prenatal cocaine exposure. On the other hand, a recent study (23) found reduced median heart rates and increased HRV in cocaine-exposed neonates relative to controls. Spectral analysis of these data has shown that the higher HRV is due to increases in spectral power across all frequency bands in quiet sleep and increases in spectral power in the low-frequency (0.03-0.1 Hz) and mid-frequency (0.1-0.2 Hz) bands in active sleep (24). These results differ somewhat from the study of Oriol et al. (21), which found a significant reduction in high-frequency power.

The apparent discrepancy among these previous findings may be due in part to differences in the development of the subject groups that were studied as well as to the difficulty of controlling and determining the amount of prenatal cocaine exposure. Another possibility is that summary statistics and linear measures of heart rate dynamics may not be sufficiently sensitive to uncover differences in cardiovascular function between infants with prenatal cocaine exposure and controls without prior prenatal cocaine exposure. There may well be more subtle differences that become detectable only through methods that have the capability of characterizing the nonlinear aspects of the underlying dynamics. A number of studies (9, 11, 29) have suggested that the irregular behavior found in HRV may be a manifestation of deterministic chaos: complex dynamics that arise from nonlinear interactions among the many mechanisms that control or influence heart rate. However, other studies (7, 14), applying more stringent mathematical tests for chaos, found no evidence to support this hypothesis, although these researchers did find significant nonlinear correlations in the data. Recent findings by some groups (6, 13) suggest further that the irregular dynamics of HRV may be due in large part to stochastic influences.

In this study, we applied a variety of computational techniques to determine whether there are differences in the nonlinear dynamics of HRV of cocaine-exposed neonates and age-matched controls. We further hypothesized that the fluctuations in R-R interval (RRI) in both groups of neonates can be modeled as the output of a dynamic deterministic feedback system with stochastic noise as a possible perturbing input. Employing this structural framework, we sought to determine whether the dynamics of the feedback system could be adequately characterized by a linear model or whether it was necessary to include nonlinear contributions. Furthermore, this model allowed us to estimate the relative contributions of deterministic versus stochastic components to overall HRV in these two groups of infants during quiet and active sleep.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Subjects. We studied 12 normal neonates and 9 infants with prenatal cocaine exposure. All infants were products of single births by vaginal delivery with birth weights of >2,500 g and Apgar scores at 5 min of >7. The birth weights of the two groups of infants [controls 3,415 ± 130 g (SE) vs. cocaine-exposed 3,528 ± 236 g] were not significantly different. All subjects were studied at 2 wk postpartum. The mothers of all the subjects were single African-American and Hispanic women. Cocaine exposure was determined by maternal self-report or neonatal toxicological urinalysis (EMIT procedure) in the perinatal period. Histories of current and past cocaine (reported in terms of frequency of use during each trimester of pregnancy), alcohol (reported as absolute ounces of alcohol per week during each trimester), and tobacco (reported as mean number of cigarettes per week during each trimester) use were taken from each woman at the time of the sleep recording according to the recommendations of Day et al. (8). The Obstetric Complications Scale (17) was completed for each participant. Mothers of the control infants tested negative for cocaine and were selected from the same population of single African-American and Hispanic women. Radioimmunoassay of the mothers' hair (3) was used to verify cocaine exposure and the lack thereof in the cocaine group and the control group, respectively. The study was approved by the institutional review board of the King-Drew Medical Center (Los Angeles, CA), and each participating mother gave written informed consent.

Measurements and data preprocessing. Four-hour daytime recordings of the electrocardiogram (ECG) were obtained from the infants during spontaneous sleep and wakefulness. These recordings were made between 0900 and 1500 hours. Each 1-min epoch (of a total of ~240 epochs) was classified as quiet sleep, active sleep, indeterminate sleep, or waking on the basis of behavioral criteria using a previously reported protocol (25). The intervals between successive R waves of the ECG (RRI) were determined with 1-ms accuracy by subtracting the time of each R wave from the time of the previous R wave. In each subject, we selected for analysis two artifact-free segments of RRI data of ~1,000 beats (8-10 min) duration each; one of these was from quiet sleep and the other was from active sleep. Care was taken to ensure that no state changes occurred within a given segment. Segments representing wakefulness were not included in our comparisons because we were unable to find a sufficiently large number of data segments (containing 1,000 contiguous beats) that were free of artifacts. During wakefulness, there was frequently crying or other behavioral activities that produced motion artifacts in the ECG recordings.

Spectral analysis. A linear interpolation algorithm was first used to convert the RRI into equally spaced measures of heart rate with a new sampling rate of 16 Hz (4). After linear detrending, the power spectrum of HRV was computed from each data segment using the prewhitened autoregressive spectral analysis method of Birch et al. (5). Briefly, this procedure involved the following steps: fitting an autoregressive model to the data; determining the residuals between the measurements and the model predictions; computing the spectrum of the residuals via fast Fourier transform; and, finally, filtering the residuals spectrum with the autoregressive model to obtain the spectrum of variations in RRI. It should be noted that a relatively high resampling rate (16 Hz) was required to obtain good frequency resolution in the resulting spectrum, because this algorithm required the fast Fourier transform of the residuals. Before the spectral computations, a preliminary analysis showed that an autoregressive model of order 10 was adequate for fitting the data.

(1)

Approximate entropy. Approximate entropy (ApEn) is defined as the logarithmic likelihood that the patterns of the data that are close to each other will remain close for the next comparison with a longer pattern. Thus ApEn provides a generalized measure of regularity. A deterministic signal with high regularity has a higher probability of remaining close for longer vectors of the series and hence has a very small ApEn value. On the other hand, a random signal has a very low regularity and produces a high ApEn value.

(2)

(3a)

(3b)

Correlation dimension. The correlation dimension (CD) describes the dimensionality of the underlying process in relation to its geometrical reconstruction in phase space. We estimated CD using an approach based on the Grassberger-Procaccia algorithm (10). From the RRI time series, the following m-dimensional time-lag vectors were first constructed

(4)

where  is the embedding lag. For each of the vectors z_i we computed the distances  z_i  z_j to all the remaining vector points z_j, excluding those that were close because of temporal correlations (15). The number of data points within a distance r in the phase space for each vector was counted, and this is counted for all z_i. The correlation integral C(r) is given by

(5a)

where _min is the average correlation time (expressed in units of the number of data points), defined as the time taken for the autocorrelation function to first decay to 1/e. (u) is the Heaviside function, defined as

(5b)

C(r) was computed for a range of values of distances r, and, subsequently, log C(r) was plotted against log r. A scaling region was chosen in which this curve was approximately linear, and CD was computed as the slope of the curve in this region

(6)

The dimension m of the reconstructed vectors was selected by applying the method of false nearest neighbors (26). For the embedding time lag , we chose the value of one beat, because this was the natural time scale of the RRI time series (2). We also repeated our analyses with embedding lags of up to six beats (corresponding to the first zero crossing of the autocorrelation functions in most of our data sets), but these did not alter the relative differences of CD among subject groups and sleep states.

Nonlinear predictability. Nonlinear predictability provides a means for detecting determinism in any given time series. Predictability is low for a stochastic time series regardless of how far in the future one tries to predict. On the other hand, a periodic signal is highly predictive. With a chaotic signal or a correlated noise sequence, predictability would be high in the immediate future, but with increasing time steps, predictive capability would decrease significantly.

(7a)

(7b)

(7c)

Detection of nonlinearity using surrogate data. The preceding measures of nonlinear dynamics can be easily corrupted by the presence of stochastic noise. The surrogate data technique provides a means for testing the statistical significance of each computed measure. We employed the amplitude-adjusted Fourier transform algorithm (32) to generate surrogate data sets from the original time series. Here, the original data set was first rescaled so that the distribution became Gaussian. Surrogate data sets of the same length were then generated from this rescaled time series by randomizing the phase components of its Fourier transform while preserving the magnitude of the spectrum. Finally, the Gaussian surrogates were rescaled back to the original amplitude distribution of the data.

(8)

Parametric models. In addition to detecting nonlinearity in the underlying dynamics of HRV, we were also interested in determining the extent of the nonlinearity present. To quantify the degree of nonlinearity, we turned to parametric modeling. Figure 1 shows two possible model structures, both assuming delayed feedback, that can produce the dynamic fluctuations observed in the RRI time series. Model A assumes that the feedback dynamics are constrained to be linear, whereas in model B the dynamics of the feedback block are nonlinear. In model A, aperiodic fluctuations in RRI can only be produced when the feedback system is driven by a stochastic noise input. In contrast, in model B, aperiodic fluctuations in RRI can arise with or without any noise perturbation if the system dynamics are chaotic. The feedback structure inherent in both models encapsulates all the deterministic correlation between the present RRI and past changes in RRI. Thus both models represent the totality of all physiological mechanisms that can affect HRV, including respiration. The stochastic noise that enters these systems represents the combined effects of random fluctuations in autonomic neural modulation of heart rate, cardiac contractility, peripheral circulatory resistance, and blood pressure as well as transient variations in sleep-wake state (microarousals).

View larger version (11K): [in this window] [in a new window]   Fig. 1.   Parametric models for underlying dynamics of heart rate variability. Top: model A, linear model; bottom: model B, nonlinear model. RRIn, change in R-R interval (RRI) at beat n from the mean RRI of the data set in question.

(9)

(10)

(11)

(12)

(13)

Statistical analysis. Two-way repeated measures analysis of variance was employed, with subject group (control vs. cocaine) as one factor and sleep state (quiet sleep vs. active sleep) as the repeated factor. A Student-Newman-Keuls test was employed for post hoc multiple pairwise comparisons if statistical significance was indicated by the analysis of variance. All statistical tests were implemented using SigmaStat/Windows (Jandel Scientific; San Rafael, CA). The level of significance was set at P = 0.05 unless otherwise stated. In addition, each of the indexes of nonlinear dynamics (i.e., ApEn, CD, and nonlinear predictability) was tested in every subject for significance by computing  in Eq. 8 and determining whether the computed value was >2.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Spectral analysis. Table 1 shows the values for the estimated parameters in the two groups of infants in quiet and active sleep. Mean RRI in both groups of infants was not significantly different. There were no group differences in SDRR, LHR, and NHFP. However, NHFP was significantly higher in quiet sleep compared with active sleep in both groups (P < 0.001). At the same time, LHR and SDRR were each higher in active sleep relative to quiet sleep (P < 0.001 for both).

Approximate entropy. ApEn was not significantly different between the two groups regardless of which combinations of parameters (m = 2 or 3 and r = 10 or 20%) were employed. However, ApEn was higher in quiet sleep compared with active sleep (P < 0.001). The results for the individual subjects for m = 2 and r = 10% are presented graphically in Fig. 2 (cocaine subjects are shown as filled circles, whereas the controls appear as open squares). For each subject, the ApEn value deduced from data (horizontal axis) has been plotted against the mean ( ±2 confidence intervals, shown as error bars; vertical axis) of the corresponding ApEn values for the 10 sets of surrogates generated from the original time series. Significance levels for the difference in ApEn between the measured RRI and the corresponding surrogate data were >2 in virtually all subjects in both groups and states. In Fig. 2, this result is represented by the fact that almost all the circles and error bars lie above the line of identity, which is the graphical correlate of the null hypothesis. Thus ApEn values were found to be lower than what would have been expected if the dynamics of the RRI fluctuations reflected only linear correlations. In other words, the time series of all subjects showed a greater degree of regularity as a consequence of nonlinear correlations in the data.

View larger version (13K): [in this window] [in a new window]   Fig. 2.   Results of approximate entropy analysis. A: quiet sleep; B: active sleep.

Correlation dimension. The estimated values of CD for the two infant groups in the two sleep states are displayed in Table 1. Group differences were not significant. However, CD was higher in quiet sleep compared with active sleep in both groups (P < 0.001). Figure 3 shows the estimated CD values plotted against their corresponding surrogate data values. The surrogate data show higher CD values than the RRI signal itself in both states. The significance factor of this difference in the original data and the surrogate data was >2 for both groups, indicating the contribution of nonlinear correlations in the signal.

View larger version (12K): [in this window] [in a new window]   Fig. 3.   Results of correlation dimension analysis. A: quiet sleep; B: active sleep.

Prediction analysis. The normalized variance of the error between the one-step-ahead predictions and their corresponding data values was higher in quiet sleep versus active sleep (P < 0.001) but was not significantly different between the cocaine and control subjects (Table 1). Individual prediction errors are plotted against the corresponding prediction errors estimated from the surrogate data sets in Fig. 4. In most of the cases, prediction error is larger in the surrogate data relative to the prediction error deduced from the original data sets, implying again that the nonlinear contribution to the underlying dynamics was significant in both groups of infants in both sleep states.

View larger version (15K): [in this window] [in a new window]   Fig. 4.   Results of nonlinear prediction analysis shown as the normalized mean squared error. A: quiet sleep; B: active sleep.

Parametric model. The number of autoregressive terms that minimized the cost function (Eq. 12) for model A ranged from 7 to 14. In model B, the "optimum" model order contained 36-43 terms, which generally included linear, second-degree, and some third-degree terms. Figure 5 shows the percent CTV for model A (top) and model B (bottom) in each individual data set. Statistical analysis showed that there were no differences in CTV between subject groups. However, CTV increased substantially (P < 0.001) from model A to model B in both sleep states. Furthermore, there was a significant state versus model interaction (P < 0.001). In quiet sleep, the linear model (model A) was able to account for 31.6 ± 5.1% of the dynamic fluctuations in RRI in the cocaine group and 32.6 ± 5.2% in the control group, whereas in active sleep the linear CTV increased to 65.7 ± 3.3 and 65.2 ± 2.1%, respectively. In quiet sleep, CTV for the nonlinear model was 81.5 ± 4.7 and 76.6 ± 3.5% in the cocaine and control groups, respectively. These contributions remained little changed in active sleep at 85.6 ± 2.5 and 89.9 ± 1.0%, respectively. Because model B also includes linear terms, these results imply that the nonlinear contribution to RRI dynamics in quiet sleep was substantially larger than that in active sleep. However, there were no differences in nonlinear contributions between infant groups.

View larger version (16K): [in this window] [in a new window]   Fig. 5.   Results of parametric models shown as percent contribution to total variance. Top: model A, linear model; bottom: model B, nonlinear model.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Two findings emerged with great consistency in this study. First, there was a substantial amount of intersubject variability in all the measures of linear and nonlinear dynamics that were estimated. Second, we could find no significant differences in ApEn, CD, nonlinear predictability, or any of the spectral measures between a group of neonates with prenatal cocaine exposure and a group of age-matched controls. The second finding could well be a consequence of the first: a large degree of intersubject variability can easily mask small differences in heart rate dynamics between the groups. It should be emphasized that the large intersubject variability was found not only in the cocaine-exposed group but also in the control group. Thus it appears that any alterations in heart rate dynamics in cocaine-exposed neonates are likely to be too subtle for detection even by nonlinear techniques unless much larger sample sizes are employed. This could explain why previous studies using spectral analysis of HRV have arrived at differing conclusions. For instance, Mehta et al. (19) reported a significantly lower LHR and higher NHFP in 21 cocaine-exposed neonates, suggesting enhanced parasympathetic activity. Oriol et al. (21) found reduced overall HRV, LFP, and HFP in a group of five cocaine-exposed neonates relative to normal controls, which suggested increased sympathetic tone. On the other hand, they did not find significant differences in ApEn between the groups. Another confounding factor could be the possible influence of nonstationarity in these previous studies. In the Mehta et al. study (19), the spectral analysis results appear to have been based on data derived from Holter recordings of 22 h or more. Comparison of spectral indexes of HRV between the cocaine infants and controls was performed without consideration of the sleep-wake state. In the Oriol et al. study (21), the data segments analyzed were similar in length (~10 min) to ours. However, only visual inspection was employed to determine stationarity of the data, whereas in our present study we applied statistical testing to rule out any nonstationary behavior.

In our previous analysis of a larger group of 15 cocaine-exposed neonates and 13 controls (inclusive of the subjects studied here), the cocaine infants showed enhanced HRV, reflected by an increase in spectral power across all frequency bands during active sleep; however, in quiet sleep, only HFP was higher (24). The discrepancy between these results and our present findings may be due in part to the slightly larger sample size employed in the previous analysis. In the present study, we were constrained to use only a subset of the overall database because, in some of the subjects, we were not able to find contiguous periods of at least eight (1 min) epochs of quiet or active sleep. The discrepancy may also be related to methodological differences between our previous and present analyses. In the previous study, spectral estimates were computed on an epoch-by-epoch basis. These epochs were classified into one of the following four states: quiet sleep, active sleep, indeterminate sleep, and wakefulness. Subsequently, the median value deduced from all epochs in each state was taken to be representative of the spectral estimate in that state for a given subject. No attention was paid to whether these median values reflected heart rate dynamics during relatively stable and extended periods of a given sleep state. In contrast, in the present study, care was taken to ensure that the data analyzed were extracted from sections in which there were 8-10 contiguous epochs (1 min) of either quiet or active sleep; furthermore, we were careful to test these data segments for stationarity. Because sleep organization is known to be altered in cocaine-exposed infants (25), it is possible that the conclusions arrived at in our previous study may have been affected by the cardiorespiratory effects of transitions between states. If this was indeed the case, our present findings would suggest that the primary effect of prenatal cocaine exposure is the disorganization of sleep architecture and that any observed differences in heart rate control are secondary to these alterations in sleep-state patterning.

Application of the surrogate data method to our estimates of ApEn, CD, and nonlinear predictability showed that there was a significant nonlinear deterministic component in the HRV of both groups of neonates during quiet and active sleep. Parametric modeling confirmed this finding and, furthermore, allowed us to quantify the relative contributions of the linear and nonlinear components of the underlying dynamics. Our results suggest that the dynamic fluctuations in RRI in both cocaine-exposed and control infants can be modeled as the output of a deterministic nonlinear delayed-feedback system with stochastic noise as the perturbing input ("model B"). In both infant groups, the nonlinear system was of relatively low order, containing lagged products of RRI up to only the third degree. We found that the combined contributions from linear and nonlinear correlations accounted for between 65 and 99% of the total variance in the data, which meant that in some cases the direct contribution from stochastic noise was as low as 1%. This suggests that, in some of the data sets, the dynamics of HRV may have been chaotic. To test this possibility further, we estimated the largest Lyapunov exponents of these data sets using the method of Rosenstein et al. (27), with the presumption that the presence of a positive characteristic exponent would indicate chaos. We found that these exponents were positive but not significantly different from the Lyapunov exponents computed from the corresponding surrogate data. Thus the determination of whether chaos was present or not was inconclusive. This highlights a major limitation of current methods of nonlinear dynamical analysis: the sensitivity to noise of estimates derived from relatively short (<10,000 points) data sets.

One feature that was continually affirmed in all our computational tests was the clear difference in heart rate dynamics between quiet sleep and active sleep regardless of infant group. Overall HRV was higher in active sleep. NHFP was lower and LHR was higher in active sleep versus quiet sleep, indicating a shift in relative dominance of LFP versus HFP in HRV. All measures of nonlinearity, such as ApEn, CD, and nonlinear prediction error, decreased in active sleep, indicating a reduction in complexity and an increase in regularity in heart rate dynamics. Furthermore, parametric modeling showed that, in percent terms, the linear contribution to HRV was approximately doubled in active sleep, whereas the nonlinear contribution was reduced by roughly one-half. These results are consistent with previous findings that sympathetic modulation of heart rate is enhanced during active or rapid-eye movement sleep, whereas parasympathetic activity is decreased; conversely, in synchronized or quiet sleep, vagal modulation predominates (12, 18). Furthermore, sympathetic modulation leads to dynamic variations of heart rate that are predominantly periodic and in the low-to-mid-frequency range; in contrast, vagal modulation of heart rate produces broadband fluctuations that are dynamically more complex and much less predictable. There are a number of possible reasons why HRV assumes the form of low-frequency periodicities during sympathetic dominance in active sleep. First, the sinoatrial node can only track changes in sympathetic activity that are slower than 0.15 Hz, whereas vagal activity can modulate heart rate to much higher frequencies (28). Sympathetic modulation of changes in vascular resistance is also slow, on the order of several seconds. With increased sympathetic gain, these delays in the baroreflex control system can lead to oscillatory activity mediated by feedback instability. A less likely, but nevertheless plausible, explanation is that sympathetic dominance during active sleep leads to a filtering out of high-frequency activity, thereby unmasking the low-frequency oscillation that is intrinsic to central rhythmic modulation of neural activity (18).