Autonomic responses to cold face stimulation in sickle cell disease: a time-varying model analysis

Sickle cell disease (SCD) is characterized by sudden onset of painful vaso-occlusive crises (VOC), which occur on top of the underlying chronic blood disorder. The mechanisms that trigger VOC remain elusive, but recent work suggests that autonomic dysfunction may be an important predisposing factor. Heart-rate variability has been employed in previous studies, but the derived indices have provided only limited univariate information about autonomic cardiovascular control in SCD. To circumvent this limitation, a time-varying modeling approach was applied to investigate the functional mechanisms relating blood pressure (BP) and respiration to heart rate and peripheral vascular resistance in healthy controls, untreated SCD subjects and SCD subjects undergoing chronic transfusion therapy. Measurements of respiration, heart rate, continuous noninvasive BP and peripheral vascular resistance were made before, during and after the application of cold face stimulation (CFS), which perturbs both the parasympathetic and sympathetic nervous systems. Cardiac baroreflex sensitivity estimated from the model was found to be impaired in nontransfused SCD subjects, but partially restored in SCD subjects undergoing transfusion therapy. Respiratory-cardiac coupling gain was decreased in SCD and remained unchanged by chronic transfusion. These results are consistent with autonomic dysfunction in the form of impaired parasympathetic control and sympathetic overactivity. As well, CFS led to a significant reduction in vascular resistance baroreflex sensitivity in the nontransfused SCD subjects but not in the other groups. This blunting of the baroreflex control of peripheral vascular resistance during elevated sympathetic drive could be a potential factor contributing to the triggering of VOC in SCD.


Introduction
In sickle cell disease (SCD), red blood cells with sickle hemoglobin (HbS) undergo shape transformation and become rigid when HbS polymerizes, following the release of oxygen to the surrounding tissues (Rees et al. 2010). The delay between oxygen release and HbS polymerization is small (Eaton and Hofrichter 1995), and thus, if transit time of these red blood cells through the microvasculature is prolonged, sickling can lead to blood flow obstruction in some capillary beds. The clinical manifestation of extensive obstruction of microvascular flow is vaso-occlusive crisis (VOC), which is characterized by episodes of pain and subsequently organ damage or even more severe consequences. However, the mechanisms by which transient regional vaso-occlusion cascades into full-blown VOC remain unclear.
Recent studies from our laboratory and others point to significant sympathovagal imbalance of the autonomic nervous system (ANS) in SCD. Relative dominance of sympathetic drive promotes peripheral vasoconstriction and, as such, prolongs the transit time of sickle cells as they traverse the microvasculature. Thus, it is quite possible that autonomic dysfunction may be a key factor that predisposes SCD patients to VOC (Alexy et al. 2010;Sangkatumvong et al. 2011;Connes and Coates 2013). Previous studies have reported that decreased parasympathetic activity in SCD and sickle cell trait (SCT) subjects, and elevated cardiac sympathetic activity in SCT subjects (Pearson et al. 2005;Connes et al. 2006;Martins Wde et al. 2012). The association of greater parasympathetic withdrawal with disease severity has also been reported (Pearson et al. 2005). Sudden death, which is known to be associated with SCD, has been attributed to autonomic dysfunction (Romero Mestre et al. 1997). We have recently shown marked parasympathetic withdrawal in response to transient hypoxia (Sangkatumvong et al. 2008) in SCD subjects. As well, sighs are more frequently associated with drops in peripheral microvascular perfusion in SCD patients, making pulmonary-vasoconstriction coupling a critical and measurable trigger of VOC (Sangkatumvong et al. 2011).
Spectral analysis of heart-rate variability (HRV) has been employed as a ubiquitous tool for assessing ANS function due to the ease and accuracy with which heart rate can be monitored in humans. All previous studies of ANS function in SCD employed noninvasive measures derived from HRV. However, there are important limitations with this methodology that are frequently overlooked. Some prominent researchers have called into question the assumptions that underlie the use of spectral analysis for delineating the sympathetic from parasympathetic components of cardiac autonomic drive (Eckberg 1997;Malpas 2002). As well, changes or differences in breathing frequency, tidal volume, or ventilatory pattern can significantly confound the interpretation of autonomic activity that one derives from HRV spectral analysis (Brown et al. 1993). To circumvent the aforementioned limitations, in this study, we have supplemented conventional HRV spectral analysis with a modelbased approach that allows us to quantify the strengths of the key functional mechanisms (such as baroreflex control of heart rate and respiratory-cardiac coupling) that contribute to HRV. This approach has also been applied to evaluate the relative contributions of the major functional mechanisms that account for fluctuations in peripheral resistance.
Since it is known that VOC can be triggered by cold (Resar and Oski 1991), we hypothesized that perturbing the ANS with cold face stimulation (CFS) would provide greater delineation between the responses in SCD subjects and those elicited in healthy controls. Cold face stimulation is a noninvasive test that activates both the cardiac parasympathetic and the peripheral sympathetic nervous system (Khurana et al. 1980). In normals, CFS generally leads to an initial bradycardia due to vagal stimulation, followed by a rise in BP due to sympathetically induced peripheral vasoconstriction. In this study, CFS was applied to healthy controls, SCD subjects participating in a chronic transfusion program, as well as SCD subjects not undergoing chronic transfusion. Subjects on chronic transfusion receive sufficient blood transfusion every 3 weeks to suppress hemoglobin S concentration below 30% and maintain a hemoglobin level above 9.5 g/dL. We hypothesized that the SCD subjects not participating in the transfusion program would exhibit altered autonomic function, but those in the chronic transfusion program would show partial recovery toward normal ANS status. To test these hypotheses, we incorporated the ability to accommodate time-varying changes in the parameters into the model, in order to take into account potential changes in the underlying dynamics resulting from transient exposure to the cold face stimulus.

Participants
All experiments were conducted at Children's Hospital Los Angeles (CHLA). The study protocol was approved by the Committee on Clinical Investigations (institutional review board of CHLA). Participants were selected from healthy controls and SCD patients who were African American and 10 years or older. Some of the SCD subjects were participants in the chronic transfusion program. Those who recently had crisis or were hospitalized for sickle-related event (less than 2 weeks prior to the study), or unable to lay flat for 30 min or move from laying to standing posture were excluded. The participants were divided into three treatment groups: healthy controls (CTL), nontransfused SCD (NTF), and chronically transfused SCD (CTF). All but two SCD subjects were SS genotype (the two subjects were SB0 thalassemia and included in NTF group). There were (male/female) 4/5, 4/7 and 5/3 subjects in each respective group. The grouped average ages (mean AE standard error) were 25.7 AE 3.5, 23.1 AE 3.2 and 19.9 AE 3.7 years, respectively. Hemoglobin count in each group (CTL, NTF, and CTF) was 13.6 AE 0.91, 9.0 AE 0.14 and 9.8 AE 0.19 g/dL. Percentage of sickle hemoglobin were 0, 72 AE 1.6 and 32 AE 1.9, respectively. Oxygen saturation level (SpO 2 ) in each group was 97.9 AE 0.3, 96.6 AE 0.7 and 96.8 AE 1.1, respectively. For CTF subjects, measurements were obtained before their routine blood transfusion.
Written informed consent/assent was obtained before participation in the study. All subjects completed a physical examination and history to evaluate for wellness prior to the exam. A data safety monitoring plan included a pre and post-test questionnaire to evaluate for the presence of pain and sickle cell related vaso-occlusive symptoms. The post-test questionnaire was performed immediately following the study and 24-72 h later. All subjects tolerated CFS and none of the subjects reported symptoms at either time point.

Protocols
The experimental protocol was designed to measure the physiological responses to the CFS. The subject rested comfortably in recumbent position with their torso elevated at a 30-degree angle. The subject was allowed to rest quietly for 10 min to allow the physiological measurements to stabilize. An ice pack was then placed on the subject's forehead for 60 sec. After the removal of the ice pack, the subject was again allowed to rest quietly for 10 min.
Respiration and electrocardiogram (ECG) were recorded continuously using Lifeshirt physiologic monitoring system (VivoMetrics, Ventura, CA). Continuous blood pressure (BP) was monitored noninvasively using Nexfin (BMEYE, Ansterdam, The Netherlands). Peripheral arterial tonometry (WatchPAT-100, Itamar Medical, Caesarea, Israel) was employed for assessment of pulse wave amplitude at the level of the finger. "PAT" noninvasively obtains a beat-to-beat plethysmographic recording of the finger arterial pulse wave amplitude. During vasoconstriction when peripheral vascular resistance increases, the pulse amplitude of PAT (PATamp) decreases (Schnall et al. 1999;Grote et al. 2003). Therefore, we used the PAT signal to detect peripheral vasoconstriction, and assumed PATamp to be inversely related to peripheral vascular resistance. Respiration, ECG and continuous BP were recorded at 200 Hz; PAT was recorded at 40 Hz.

Data preprocessing
The beat-to-beat interval of the R-to-R wave on the ECG (RRI) and tidal volume (V T ) were derived from the ECG and respiration recordings using VivoLogic software (VivoMetrics, Ventura, CA). Beat-to-beat systolic and diastolic blood pressure (SBP and DBP) were extracted from the peak and nadir of each beat on the continuous BP recording. Beat-to-beat mean arterial pressure (MAP) was approximated as a sum of 1/3 of SBP and 2/3 of DBP (Levy et al. 2007). Beat-to-beat PATamp was defined as the difference between the peak and nadir of each beat on the PAT recording. Because the PAT signal is given in arbitrary units, the beat-to-beat values of PATamp were normalized by the mean and standard deviation (STD) of the PATamp time-series during baseline (1-5 min of quiet data recording before the CFS) to allow comparison across subjects. Normalized PATamp (PATampN) was computed as follows: The above equation yields a measure of peripheral vasoconstriction/vasodilation relative to the baseline PAT amplitude in each subject, similar to a z-score.
All beat-to-beat signals were interpolated and then uniformly sampled at 2 Hz. The V T signal was also down-sampled to 2 Hz. R-R interval, V T , SBP, MAP, and PATampN were aligned such that time 0 indicates the onset of the CFS. All signals were subjected to 5th-order polynomial trend removal such that each signal reflected fluctuations around zero. Further, a lowpass filter using Kaiser window (Oppenheim et al. 1999) with a passband of 0-0.5 Hz and a stopband at 0.85-1 Hz was applied to remove high-frequency (HF) noise.

Modeling and parameter estimation
To determine the effect of fluctuations in V T and BP on changes in RRI and peripheral resistance, reflected by PA-TampN, in response to CFS, we employed two linear time-varying models whose structures were similar to that published by our group previously (Chaicharn et al. 2009). For both models, any effects of autonomic activity on RRI and PATampN were assumed to have been reflected in V T and BP fluctuations. The first model was the model of HRV. Fluctuations in RRI (DRRI) were assumed to be governed by two functional autonomic mechanisms: the arterial baroreflex (ABR) and the respiratory-cardiac coupling (RCC). Arterial baroreflex relates fluctuations in SBP (DSBP) to DRRI and RCC relates fluctuations in V T (DV T ) to DRRI. e RRI represents extraneous influences that are not captured by the ABR and RCC. The mathematical representation of the HRV model is as follows: A second model was proposed for characterizing the key factors influencing peripheral vascular resistance variability. As peripheral resistance is modulated by the baroreflexes, we assumed a dynamic relationship between fluctuations in PATampN and fluctuations in BP. It is also well established from peroneal nerve activity measurements that respiration modulates sympathetic neural outflow (Eckberg 1995). In turn, the respiratory-related oscillations in sympathetic activity exert a strong modulatory influence on peripheral vascular tone (Malpas 2002).
Previous studies have also shown that deep breaths/sighs trigger a peripheral vasoconstriction response (Bolton et al. 1936;Sangkatumvong et al. 2011). Thus, fluctuations in PATampN (DPATampN) were assumed to be produced through two functional mechanisms (Chalacheva and Khoo 2013): (1) baroreflex control of peripheral vascular conductance (BPC), which relates fluctuations in MAP (DMAP) to DPATampN; and (2) respiratoryperipheral vascular conductance coupling (RPC), which relates fluctuations in DV T to DPATampN. e PATampN represents extraneous influences in DPATampN that are not captured by the BPC and RPC. The mathematical representation of the peripheral resistance variability model is as follows: In equation (2) and (3), h(t,i) represents the timevarying impulse response of each mechanism. An impulse response quantifies the dynamics of the changes in the output resulting from a brief and abrupt unit increase in the input. For example, h ABR quantifies the dynamics of the DRRI resulting from an abrupt increase in SBP of 1 mmHg. The impulse responses in this study were allowed to vary with time, thus time-varying, in order to capture changes in the system in different phases of the CFS. M represents the memory of the system. Based on the nature of the dynamics of interest, each impulse response was assumed to persist for 50 sampling intervals where each sampling interval was 0.5 sec (i.e., 25 sec).
It should be noted that equations (2) and (3) provide open-loop representations of the effects of BP and respiration on DRRI and DPATampN, respectively. DRRI and changes in peripheral resistance, in turn, are expected to subsequently affect BP, as the measurements were made with these systems operating in closed-loop mode. In order to delineate the feedback from feedforward dynamics of these closed-loop systems, it was necessary to explicitly impose causality on the models by introducing delays, s, into each of the proposed mechanism. Based on the physiological considerations, s ABR and s RCC were assumed to range from 0.5 to 3 sec (Belozeroff et al. 2002(Belozeroff et al. , 2003 and from À3 to 0 sec (Mullen et al. 1997;Belozeroff et al. 2002), respectively. s BPC and s RPC were assumed to fall within 3.5-8 and 0-3 sec, respectively. The final estimates for each of these delays in each dataset analyzed were obtained by selecting the candidate solution that yielded the lowest criterion function value (see paragraph below).
To estimate the impulse responses, we employed a kernel expansion technique where each impulse response was represented as a weighted sum of a set of basic functions. This approach was chosen because it greatly reduces the number of parameters to be estimated and thus improves the estimation accuracy even when applied on a relatively short and noise-contaminated data segment (Marmarelis 1993). Meixner basis functions (MBF) were employed as the basis functions for the impulse responses because of their exponential decaying characteristics. In addition, they have a parameter, called the order of generalization, which allows us to have control over the rise time of the functions, making the approximation of impulse responses with slow start more accurate (Asyali and Juusola 2005). In this study, the allowed ranges of order of generalization for each mechanism were ABR: 1-5, RCC: 0-5, BPC: 2-8, and RPC: 2-8. The allowed ranges of Meixner function order (the total number of MBFs used in the expansion of an impulse response) for each mechanism were ABR: 3-6, RCC: 3-6, BPC: 2-6, and RPC: 2-6. The unknown expansion coefficients, that is the weights of the MBFs, were estimated using least-squares minimization. This least-squares minimization process was repeated for each combination of the delays, order of generalization, and Meixner function order. Minimum description length (MDL) was employed as a measure of the balance between quality of data fitting and complexity of the model (Rissanen 1982). The optimal model was selected based on the global minimum MDL.
To estimate the time-varying coefficients of the impulse responses, a recursive least-squares (RLS) algorithm was employed (Haykin 1996). Recursive least squares minimizes the weighted squared error between the model prediction and the data. In this case, the forgetting factor method was chosen, in which the weighting coefficients take the form of an exponential function where the decay rate is controlled by a parameter called the forgetting factor (k). k reflects the memory of the adaptive filter and can vary from 0 to 1. When k = 1, all data before the present time are used in the estimation of the current parameter. Smaller values of k imply that the most recent data points are weighted more heavily; however, the tradeoff is that there is greater variability in the parameter estimates as there is effectively less averaging of data. In this study, we allowed k to vary between 0.85 and 0.97. The value of k that minimized the squared error between the model prediction and the data was selected.
To reduce the variability in parameter estimates, we employed the following technique. Upon obtaining the initial estimate of the impulse responses using the method described above, the model prediction and the residuals (difference between data and model prediction) were calculated. Surrogate residuals were generated from the original residuals time-series by applying the Amplitude-Adjusted Fourier Transform (AAFT) technique (Theiler et al. 1992;Garde et al. 2001). This AAFT technique allowed us to reshuffle the original residuals time-series such that its Fourier Transform magnitude spectrum was preserved but the phase component was shuffled randomly. The surrogate residuals were subsequently added back to the original predicted output to generate a new output time-series. The estimation process was repeated by treating the new output time-series as if it were a new measurement of the output. The parameter estimation described above was applied to the "new" output timeseries. This process was repeated 50 timeseach time, the impulse response of the model was estimated. The median impulse response, computed from the 50 estimated responses, was taken to represent the final estimated impulse response.
Once the impulse response was obtained, its transfer function was computed by applying fast Fourier transform to the impulse response at each time step. Compact descriptors were then extracted from the transfer function at each time step. These descriptors were: (1) the low-frequency (LF) gain, the average transfer function magnitude between 0.04 and 0.15 Hz; (2) the HF gain, the average transfer function magnitude between 0.15 and 0.4 Hz; and (3) the overall dynamic gain, the average transfer function magnitude between 0.04 and 0.4 Hz. Based on results obtained from our previous work (Belozeroff et al. 2002(Belozeroff et al. , 2003Chaicharn et al. 2009) in other applications, we focused attention on the following estimated model parameters: (a) Baroreflex control of heart rate: the average gain of the ABR transfer function in the LF region (|H ABR | LF ). Decreased |H ABR | LF would be expected to reflect reduced vagal tone and/or increased sympathetic tone. (b) Respiratory-cardiac coupling: the average gains of the RCC transfer function in the LF (|H RCC | LF ), HF (|H RCC | HF ) and overall frequency range (|H RCC | Overall ) were taken to reflect vagal modulation of heart rate. (c) Baroreflex control of peripheral resistance: the average gain of the BPC transfer function in the LF region (|H BPC | LF ) was taken to represent sympathetic modu-lation of peripheral vascular conductance. However, as PATampN does not provide an absolute measure of peripheral vascular conductance, we were able only to determine the relative changes in |H BPC | LF from baseline during and after CFS in each subject group. (d) Respiratory-peripheral vascular conductance coupling: the average gain of the RPC transfer function in the HF region (|H RPC | HF ) was taken to reflect the strength of the respiratory-vasoconstriction reflex. As in (c), since PATampN provides only a relative measure of peripheral vascular conductance, we were able only to determine how CFS may have elicited relative changes in |H RPC | HF from baseline in each subject group.

Other calculations
For each dataset, the median values of physiological measurements, namely RRI, SBP, DBP, and PATampN, were computed during baseline (1-4 min before the cold face onset) and every half-minute interval from the onset of the CFS (0 min) for 2 min. The time-varying power spectral density of RRI, SBP, and PATampN were also computed using the RLS algorithm with parameter error reduction technique described in the previous section. The spectral indices: low-and high-frequency power (denoted LFP and HFP, respectively) were extracted from the 0.04-0.15 Hz and 0.15-0.4 Hz frequency bands of the power spectral densities of RRI and SBP. For the RRI, the LF/HF ratio (LHR RRI ) was also computed as another spectral index to reflect the "sympathovagal balance" (Task Force, 1996). Then the mean values of these spectral indices were computed during baseline and every half-minute interval from the onset of the CFS for 2 min.

Statistical tests
Mixed model repeated measures analysis of variance (mixed model ANOVA) was applied to each of the physiological measurements, spectral indices and modelderived descriptors -LF, HF and overall dynamic gains. The model-derived descriptors were log transformed to satisfy the assumptions of ANOVA (normality and equal variance). The mixed model ANOVA tested for significant main effects for treatment group (CTL, NTF, and CTF) and time (baseline, 0-0.5, 0.5-1, 1-1.5 and 1.5-2 min), and significant interaction effects. If significant differences were detected in the main effects and/or interactions, post hoc pairwise comparisons were performed using the Tukey's honestly significant difference test. All statistical analyses were performed using JMP statistical software, version 11.0 (SAS Institute Inc., Cary, NC).

Physiological measurements
Changes in physiological measurements, namely heart rate, BP, and peripheral resistance (represented by PA-TampN), in response to the CFS are summarized in Table 1. Representative data of RRI, SBP, and PATampN are shown in Figure 1. Only CTL subjects showed apparent increased RRI during the CFS (0-1 min). However, this increase in RRI was not enough to differentiate CTL subjects from SCD subjects. For all treatment groups, both SBP and DBP were maintained during the first half of the stimulation then increased during the second half (Table 1). Post hoc comparisons showed that SBP during the second half of the CFS was overall significantly higher than baseline and 1.5-2 min after the cold face onset. Likewise, DBP during the second half of the stimulation was overall significantly higher than all time-bins being tested. The CFS caused PATampN to decrease, reflecting  vasoconstriction, during and after the application of ice pack as indicated in significant time factor (P < 0.001). Post hoc comparisons on the time effect showed that the decreased PATampN occurred from the onset of the stimulation to 1.5 min after relative to the baseline.

Spectral indices of cardiovascular variability
Changes in spectral indices in response to the CFS are detailed in Table 2. LFP RRI , reflecting a combination of cardiac sympathetic and parasympathetic activity, was   significantly different across subject groups (P = 0.015), being substantially lower in NTF compared to CTL. The same tendency was displayed by HFP RRI (generally taken to reflect cardiac parasympathetic activity) but this did not achieve statistical significance. As well, both LFP RRI and HFP RRI in the CTF subjects tended to fall between their corresponding values in the CTL and NTF groups at each of the time-points considered. However, LHR RRI (reflecting "sympathovagal balance") displayed no significant differences among treatment groups. LFP SBP , which reflects BP oscillations that could be mediated sympathetically through the baroreflexes, displayed a tendency to increase during and after CFS in the CTL subjects, but this change did not attain statistical significance. There were also no significant differences in the responses among treatment groups.

Estimated model parameters
The responses to CFS of the heart-rate baroreflex gain (| H ABR | LF ) derived from the HRV model are shown in Figure 2. There was a tendency for |H ABR | LF during baseline to be highest in CTL, followed by CTF and lastly NTF subjects, but these differences were not statistically significant. However, application of the CFS accentuated these group differences, in that CFS led to an increase in |H ABR | LF in CTL subjects but not in the NTF or the CTF subjects. In the CTF subjects, there was a small tendency for |H ABR | LF to increase only after a delay, following removal of the CFS. The mixed model ANOVA showed a significant main effect for treatment group (P = 0.003). Pairwise comparisons revealed significant differences between CTL and NTF at various timepoints during and after CFS application, whereas there was no significant difference between CTL and CTF. Figure 3 displays the effect of CFS on RCC gain, |H RCC | LF . |H RCC | LF in CTL was higher overall compared to both SCD groups during baseline, CFS and following CFS application (P = 0.009). Thus, unlike the heart-rate baroreflex gain, RCC gain in CTL did not change during CFS, although it displayed a tendency (not significant) to be slightly decreased following removal of the cold face stimulus. |H RCC | LF did not change with CFS in both NTF and CTF subjects.
The response to CFS of the gain representing baroreflex modulation of peripheral vascular resistance (|H BPC | LF ) is shown in Figure 4. |H BPC | LF in all treatment groups significantly decreased relative to baseline from the onset of CFS until the half minute following CFS (P < 0.001). However, pairwise comparison between baseline and the last half-minute of CFS confirmed that the NTF subjects displayed the largest drop in |H BPC | LF during CFS, as illustrated in Figure 4. Figure 5 shows that the gain representing the respiratory-peripheral vasoconstriction reflex, |H RPC | HF , was transiently depressed by CFS in all the treatment groups (P < 0.001). This was the case for the LF and overall (0.04-0.4 Hz) frequency ranges as well. Post hoc comparisons showed that |H RPC | HF during the second half of the CFS and half a minute after the stimulation were significantly lower than the baseline. The trend of changes over time in all treatment groups appeared to be similar: pairwise comparisons did not reveal any significant differences in |H RPC | HF among treatment groups or interactions between time and treatment groups.

Cold face stimulation in autonomic testing
Cold face stimulation has been employed in previous studies as a noninvasive test of autonomic function since it activates both the cardiac parasympathetic and the peripheral sympathetic nervous systems. The characteristic responses induced by the combined activation include bradycardia and peripheral vasoconstriction with subsequent increase in BP (Khurana et al. 1980;Hilz et al. 1999). Bradycardia is mediated by efferent cardiac parasympathetic pathway while peripheral vasoconstriction is mediated by efferent sympathetic pathway (Heistad et al. 1968;Hilz and Dutsch 2006). Cold face stimulation is unique in that it elicits increases in both parasympathetic and sympathetic output, whereas most other common autonomic tests, such as tilt, produce changes in parasympathetic and sympathetic activity that are opposite in direction. It was reported that subjects with a disturbance in the integrity of the trigeminalbrainstem-vagal reflex arc had absent or depressed bradycardia response to CFS (Khurana et al. 1980). Cold face stimulation is a modification of the diving reflex, which occurs with immersion of the face in water. Previous studies demonstrated that CFS could reproduce cardiac parasympathetic response with similar sensitivity to the diving reflex (Khurana et al. 1980;Hilz et al. 1999;Hilz and Dutsch 2006) and thus could be used to evaluate parasympathetic function in patients with limited ability to cooperate (Hilz et al. 1999). As such, CFS has been employed as a means to perturb the ANS in various pathological conditions such as in diabetes mellitus, brainstem stroke, Shy-Drager syndrome, familial dysautonomia, obstructive sleep apnea, or SCT/anemia (Khurana et al. 1980;Hilz et al. 1999;Chaicharn et al. 2009;Martins Wde et al. 2012).

Model-based analysis
Using a model-based approach, we were able to detect differences in autonomic function among CTL, NTF, and CTF subjects using only noninvasive measurements of respiration, heart rate, and peripheral vascular perfusion. Cold face stimulation served the purpose of accentuating the autonomic differences across these subject groups. The model complemented the conventional spectral analyses of HRV and BP that were also applied, enabling us to detect the autonomic differences between the CTL and SCD subjects, as well as the autonomic effects of transfusion therapy on SCD subjects (Figs. 2, 3) that were not detectable in the raw data (Table 1). Since the model allowed for a decomposition of the complex closed-loop mechanisms regulating heart rate and BP into simpler functional input-output components, we were able to quantify the strength of each investigated mechanism while adjusting for other contributing factors. For instance, CFS is also known to affect breathing pattern and respiratory rate, and these changes can act to con-found the autonomic effects of CFS. In our model, we incorporated respiration as one of the key inputs that can affect heart rate and peripheral vascular resistance. As well, we allowed for the likelihood that application and subsequent removal of the cold stimulus would affect the parameters (e.g., gains) characterizing the model-thus, the model equations were formulated to be time-varying in the parameters. Adaptive filtering methodology was employed to track the changes in the model parameters during and after CFS (Haykin 1996). We also introduced a new technique of reducing the error in the estimates of the time-varying parameters while enabling the RLS estimation algorithm to track relatively fast changes in the system.

Autonomic dysfunction in SCD
The present work is the first study to employ a rigorous quantitative framework for analyzing the dynamic autonomic responses of SCD subjects to CFS, and to compare the responses of those patients who have not undergone transfusion therapy (NTF) against those subjects who are undergoing transfusion therapy (CTF). Heart-rate variability, in particular LFP RRI , was substantially lower in NTF subjects compared to the CTL and CTF groups ( Table 2). These results, together with the seemingly lower HFP RRI in NTF group (although not significant), imply that CTL and CTF had higher parasympathetic activity compared to NTF, suggesting that transfusion may have partially corrected the ANS defect seen in SCD. Our inability to detect significant differences in HFP RRI across treatment groups was likely due to the confounding effect of respiration. Heart-rate variability is known to be most reliable when applied under conditions of controlled breathing, whereas the CFS tended to affect the subject's breathing pattern. To minimize the confounding effect of respiration, we used the model to estimate the RCC gain, |H RCC | LF (Belozeroff et al. 2002(Belozeroff et al. , 2003Chaicharn et al. 2009). Although in all subject groups respiratory-coupling gain did not change during CFS, |H RCC | LF was significantly lower in both SCD groups compared with CTL (Fig. 3). These findings underscore the central point that SCD subjects exhibit impaired parasympathetic modulation of heart rate. This is consistent with previous studies that found lowered parasympathetic activity in SCT subjects (Connes et al. 2006;Hedreville et al. 2008;Sangkatumvong et al. 2011 to healthy volunteers, suggesting lower parasympathetic reserves in the former groups. The gain of the baroreflex control of heart rate was highest in CTL, followed by CTF, and lastly NTF. During CFS, baroreflex gain increased in CTL, but remained unchanged in NTF. The increase in heart-rate baroreflex gain with CFS in normals has been documented in previous studies (Hilz et al. 1999;Chaicharn et al. 2009). In CTF subjects, there was some tendency for baroreflex gain to increase but only after a delay. Taken together with the results for |H RCC | LF , these observations suggest impaired parasympathetic modulation of heart rate in NTF and CTF relative to CTL, but partial restoration of baroreflex control toward normality in the CTF group. However, the physiological mechanisms for the improved baroreflex gain following chronic transfusion therapy remain unknown at this time. Our findings are consistent with a previous study (Martins Wde et al. 2012), that reported reduction in baroreflex sensitivity in both SCD and SCT subjects, as assessed by the more invasive technique of administration of nitroglycerine and phenylephrine. Martins et al. also found that subjects with iron deficiency anemia had normal baroreflex sensitivity, suggesting that the reduction in baroreflex sensitivity in SCD was likely not the result of chronic anemia alone. On the other hand, another study reported normal baroreflex sensitivity in SCD subjects during supine and standing (Kim et al. 2009).
Although "sympathovagal balance", as characterized by LHR RRI , was not different across the subject groups, the substantially lower heart-rate baroreflex gain in the SCD subjects relative to CTL suggests that the SCD subjects had elevated sympathetic tone. An inverse association between baroreflex sensitivity and sympathetic tone in a variety of conditions (e.g., hypertension, heart failure, postmyocardial infarction) has been well documented in past studies (Grassi et al. 1995;La Rovere et al. 2008). To account for this association, it has been argued that impaired baroreflex sensitivity leads to reduction in reflex sympathetic restraint, which in turns allows sympathetic overactivity to occur. In the CTL subjects of this study, application of CFS led to an increase in vagal tone (hence the bradycardia), which also increased the gain of baroreflex control of heart rate. This CFS-induced increase in heart-rate baroreflex gain did not occur in the NTF subjects; but the CTF subjects displayed a delayed, small increase.
Peripheral vasoconstriction in response to CFS occurred in all three subject groups, as evidenced by statistically significant decreases in PATampN and corresponding increases in SBP (Table 1). There were no differences in the time-course of these changes in PA-TampN and SBP across subject groups. Reduction in PATampN was expected as CFS has been shown to induce peripheral vasoconstriction through the efferent sympathetic pathway (Khurana et al. 1980;Hilz et al. 1999). Thus, the fact that both SCD groups of subjects exhibited similar time-courses of PATampN reduction indicates that sympathetic modulation of peripheral vascular resistance remains intact in SCD. However, an important limitation of our study is that peripheral arterial tonometry does not provide an absolute measure of peripheral vascular resistance. Based on our measurements of PATampN, it was possible to only compare changes (elicited by some stimulus) relative to baseline levels. Thus, we were not able to make direct inferences about sympathetic tone in the different subject groups.
The gain corresponding to baroreflex control of peripheral vascular conductance decreased transiently following CFS in all subject groups, again indicating intact sympathetic modulation of peripheral vasoconstriction in spite of the presence of SCD (Fig. 4). The time-courses for the CFS-induced reduction in baroreflex gain were largely similar between CTL and CTF subjects; however, the NTF subjects displayed reductions in |H BPC | LF from baseline that were distinctively larger. Although the reason for this difference remains unclear at this point, we speculate that there may be an ominous clinical implication to this finding. A large reduction in the vascular resistance baroreflex gain during cold face exposure, coincident with rising BP means that there may be insufficient negative feedback from the baroreflexes to effectively counter the progressive peripheral vasoconstriction. The reduction in peripheral flow due to this transient blunting of the vascular resistance baroreflex in vulnerable SCD subjects could greatly enhance the chances for triggering local VOC and ultimately a VOC cascade.
A number of limitations in this study should be addressed. First, male and female subjects were combined for analysis of their autonomic responses to CFS because the sample size was not large enough to include gender as an independent factor although autonomic function can differ between genders. A separate analysis (not reported here) was performed to investigate the effect of gender difference. We found that there was no statistical difference in response to CFS between genders. Another limitation was the age of the recruited subjects (youngest subject was 10 years old). While age could be a strong factor affecting autonomic function, we found that exclusion of these young subjects did not change our findings in this study. Lastly, two SB0 thalassemia subjects were included in the NTF group while the rest of the SCD subjects were SS genotype. We found that the findings remained unaltered after excluding these two subjects. Thus, although differences in SCD genotype could contribute to differences in autonomic function, future studies involving larger numbers of subjects with various genotypes would be required to validate this possibility.
In summary, this study investigated the effect of CFS on SCD subjects who had and had not gone through transfusion therapy program using a time-varying modeling approach as a complement to conventional spectral analyses of heart rate and BP variability. Using the modeling approach, we found that regardless of their participation in the transfusion program, SCD subjects had reduced respiratory sinus arrhythmia, suggesting that their low parasympathetic activity is not corrected by transfusion. On the other hand, only nontransfused SCD subjects had impaired cardiac baroreflex compared to controls and chronically transfused SCD subjects, suggesting improvement in the baroreflex sensitivity following transfusion therapy. In addition, we found that the lungperipheral vascular conductance coupling appeared to be similar in all treatment groups. However, the sensitivity of the baroreflex control of vascular resistance in nontransfused SCD subjects was distinctly lower during CFS compared to its baseline. This blunting of the baroreflex during elevated sympathetic drive (such as induced by CFS) could be a potential factor that contributes to the triggering of vaso-occlusion in SCD.