Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Dev Psychobiol. Author manuscript; available in PMC 2010 December 1.
Published in final edited form as:
PMCID: PMC2936240

Heart Rate Variability in Response to Pain Stimulus in VLBW Infants Followed Longitudinally During NICU Stay


The objective of this longitudinal study, conducted in a neonatal intensive care unit, was to characterize the response to pain of high-risk very low birth weight infants (< 1500 g) from 23 through 38 weeks post-menstrual age (PMA) by measuring heart rate variability (HRV). Heart period data were recorded before, during, and after a heel lanced or wrist venipunctured blood draw for routine clinical evaluation. Pain response to the blood draw procedure and age-related changes of HRV in low-frequency and high-frequency bands were modeled with linear mixed-effects models. HRV in both bands decreased during pain, followed by a recovery to near-baseline levels. Venipuncture and mechanical ventilation were factors that attenuated the HRV response to pain. HRV at the baseline increased with post-menstrual age but the growth rate of high-frequency power was reduced in mechanically ventilated infants. There was some evidence that low-frequency HRV response to pain improved with advancing PMA.

Keywords: pain, respiratory sinus arrhythmia, vagal tone, very low birth weight infants, spectral analysis, mechanical ventilation, growth curve, heel lance, venipuncture, frequency domain measures


Painful procedures are suspected to contribute to poor neurodevelopmental outcomes in very low birth weight (VLBW) neonates (Grunau, 2002). Clinical scales have been developed to assess the behavioral expression of pain in preterm infants (e.g. Stevens, Johnston, Petryshen, & Taddio, 1996), but there is as yet no gold standard for measuring pain in these neonates. Behavioral responses to painful stimuli are less robust in VLBW neonates than in their full-term counterparts (Johnston, Stevens, Yang, & Horton, 1995; Grunau, Holsti, Haley, Oberlander, Weinberg, Solimano, Whitfield, Fitzgerald, & Yu, 2005), and may be more effectively used as a component of a pain scale based on multidimensional measures. Heart rate variability (HRV) provides measures of autonomic state and has been increasingly used as a measure of stress and of neonatal well-being (Porges, 1992; Rosenstock, Cassuto, & Zmora, 1999; Pumprla, Howorka, Groves, Chester, & Nolan, 2002), which suggests that it may be an informative measure of reactivity to painful events in ill preterm infants who are able to mount only nonspecific and ambiguous behavioral distress signals. Greater behavioral reactivity to various stimuli is associated with higher vagal tone, as measured by HRV, in infants and children (Arditi, Feldman, & Eidelman, 2006; Huffman, Bryan, del Carmen, Pederson, Doussard-Roosevelt, & Porges, 1998; Porges, Doussard-Roosevelt, Portales, & Suess, 1994; Fox, 1989) as well as in fetuses (Smith, Dmochowski, Muir, Kisilevsky, 2007).

Time and frequency domain measures of HRV have been studied in newborns during routine heel lancing. A clear stress response is provoked when the heel was pricked and squeezed (Johnston et al., 1995; McIntosh, van Veen, & Brameyer, 1993). The heart rate (HR) was observed to increase and low-frequency (LF) power decreased in response to pain (Lindh, Wiklund, Sandman, & Hakansson, 1997; Grunau, Oberlander, Whitfield, Fitzgerald, & Lee, 2001). Vagal tone has been reported to decrease during pain (Arditi et al., 2006). Cumulative exposure to pain has been linked to cortisol levels and facial responses, but not to HRV (Grunau et al., 2005). Long-term changes in the pain responses of formerly preterm neonates have been reported (Grunau, Whitfield, & Petrie, 1994; Porter, Grunau, & Anand, 1999; Oberlander, Grunau, Whitfield, Fitzgerald, Pitfield, & Saul, 2000; Buskila, Neumann, Zmora, Feldman, Bolotin, & Press, 2003). The cross-sectional studies reviewed above confound age and health status. The youngest infants tested also tend to be the sickest. Longitudinal studies have measured behavioral pain responses. Reduction in flexion reflex threshold with repeated heel lances has been reported (Fitzgerald, Millard, & McIntosh, 1989) and increasing maturational response to heelstick pain has been measured by behavioral indicators (Johnston & Stevens, 1996; Williams, Khattak, Garza, & Lasky, 2009). Infants in a NICU anticipated painful events with increased HR between baseline and the pickup of the leg at the beginning of a heelstick procedure (Goubet, Clifton, & Shah, 2001). It is evident that not all conclusions of the various studies can be reconciled with each other. Mechanical ventilation has been observed to induce respiratory sinus arrhythmia in preterm infants (van Ravenswaaij-Arts, Hopman, Kollee, Stoelinga, & van Geijn, 1995) but its impact on HRV response to pain is not known.

The objective of our study was to assess the HR and HRV responses to painful procedures in VLBW infants when followed longitudinally during their NICU stay. This is the first longitudinal study of HRV as a measure of pain in infants in a NICU setting. The effect of mechanical ventilation on post procedural HRV is unknown, and forms an important part of our considerations for this study. Differences were tracked between the pain response to two types of blood draw procedures used in this study, heel lance and venipuncture. As a post-hoc objective, we sought to characterize correlations between physiologic measures and behavioral pain responses.


Data Collection

This was a prospective cohort study following VLBW infants over several weeks in the NICU, ranging from 23 to 38 weeks post-menstrual age (PMA). Infants admitted to Children’s Memorial Hermann Hospital NICU <30 weeks PMA at birth were eligible for the study except those with major congenital anomalies (Table 1). Following parental informed consent, infants were enrolled in the study within 72 hours of birth and then longitudinally followed weekly through 30 weeks and biweekly thereafter. Only one infant dropped out during the study. Infants with major congenital anomalies were excluded due to short life expectancy that could interfere with the longitudinal design. Infants were studied between 2 to 6 am before, during, and after a heel lanced or wrist venipunctured blood draw for routine clinical evaluation while lying in a supine position. While most blood draws were done by heel lances, venipunctures were done in a few infants. The wrist was squeezed during the venipunctures to mimic heel lances as closely as possible.

Table 1
Patient Demographics

Infants were videotaped during the entire recording session to document behavioral responses to the painful stimulus. Heart rate and other vital sign data were recorded from the electrodes and transducers used for routine clinical monitoring in the NICU. Data collection began approximately 10 minutes before the painful procedure and continued for up to 10 minutes after the end of the painful procedure. The data were digitized and stored using PC-based custom instrumentation (Quantum Automation Inc., Houston, Texas). Heel warmers were used in the baseline period and skin was prepared with alcohol during the seconds before the start of the pain phase. Times for heel warming and skin preparation were similar across the study, while the time for heel squeeze varied depending on how fast blood was obtained. Clinical status and medications were recorded in order to assess their effect on HRV. The Institutional Review Board approved all studies.

There were 153 blood draw sessions from 38 VLBW infants. Precise timing of the blood draw was determined from the video. Thirty-nine of the datasets were excluded because the time of blood draw could not be determined with certainty from the videos, leaving 114 datasets from 35 VLBW neonates for the analysis. The duration of the painful procedure for these 114 datasets ranged between 102 and 762 seconds, with a median value of 324 seconds. The 114 datasets included 43 sets that were collected when infants were mechanically ventilated. A small fraction (10%) of blood draws were venipunctures through the wrist while the rest were heel lances with snap-loaded lancets. Due to the paucity of literature on this topic, realistic effect sizes could not be determined a priori; hence power calculations were not done for this study.

ECG Signal and Power Spectrum

The signal from the electrocardiograph (ECG) was processed by a National Instruments data acquisition module that digitized the signal at a rate of 1 kHz. Custom software programmed in LabVIEW (National Instruments Inc., Austin, Texas) identified the R-wave of the QRS complexes of the ECG signal (18) and generated a vector of interbeat intervals henceforth referred to as the RR-series. HRV is studied by analysis of the RR-series.

Natural variability of the heart period resulted in RR intervals that were non-uniformly sampled in time. Lomb’s algorithm (Lomb, 1976; Press, Teukolsky, Vetterling, & Flannery, 1992) is appropriate for computing the power spectrum for such non-uniformly sampled data sets since it avoids the low-pass effect inherent in uniform resampling techniques (Laguna & Moody, 1998; Chang, Monahan, Griffin, et al., 2001) The Lomb power spectrum was computed on segments of 128 consecutive heart periods, with 50% overlap between consecutive segments. A segment spanned a time interval of approximately 50-seconds, depending on each infant's average heart rate. Longer segment lengths yield more reliable estimates of low-frequency power, but also reduce the likelihood of achieving stationarity in the segment, which is discussed in more detail below.

Stationarity requires that the mean and variance of the time series, or more generally its probability distribution, be invariant in time (Rangayyan, 2002). Spectral measures are reliable only when computed on stationary data segments. Large errors in estimates of the power are introduced by artifacts such as missed QRS peaks, or false detections. Such errors may also be avoided by selecting segments that satisfy stationarity. We employed the Kolmogorov-Smirnov test of equivalence of distributions on two partitions of each segment (Shiavi, 1999). In our implementation of the stationarity test, the segment was partitioned into 2 subsets and the alpha value was set at 0.01. Testing showed that over 60% of segments passed the stationarity test in moderately noisy RR-series and that the mean spectral power in either low-frequency (LF) or high-frequency (HF) band was lower on segments that pass the stationarity test than in those that fail the test. Cubic polynomials were used to detrend segments before testing for stationarity and computing the spectral power.

The mean Nyquist frequency of the RR-series varied from segment-to-segment due to variation in mean heart rate. Since the Nyquist frequency is an upper limit on reliable information content of the power spectrum (CITE), it was required of each contributing RR-series segment that the average Nyquist frequency exceed the upper limit of the HF-power band. The power spectrum was averaged on all segments that passed the checks on stationarity and Nyquist frequency, and summed within a LF-power band (0.05 to 0.25 Hz) and a HF-power band (0.25 to 1.00 Hz) that is associated with respiratory sinus arrhythmia. This yielded estimates of the heart rate variability in timescales between 4 and 20 seconds (LF band) and between 1 and 4 seconds (HF band), approximately. The band-powers Pi were converted to decibels (dB) using a reference value of 0.02 Hz ms2: 10 log(Pi /0.02). This logarithmic conversion of units reduced the variability in the raw band powers and resulted in more normal distributions. This process was repeated for each data recording session for all infants in the study. All processing of the RR-series was done in MatLab 7.0/R14 (The Mathworks, Inc., Lowell, Mass.).

Repeated Measures

Repeated measures within a single blood draw session were termed phases and were treated as a categorical variable, while the post-menstrual age varied between repeated recordings and was treated as a continuous numerical variable. Four repeated measures of HR and HRV were extracted from each recording (Fig. 1), corresponding to four phases: baseline (3 to 1 minutes before blood draw), pain (0 to 2 minutes after blood draw), early recovery (2 to 4 minutes after the blood draw was completed), and late recovery (5 to 7 minutes after the blood draw was completed). Recovery was measured in two windows because it was not clear a priori which times should be used for monitoring HRV during recovery. Spectral powers within each band were averaged over all qualifying segments of RR-series that were contained within the phase-windows. It is assured that at least one segment of approximately 50-second duration contributes to a 2-minute phase. Since heel-squeezing too has an impact on HRV (Lindh, Wiklund, & Hakansson, 1999) the “end” of the painful procedure here was defined to be the end of the entire procedure, when the person administering the procedure departed from the vicinity of the infant. Heel lancing and squeezing occurred within seconds of each other and were treated as a single event in this study. Infants often responded with muscle tension and movement during the procedure, causing the RR-series to be noisy, especially during the blood draw. In such instances stationarity was violated and the data for a phase were not available. Of 456 repeated measures for the 114 data sets included in this study, 377 (83%) qualified for inclusion in the models presented below. For the baseline, pain, early recovery, and late recovery, the numbers of missing values were 17, 30, 11, and 21, respectively. Missing values were not replaced. The linear mixed models used for data analysis incorporate empirical Bayes estimation technique for handling missing data.

Figure 1
Example of an individual heart rate and spectral power of heart rate variability and the four phases or windows utilized for longitudinal study. Relative to the time of the pain stimulus, the baseline phase extended from −3 to −1 minutes ...

The very nature of this study meant that the heart rate and its variability were dynamically changing in time in response to the pain. The limited time width of the phases allowed us to analyze the response. Sensitivity of the results to the width of phase-windows was assessed by repeating the data analysis with 3-minute phase-windows instead of 2-minute phase-windows. Compared to the 2-minute phases defined previously, the wider baseline phase began 1-minute sooner and other phases ended 1-minute later. Wider phase-windows lead to fewer missing data points since there is a greater chance that the window contains at least one subinterval that satisfies the stationarity condition. The competing effect is that time resolution of the response to pain is lowered with wide windows. With 3-minute phase-windows the total number of missing values dropped from 17% to 13%; by phase the number of missing values were 10, 21, 7, 19 for the baseline, pain, early recovery, and late recovery, respectively. The wider baseline phase was between 4 and 1 minutes before the beginning of the blood draw, and the wider pain phase was during 0–3 minutes from the beginning of the painful procedure. The wider early and late recovery phases were from 2–5 and 5–8 minutes after the end of the painful procedure.

Linear Mixed Models

Age related changes in the LF and HF-power bands and HR were analyzed with linear mixed-effects models (Pinheiro & Bates, 2000; Diggle, Heagerty, Liang, & Zeger, 2002), with fixed-effects that were representative of the population and random-effects that varied between individuals. The longitudinal growth curve models were stratified by the type of blood draw, phase, and mechanical ventilation status. All models that are described here are for the case that had 2-minute phase-windows. Secondary models involving 3-minute phase-windows are briefly described in the Results section. The following model, computed with the method of maximum likelihood, was found to fit the data well for the LF power:


where yij is the HRV power of infant i at observation j, pij is the phase (baseline, during pain, early recovery, late recovery), vij is the type of blood draw (heel lance, venipuncture), zij is the mechanical ventilation status (no, yes), tij is the PMA that was centered at 30 weeks, β’s denote coefficients of fixed effects, b0i and b1i are random-effect intercept and slope of infant i, and εij is the error term. The reference level for measuring contrasts was set to heel lanced, non-ventilated, and baseline phase. The type of blood draw vij was bivalued, (0 when infants were heel lanced; 1 when venipunctured). The ventilation variable zij was also bivalued (0 when infants were not mechanically ventilated during a recording; 1 otherwise). The phase pij had four levels, so each term involving pij in the equation represents three terms (0 for the baseline phase; 1 for each of the other three phases).

Model exploration began with the full model involving 3-way interaction among the fixed effects, p * z * t, and was systematically reduced to retain essential terms only. Differences in pain response to heel lances and to venipuncture were modeled with an interaction term v * p. Polynomial dependences on t were explored up to rank five. A cubic dependence was found that served to modulate the linear growth rate at the ends of our observation range. The possibility of linear dependence of the power in recovery on the duration of the painful procedure was also investigated, but found to be unnecessary. The Aikaike Information Criterion (AIC) was used for model selection. The need for heterogenous variances in the strata and an autocorrelation structure nested within phase were also explored with AIC, and were found to be unnecessary. Model diagnostics satisfied linearity and homogeneity assumptions.

The procedures were repeated in constructing models for the HF-power and for the HR. The best model for HF-power was similar to that for the LF-power, but the terms involving the type of blood draw (v) and the interaction between age and phase (t * p) were not significant. Instead, the interaction term between age and ventilation (t * z) was present in the model for HF power:


Six outliers were removed from the models for HF power. Each of the six observations was confirmed to have arisen from RR-series with high-frequency noise and HF power values exceeded the box defined by 1.5 interquartile ranges above the 3rd quartile of the distribution. HR data was obtained from stationary segments of the RR-series and the following model was effective:


where Yij is the HR, b2i and b3i are random-effect coefficients corresponding to ventilation and interaction of age with ventilation, and other coefficients and variables are as described previously. Degrees-of-freedom per subject were more than sufficient to allow the additional random-effects. The AIC statistic indicated that within-subject autoregressive AR(1) structure stratified by phase was of benefit to the model. Plots of residuals indicated that linearity and homogeneity assumptions were met and there were no outliers. Statistical analysis was done in S-PLUS 7.0 (Insightful Corp., Seattle, Washington).

Behavioral Pain Response

Video recordings were analyzed for behavioral pain responses in a subset of blood draw sessions (n=21) for which infants’ faces and bodies were clearly visible in the video, and HR, HRV were available for all phases. The Neonatal Infant Pain Scale (NIPS) was utilized to score behavioral responses (Lawrence, Alcock, McGrath, Kay, MacMurray, & Dulberg, 1993). NIPS behaviors were scored in 1-minute intervals throughout each phase of the blood draw using Noldus Observer 5.0 (Noldus Information Technology, The Netherlands). The behaviors included facial expression (grimace), crying, respiration (irregular or regular), arm movement, leg movement, and state of arousal. Each component behavior was scored as the highest level reached by the newborn for each minute. The total NIPS score for each minute was the sum of individual behavior scores, and ranged from 0 to 7 (all behaviors except crying were scored as 0 or 1, crying was scored as 0, 1, or 2). Inter-rater reliability computed on ten randomly chosen blood draw sessions showed agreement ranging from 78.9% for facial expression to 95.0% for respiration. NIPS measurement methods are discussed in more detail in Williams et al (2009).

We restrict our attention here to studying the relationship of NIPS scores to HR and HRV data. Principal components analysis with varimax rotation (Mardia, Kent, & Bibby, 1979) was used to measure the relationships between the variables. For improved linearity and isolation of the response to pain, each variable was transformed to represent the change from baseline measurement.


In response to the painful stimulus, the LF and HF-powers decreased significantly during the heel lanced blood draw procedure. The −5.6 dB drop in LF-power (Table 2, Fig. 2) was sharper than the −1.8 dB drop in HF-power (Table 3, Fig. 3). The powers recovered to near-baseline levels in early and late recovery phases. The HRV response was mirrored in the HR data by an increase of 13.6 bpm after the painful stimulus (Table 4, Fig. 4). Subsequently, the HR dropped to 3.7 bpm above the baseline rate in early recovery and was maintained at 3.4 bpm above baseline in the late recovery phase. LF-power was 6.8 dB higher and HR was −5.8 bpm lower (Tables 2 and and4,4, Figs. 2 and and4)4) during the pain phase in venipunctured infants relative to the heel lanced infants.

Figure 2
Longitudinal model estimates of fixed-effects at 30 weeks post-menstrual age show changes in low-frequency (LF) power before, during, and after the painful procedure. Two curves are shown for ventilated and non-ventilated infants with heel lanced blood ...
Figure 3
Longitudinal model estimates of fixed-effects at 30 weeks post-menstrual age show changes in high-frequency (HF) power before, during, and after the painful procedure in ventilated and non-ventilated infants.
Figure 4
Longitudinal model estimates of fixed-effects at 30 weeks post-menstrual age show changes in heart rate before, during, and after the painful procedure. Two curves are shown for ventilated and non-ventilated infants with heel lanced blood draws. A third ...
Table 2
Parameters for longitudinal model of low-frequency HRV* and 95% confidence intervals (CI).
Table 3
Parameters for longitudinal model of high-frequency HRV* and 95% confidence intervals (CI).
Table 4
Parameters for longitudinal model of heart rate* and 95% confidence intervals (CI).

The response to pain differed sharply by mechanical ventilation status as noted by significant interactions between ventilation and phase in each longitudinal model (Tables 2, ,3,3, and and4).4). In mechanically ventilated infants, the LF, HF-powers and HR were almost unchanged from phase to phase (Figs. 2, ,3,3, and and4).4). The interaction terms between ventilation and the pain phase were of opposite sign as terms for the main effect of the pain, effectively canceling out the drops in powers and increase in heart rate observed in the spontaneously breathing infant. HF-power increased by 2.9 dB during the pain phase in ventilated newborns compared to a −1.8 dB decrease in non-ventilated newborns.

LF and HF-power showed growth with PMA of the infant (Fig. 5). The linear growth was modulated by a cubic term that reduced the growth rate at the edges of our observation region from 23 to 38 weeks PMA. In non-ventilated infants, the linear portion of the growth rate of HRV power was 1.0 dB/week in the LF band (Table 2) and 0.9 dB/week in the HF band (Table 3). The growth rate of HF-power was reduced substantially, by −0.5 dB/week, in infants receiving mechanical ventilation. At PMA 30 weeks, the expected power in the LF band was 28.1 dB, and in the HF band it was 22.4 dB. The expected HR at PMA 30 weeks was 165.5 bpm. HR had a quadratic dependence on PMA (Table 4), increasing from 149.8 bpm at 24 weeks to 167.2 bpm at 32.7 weeks, after which it declined with age, reaching 164.8 bpm at 36 weeks.

Figure 5
Longitudinal model estimates of fixed-effects show growth of baseline low-frequency (LF) power and high-frequency (HF) power differentiated by mechanical ventilation status.

Small changes in the pain response were observed as the infant matured. The significant interaction term between PMA and the late recovery phase in the models for LF-power (Table 2) and HR (Table 4) show that LF-power differential in late recovery increased at the rate of 0.5 dB/week and HR differential decreased by −1.0 bpm/week (Fig. 6).

Figure 6
Maturational changes in the response to pain in late recovery (5 to 7 minutes after last contact). The difference of low-frequency (LF) power between the late recovery phase and the baseline phase increased with neonatal age (top panel) while the difference ...

Random intercepts and PMA slopes reduced AIC in each model. Moreover, the standard deviation of the random intercepts in the LF and HF-power models exceeded the residual standard deviation (Tables 2 and and3).3). Taking account of the variation of levels and growth rates between individuals was therefore important in the modeling efforts. The HR model improved by addition of the interaction between ventilation and PMA as a random effect.

Longitudinal models based on data from 3-minute phase-windows yielded results that were similar to those obtained from 2-minute phase-windows described above. The type of blood draw and mechanical ventilation status continued to be statistically significant factors in the pain response. Among the differences were small but significant increases in the LF and HF powers in early recovery that were negated in the ventilated infants. The cubic modulation of growth rate was found to be unnecessary in the HF model. The heart rate model differed in the following respects: the heart rate was not significantly elevated over baseline in the early and late recovery phases; the reduction of heart rate in late recovery with increasing PMA was not significant, but increase in heart rate in the pain phase with increasing PMA acquired significance.

As another test, the data with 2-minute phase-windows were analyzed after retaining only 61 data sets that had not a single missing spectral power in any of the phases. The model for HF-power changed most, showing no significant pain response. However, the tendencies remained identical to the results obtained from the complete data, suggesting that there was insufficient statistical power in the model for the reduced set. Both LF and HF powers increased with age, but the growth rates were lower and cubic modulation was unnecessary. The LF power and heart rate continued to show a distinct response to pain of approximately the same magnitudes as described previously, and mechanical ventilation continued to be an important factor that differentiated the responses to pain. The pain-response results for the HF power presented herein may therefore be considered less robust than the responses of LF power and heart rate. There were too few venipunctured blood draws in this reduced set for a meaningful study of differences between responses to heel lance and to venipuncture.

It was observed from the principal components analysis (Table 5) that change of NIPS scores from the baseline were loaded strongly on the second principal component while HR and HRV variables were loaded more strongly on the first component. The first two principal components accounted for more than 75% of the variance.

Table 5
Principal component loadings on NIPS, HR, and HRV variables with summary of component variances. Each variable was represented as change from baseline for improved capture of pain response.


Reduction of HRV was observed in both bands for non-ventilated neonates during the heel lance procedure, accompanied by an increase in HR. The HF-power reduced significantly during the pain phase, although the −1.8 dB magnitude of the reduction in HF-power was smaller than the −5.6 dB reduction in LF-power. HR increased by 13.6 bpm in the 2 minutes following the pain stimulus. These observations are in qualitative agreement with prior studies (Lindh et al., 1997; Grunau et al., 2001) and in agreement with reports of reduced vagal tone when the infant is distressed by illness or pain (Arditi et al., 2006; Porges, 1992; Porter, Porges, & Marshall, 1988; van Ravenswaaij-Arts et al., 1995). It has also been suggested that infants might react to a painful stimulus with a fear paralysis reflex that is characterized by bradycardia and temporary inhibition of the sympathetic nervous system (Lagercrantz, Edner, Runold, & Yuan, 1995; Lindh et al., 1999).

It has been reported (Lindh et al., 1997) that HR increased and LF-power decreased in response to pain in preterm infants. Arditi et al (2006) found an increase in HR and decrease in vagal tone (HF-power) in response to a heel-prick procedure in healthy neonates. In a cross-sectional design, Grunau et al (2001) studied pain in 136 VLBW infants of 32 weeks PMA and found that the LF power and HF power decreased during blood collection followed by an increase in the recovery period. They used 200-s windows for data collection, and the recovery period was defined to begin immediately after the last contact. Although we found a tendency for increased powers in the early recovery period, the increases were not statistically significant, likely due to the later timing of our early recovery phase, from 2 to 4 minutes after last contact. Mean HR was reported by Grunau et al (2001) to increase from 160.4 bpm at baseline to 184.4 during squeeze, and easing to 171.4 bpm during recovery. An early study involving 35 preterm infants (McIntosh et al., 1993) reported an increase in HRV during a 5-minute period that included a heel-prick, which is in contrast to three later studies, including the present one. However, it may possibly be explained by the wide window that was used that included the recovery period in which HRV tends to increase (Grunau et al., 2001). Note that despite the use of common terminology such as LF and HF bands, the frequency limits of the bands varied between studies. Furthermore, differences in factors and units of measurement of band powers make quantitative comparisons of HRV between studies a difficult task. Despite the differences in study design, spectrum calculation methods, and models of data analysis, it is clear that our findings of HR and HRV response to pain are in qualitative agreement with previous studies by Lindh et al (1997) and Grunau et al (2001).

The attenuated response to pain from venipuncture that was found in LF-power and HR in this study parallels the behavioral measurements of reduced pain reported in other studies (Shah, Taddio, Bennett, & Speidel, 1997; Larsson, Tannfeldt, Lagercrantz, & Olsson, 1998; Eriksson, Gradin, & Schollin, 1999; Logan, 1999; Shah & Ohlsson, 2004; Ogawa, Ogihara, Fujiwara, Ito, Nakano, Nakayama, Hachiya, Fujimoto, Abe, Ban, Ikeda, & Tamai, 2005). It should be noted that our study was not designed to probe differences between pain levels in heel lance and venipuncture. Slightly less than 10% of blood draws were by venipuncture, therefore a more balanced design aimed toward studying HRV differences between the types of blood draw would be desirable to bolster our incidental finding.

Mechanical ventilation status, which has not received attention in previous studies, was found to be an important discriminating factor. The stress response to pain, as measured by HR and HRV, was substantially reduced in mechanically ventilated neonates. On the other hand, mechanical ventilation was associated with lower growth rate of HF-power. This result parallels a report that mechanically ventilated VLBW infants had lower growth of HF power than LF power with advancing age (Smith, Doig, & Dudley, 2004). In that study neither growth rate achieved statistical significance, while we observed statistically significant growth of LF power in mechanically ventilated infants. The LF and HF-powers increased in time, particularly for infants that were not mechanically ventilated, and the intersection of curves (Fig. 5) suggests that any advantage offered by mechanical ventilation to boost HRV may be prior to 33 weeks PMA. It is not clear whether mechanical ventilation served to reduce the stress response or whether the lowered stress response may be attributed to the sicker state of the ventilated infant. Usage of pain medications was low (Table 1) and it may be ruled out as a factor in the reduced pain response in ventilated infants. The low use of pain medications in our study population was due to a practice in our NICU that is guided by evidence of longer need for mechanical ventilation in infants who are on narcotics (Bhandari, Bergqvist, Kronsberg, Barton, & Anand; NEOPAIN Trial Investigators Group, 2005). Pain medications were only used for major procedures but not used routinely in ventilated infants. Non-pharmacological pain relief measures such as containment were equally used in both groups. Another possible interpretation of the findings is that mechanical ventilation affects HRV parameters strongly enough to mask stress response.

Growth of HRV with increasing PMA has been reported earlier (Doussard-Roosevelt, Porges, & McClenny, 1996; Sahni, Schulze, Kashyap, Ohira-Kist, Myers, & Fifer, 1999; Khattak, Padhye, Williams, Lasky, Moya, & Verklan, 2007) and is in agreement with our results. Our models also incorporated a significant cubic modulation term, suggesting that growth was reduced near both ends of our observation range from 23 to 38 weeks. One of the questions of interest was to measure changes in the response to pain as the neonate matures. In this regard, a tendency toward increased HR in older neonates in the pain phase was noted, however it was statistically significant only when the wider 3-minute phase-window was employed. Age-related changes were present in the late recovery phase, 5 to 7 minutes after the end of the painful procedure. The response of LF-power to the pain stimulus, as measured by the difference between LF-powers in the late recovery and baseline phases, increased at the rate of 0.5 dB/week. Similarly, the response of the HR to the pain stimulus, as measured by the difference between HR in the late recovery and baseline phases, decreased by −1.0 bpm/week. The increased LF-power and decreased HR in late recovery with increasing PMA suggests that sympathetic inputs are restored more quickly as the infant matures. Changes in HR and HRV response to pain over the course of weeks have not been reported in previous studies, to the best of our knowledge.

Baseline HR increased from 24 to 32.7 weeks PMA and decreased thereafter. This departure from the expectation of monotonically decreasing HR may be due to infants having learned to anticipate the painful blood draw in the minutes before the procedure began. Although generally asleep, a cue for anticipation may have been provided by the heel-warmer that was applied at the beginning of data collection, approximately 6–7 minutes before the beginning of the time-window of our baseline phase.

Lindh et al (1999) made a distinction between lancing and squeezing and found that LF power increased during lancing but decreased during squeezing. A limitation of our study is that it did not distinguish between lancing and squeezing. The two events occurred within seconds of each other and were treated as a single pain stimulus. Another limitation is that separation between effects of maturation and habituation on the pain response was not possible in our study design. The growth of HF and LF power with advancing age in the baseline phase lends support to the view that we observed maturational changes in the response to pain, however the impact of habituation cannot be ruled out. Measuring the maturation of pain responses of neonates continues to be an important field of investigation that may lead to improvements in the neonates’ experiences in the NICU. Venipuncture contributed to reduced pain over heel lance, as measured by the response of HRV and HR. It appears that mechanical ventilation is an important factor in the neonatal pain response that needs further study.

Behavioral information about pain response from the NIPS scores was loaded strongly on the second principal component while the physiological variables of HR and HRV were loaded more strongly on the first component. These components are orthogonal, which indicates that the behavioral information content was to a large degree independent of physiological information. This provides a basis to expect that multidimensional models involving behavioral and physiological variables (e.g. Gibbins, Stevens, McGrath, Yamada, Beyene, Breau, Camfield, Finley, Franck, Johnston, Howlett, McKeever, O’Brien, & Ohlsson, 2008) may provide a more complete picture of pain responses than either component can on its own.


This research was supported by NIHCD Grant R01 HD42639 from the National Institute for Child Health and Human Development to Dr. Lasky. Authors are grateful to Dr. M. Terese Verklan for her support during the data collection phase of this study.

Contributor Information

Nikhil S Padhye, The University of Texas School of Nursing at Houston, Center for Nursing Research, Houston, Texas.

Amber L Williams, The University of Texas Medical School at Houston, Center for Clinical Research and Evidence-Based Medicine, Houston, Texas.

Asif Z Khattak, The University of Texas Medical School at Houston, Department of Pediatrics, & Baylor University Medical Center at Dallas, Department of Neonatology, Dallas, Texas.

Robert E Lasky, The University of Texas Medical School at Houston, Center for Clinical Research and Evidence-Based Medicine, Houston, Texas.


  • Arditi H, Feldman R, Eidelman AI. Effects of human contact and vagal regulation on pain reactivity and visual attention in newborns. Developmental Psychobiology. 2006;48:561–573. [PubMed]
  • Bhandari V, Bergqvist LL, Kronsberg SS, Barton BA, Anand KJ. NEOPAIN Trial Investigators Group. Morphine administration and short-term pulmonary outcomes among ventilated preterm infants. Pediatrics. 2005;116:352–359. [PubMed]
  • Buskila D, Neumann L, Zmora E, Feldman M, Bolotin A, Press J. Pain sensitivity in prematurely born adolescents. Archives of Pediatrics and Adolescent Medicine. 2003;157:1079–1082. [PubMed]
  • Chang KL, Monahan KJ, Griffin MP, Lake DE, Moorman JR. Comparison and clinical application of frequency domain methods in analysis of neonatal heart rate time series. Annals of Biomedical Engineering. 2001;29:764–774. [PubMed]
  • Diggle PJ, Heagerty P, Liang K-Y, Zeger S. Analysis of longitudinal data. Oxford, UK: Oxford University Press; 2002.
  • Doussard-Roosevelt J, Porges SW, McClenny BD. Behavioural sleep states in very low birth weight preterm neonatal health and vagal maturation. Journal of Pediatric Psychology. 1996;21:785–802. [PubMed]
  • Eriksson M, Gradin M, Schollin J. Oral glucose and venepuncture reduce blood sampling pain in newborns. Early Human Development. 1999;55:211–218. [PubMed]
  • Fitzgerald M, Millard C, McIntosh N. Cutaneous hypersensitivity following peripheral tissue damage in newborn infants and its reversal with tropical anaesthesia. Pain. 1989;39:31–36. [PubMed]
  • Fox NA. Psychophysiological correlates of emotional reactivity during the first year of life. Developmental Psychobiology. 1989;25:364–372.
  • Gibbins S, Stevens B, McGrath PJ, Yamada J, Beyene J, Breau L, Camfield C, Finley A, Franck L, Johnston C, Howlett A, McKeever P, O’Brien K, Ohlsson A. Comparison of pain responses in infants of different gestational ages. Neonatology. 2008;93:10–18. [PubMed]
  • Goubet N, Clifton RK, Shah B. Learning about pain in preterm newborns. Journal of Developmental and Behavioral Pediatrics. 2001;22:418–424. [PubMed]
  • Grunau RV, Whitfield MF, Petrie JH. Pain sensitivity and temperament in extremely low-birth-weight premature toddlers and preterm and full-term controls. Pain. 1994;58:341–346. [PubMed]
  • Grunau RE, Oberlander TF, Whitfield MF, Fitzgerald C, Lee SK. Demographic and therapeutic determinants of pain reactivity in very low birth weight neonates at 32 weeks postconceptional age. Pediatrics. 2001;107:105–112. [PubMed]
  • Grunau RE. Early pain in preterm infants : a model of long-term effects. Clinics in Perinatology. 2002;29:373–394. [PubMed]
  • Grunau RE, Holsti L, Haley DW, Oberlander T, Weinberg J, Solimano A, Whitfield MF, Fitzgerald C, Yu W. Neonatal procedural pain exposure predicts lower cortisol and behavioral reactivity in preterm infants in the NICU. Pain. 2005;113:293–300. [PMC free article] [PubMed]
  • Huffman LC, Bryan Y, del Carmen R, Pederson F, Doussard-Roosevelt J, Porges S. Infant temperament and cardiac vagal tone: Assessments at twelve weeks of age. Child Development. 1998;69:624–635. [PubMed]
  • Johnston CC, Stevens BJ, Yang F, Horton L. Differential response to pain by very premature neonates. Pain. 1995;61:471–479. [PubMed]
  • Johnston CC, Stevens BJ. Experience in a neonatal intensive care unit affects pain response. Pediatrics. 1996;98:925–930. [PubMed]
  • Khattak AZ, Padhye NS, Williams AL, Lasky RE, Moya FR, Verklan MT. Longitudinal assessment of heart rate variability in very low birth weight infants during their NICU stay. Early Human Development. 2007;83:361–366. [PubMed]
  • Lagercrantz H, Edner A, Runold M, Yuan SZ. Development of autonomic control in infants. In: Rognum TO, editor. Sudden infant death syndrome. New trends in the nineties. Oslo: Scandinavium University Press; 1995. pp. 214–217.
  • Laguna P, Moody GB. Power spectral density of unevenly sampled data by least-square analysis: performance and application to heart rate signals. IEEE Transactions on Biomedical Engineering. 1998;45:698–715. [PubMed]
  • Larsson BA, Tannfeldt G, Lagercrantz H, Olsson GL. Venipuncture is more effective and less painful than heel lancing for blood tests in neonates. Pediatrics. 1998;101:882–886. [PubMed]
  • Lawrence J, Alcock D, McGrath P, Kay J, MacMurray SB, Dulberg C. The development of a tool to assess neonatal pain. Neonatal Network. 1993;12:59–66. [PubMed]
  • Lindh V, Wiklund U, Sandman PO, Hakansson S. Assessment of acute pain in preterm infants by evaluation of facial expression and frequency domain analysis of heart rate variability. Early Human Development. 1997;48:131–142. [PubMed]
  • Lindh V, Wiklund U, Hakansson S. Heel lancing in term new-born infants: an evaluation of pain by frequency domain analysis of heart rate variability. Pain. 1999;80:143–148. [PubMed]
  • Logan PW. Venepuncture versus heel prick for the collection of the newborn screening test. Australian Journal of Advanced Nursing. 1999;17:30–36. [PubMed]
  • Lomb NR. Least-squares frequency analysis of unequally spaced data. Astrophysics and Space Science. 1976;39:447–462.
  • Mardia KV, Kent JT, Bibby JM. Multivariate Analysis. London: Academic Press; 1979.
  • McIntosh N, Van Veen L, Brameyer H. The pain of heel prick and its measurement in preterm infants. Pain. 1993;52:71–74. [PubMed]
  • Oberlander TF, Grunau RE, Whitfield MF, Fitzgerald C, Pitfield S, Saul JP. Biobehavioral pain responses in formerly extremely low birth weight infants at four months corrected age. Pediatrics. 2000;105:e6. [PubMed]
  • Ogawa S, Ogihara T, Fujiwara E, Ito K, Nakano M, Nakayama S, Hachiya T, Fujimoto N, Abe H, Ban S, Ikeda E, Tamai H. Venepuncture is preferable to heel lance for blood sampling in term neonates. Archives of Disease in Childhood, Fetal and Neonatal Edition. 2005;90:F432–F436. [PMC free article] [PubMed]
  • Pan J, Tompkins WJ. A real-time QRS detection algorithm. IEEE Transactions on Biomedical Engineering. 1985;32:230–236. [PubMed]
  • Pinheiro JC, Bates DM. Mixed-effects models in S and S-PLUS. New York: Springer; 2000.
  • Porges SW. Vagal tone: A physiological marker of stress vulnerability. Pediatrics. 1992;90:498–504. [PubMed]
  • Porges SW, Doussard-Roosevelt JA, Portales AL, Suess PE. Cardiac vagal tone: Stability and relation to difficultness in infants and 3-year-Olds. Developmental Psychobiology. 1994;27:289–300. [PubMed]
  • Porter FL, Porges SW, Marshall RE. New-born pain cries and vagal tone – parallel changes in response to circumcision. Child Development. 1988;59:495–505. [PubMed]
  • Porter FL, Grunau RE, Anand KJS. Long-term effects of pain in infants. Journal of Developmental and Behavioral Pediatrics. 1999;20:253–261. [PubMed]
  • Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C. Cambridge: Cambridge University Press; 1992. pp. 575–584.
  • Pumprla J, Howorka K, Groves D, Chester M, Nolan J. Functional assessment of heart rate variability: Physiological basis and practical applications. International Journal of Cardiology. 2002;84:1–14. [PubMed]
  • Rangayyan RM. Biomedical signal analysis. A case-study approach. New York: John Wiley & Sons, Inc.; 2002.
  • Rosenstock EG, Cassuto Y, Zmora E. Heart rate variability in the neonate and infant: Analytical methods, physiological and clinical observations. Acta Paediatrica. 1999;88:477–482. [PubMed]
  • Sahni R, Schulze KF, Kashyap S, Ohira-Kist K, Myers MM, Fifer WP. Body position, sleep states and cardiorespiratory activity in developing low birth weight infants. Early Human Development. 1999;54:197–206. [PubMed]
  • Shah VS, Taddio A, Bennett S, Speidel BD. Neonatal pain response to heel stick vs. venepuncture for routine blood sampling. Archives of Disease in Childhood, Fetal and Neonatal Edition. 1997;77:F143–F144. [PMC free article] [PubMed]
  • Shah V, Ohlsson A. Venepuncture versus heel lance for blood sampling in term neonates. Cochrane Database Systematic Reviews. 2004 DOI: 10.1002/14651858.CD001452.pub2. [PubMed]
  • Shiavi R. Introduction to applied statistical signal analysis. San Diego: Academic Press; 1999.
  • Smith LS, Dmochowski PA, Muir DW, Kisilevsky BS. Estimated cardiac vagal tone predicts fetal responses to mother's and stranger's voices. Developmental Psychobiology. 2007;49:543–547. [PubMed]
  • Smith SL, Doig AK, Dudley WN. Characteristics of heart period variability in intubated very low birth weight infants with respiratory disease. Biology of the Neonate. 2004;86:269–274. [PubMed]
  • Stevens B, Johnston C, Petryshen P, Taddio A. Premature infant pain profile: development and initial validation. Clinical Journal of Pain. 1996;12:13–22. [PubMed]
  • Van Ravenswaaij-Arts CM, Hopman JC, Kollée LA, Stoelinga GB, Van Geijn HP. The influence of artificial ventilation on heart rate variability in very preterm infants. Pediatric Research. 1995;37:124–130. [PubMed]
  • Williams AL, Khattak AZ, Garza CN, Lasky RE. The behavioral pain response to heelstick in preterm neonates studied longitudinally: description, development, determinants, and components. Early Human Development. (article in press) [PubMed]