# Detection of Impaired Sympathetic Cerebrovascular Control Using Functional Biomarkers Based on Principal Dynamic Mode Analysis

^{1}Department of Electrical Engineering, COMSATS Institute of Information Technology, Sahiwal, Pakistan^{2}Wellington Medical Technology Group, Centre for Translational Physiology, University of Otago, Wellington, New Zealand^{3}School of Engineering and Computer Science, Victoria University of Wellington, Wellington, New Zealand

This study sought to determine whether models of cerebrovascular function based on Laguerre-Volterra kernels that account for nonlinear cerebral blood flow (CBF) dynamics can detect the effects of functional cerebral sympathetic blockade. We retrospectively analyzed continuous beat-to-beat blood pressure, middle cerebral blood velocity, and partial-pressure of end-tidal CO_{2} (P_{ET}CO_{2}) recordings from eighteen healthy individuals who were treated with either an oral dose of the α_{1}-adrenergic receptor blocker Prazosin or a placebo treatment. The global principal dynamic modes (PDMs) were analyzed using Laguerre-Volterra kernels to examine the nonlinear system dynamics. Our principal findings were: (1) very low frequency (<0.03 Hz) linear components of first-order kernels for BP and P_{ET}CO_{2} are mutually coupled to CBF dynamics with the ability to separate individuals between control and blockade conditions, and (2) the gains of the nonlinear functions associated with low-pass and ≈0.03 Hz global PDMs for the BP are sensitive to sympathetic blockade. Collectively these results suggest that very low frequency global PDMs for BP may have potential utility as functional biomarkers of sympathetic neurovascular dysfunction which can occur in conditions like autonomic failure, stroke and traumatic brain injury.

## 1. Introduction

Biomarkers of sympathetic cerebrovascular dysfunction have proved difficult to identify owing to the linear and stationary assumptions of conventional haemodynamic modeling techniques (Zhang et al., 1998; Hamner et al., 2010). To address this issue, different approaches have recently been proposed that can either identify or explicitly account for nonlinear (Mitsis et al., 2009) and nonstationary (Latka et al., 2005) characteristics of cerebral haemodynamics (Peng et al., 2008; Saleem et al., 2016). An important class of these methods include the novel approach of principal dynamic modes (PDMs), in connection with the kernels of Laguerre-Volterra models, that was introduced by Marmarelis (1997, 2004) in order to quantify the nonlinear dynamics of relationships between blood pressure (BP), cerebral blood flow (CBF) and other covariate inputs e.g., partial-pressure of end-tidal CO_{2} (P_{ET}CO_{2}) (Mitsis et al., 2002, 2004, 2005).

Several recent studies have utilized these models (i.e., PDMs and Laguerre-Volterra kernels Marmarelis, 1997, 2004) as methods for detecting states of physiological impairment (Marmarelis et al., 1993; Chon et al., 1994; Mitsis et al., 2006; Zanos et al., 2008; Kang et al., 2015). For example, Kang et al. (2015) showed differential temporal relationships between frontal and occipital electroencephalogram (EEG) signals in Alzheimer's patients compared to healthy individuals. In contrast, Mitsis et al. (2009) modeled CBF responses to BP and P_{ET}CO_{2} fluctuations in healthy individuals using Laguerre-Volterra kernels before and after ganglionic blockade. A significant increase in spectral power of the first-order kernel in the very low frequency range (<0.04 Hz) was detected, suggesting that the autonomic nervous system is modulator of very slow components of CBF regulation. However, that study focused only on the Volterra kernels and did not examine PDMs. Moreover, it has been hypothesized that the spectral characteristics of the global PDMs and the gains of subject-specific associated nonlinear functions (ANFs) may reflect information about specific autoregulatory mechanisms (Marmarelis et al., 2012, 2013a, 2014). For example, the blunted vasomotor reactivity in mild cognitive impaired individuals affected the slope of associated nonlinear function of PDM having low-pass characteristics (Marmarelis et al., 2014). But whether these model-based indices can function as targeted biomarkers of cerebrovascular sympathetic control remains unclear due to the lack of data collected under controlled experimental conditions.

Here we seek to determine whether global PDMs can detect the attenuation of cerebrovascular control induced by targeted functional sympathetic blockade. To achieve this we analyzed spontaneous BP, P_{ET}CO_{2}, and CBF signals recorded from healthy individuals before and after pharmacological α_{1}-adrenergic receptor blockade. The intrinsic nonlinear dynamics of the underlying relations were examined using Laguerre-Volterra based first- and second-order kernels. Further, we extended the kernel approaches by applying the notion of global PDMs that can be passed through ANFs (Marmarelis et al., 2012, 2013a, 2014). Collectively this approach allowed us to explicitly characterize nonlinear dynamics which is not possible using methods like transfer function and wavelet decomposition analysis (Saleem et al., 2016). We hypothesized that (1) very low frequency linear components (defined by the first-order kernels) of the BP and P_{ET}CO_{2} are mutually coupled, and (2) the gains of the ANFs of very low frequency global PDMs for the BP can discriminate between intact and impaired states of sympathetic neurovascular control. Our findings suggest that global PDMs and ANFs gains may have potential utility as functional biomarkers of neurovascular impairment associated with adrenergic dysfunction; a key pathogenic mechanism for poor outcomes in conditions like autonomic failure, stroke (Korpelainen et al., 1999), and traumatic brain injury (Takahashi et al., 2014).

## 2. Data Collection and Processing Methods

All experimental protocols were performed according to the standards set by the Declaration of Helsinki and were approved by the New Zealand Central Regional Ethics Committee [approval number: CEN/09/10/072].

### 2.1. Participants

We performed retrospective analysis of data collected from 18 subjects after getting a written informed consent from each subject to participate in this study. None of the subjects had taken caffeine-containing products and heavy exercise for at least 12 h prior to the study, and all were instructed to have a light breakfast at least 2 h before the study. Only healthy subjects with no regular medication and having no known history of any cardiovascular, respiratory or neurological diseases were recruited for this study.

### 2.2. Measurements

As reported previously (Saleem et al., 2016), the continuous electrocardiogram, non-invasive blood pressure via finger photoplethysmography (Finometer MIDI, MLE1054-V, Finapres Medical Systems, Amsterdam, The Netherlands), right middle cerebral artery blood flow velocity (MCAv; 2 MHz pulsed Doppler ultrasound, ST3 Digital Transcranial Doppler (TCD) System, Spencer Technologies, Seattle, WA), and partial-pressure of end-tidal CO_{2} sampled from a face mask (gas analyzer model ML206, ADInstruments, Colorado Springs, CO, USA) were recorded continuously at 1 kHz per channel via an analog-to-digital converter (PowerLab/16SP ML795; ADInstruments) interfaced with a computer. The recorded signals were stored on a computer for off-line analysis using LabChart 7 Pro (ADInstruments, Colorado Springs, CO, USA).

### 2.3. Experimental Protocols

All experimental protocols were performed with participants lying in the supine position. The experiments were performed between 9 am and 12 noon in a quiet, temperature- and humidity-controlled laboratory (22–23°C). After acclimatization to the laboratory environment, 6 min of baseline resting data were recorded before the administration of placebo pill (9 subjects with 6 females, age: 21–26 years) or 0.05 mg/Kg oral Prazosin (9 subjects with 5 females, age: 21–26 years). To coincide with the peak pharmacodynamic activity on the cerebrovasculature, the post-ingestion protocols were repeated after 120 min which have previously (Jaillon, 1980) been observed to block ≈80% of the pressor response to phenylephrine in healthy normotensive subjects of similar age and body mass index to those of the current study. Participants were free to move around the laboratory, but were not allowed to eat or drink during the post-ingestion period of 120 min.

### 2.4. Data Pre-processing

Beat-to-beat mean arterial BP, mean MCAv values and breath-to-breath P_{ET}CO_{2} time series were determined from the recorded continuous BP, MCAv (an index of CBF) and P_{ET}CO_{2} waveforms. These beat-to-beat signals were then interpolated using cubic spline before decimating to 1 Hz. Occasional measurement artifacts were identified and removed using Grubb's test. Next, a fifth-order Butterworth filter with a cut-off frequency of 0.005 Hz was employed to remove the very slow dynamics and baseline drift of the time series. To remove transient artifacts associated with data filtering, the initial and last segments (equal to the filter length which is 5 in the current study) of data were ignored for subsequent analyses. All data pre-processing and analyses were performed using a custom written software in MATLAB (version R2014b; Mathworks).

### 2.5. Statistics

All values are presented as mean ± standard error, unless otherwise stated. The normality of the data was determined using the Shapiro-Wilk's test. For ease of interpretation and analysis, all values and confidence intervals are presented in standard units. Linear mixed-effects models were employed to determine differences between placebo and sympathetic blockade groups by specifying a main effect for the condition (pre- vs. post-blockade), a main effect for the group (placebo vs. sympathetic) and a condition × group interaction. Statistically significant two-way (condition × group) interactions show that the sympathetic blockade effects were different from the placebo treatment. Paired-*t* tests were performed for *post-hoc* pairwise comparisons to determine the significant difference between pre- and post-blockade responses. All data were pooled for statistical analysis because there were no gender interactions. The significance level for Grubb's test was set at 0.95. Statistical significance was set a priori at *p* < 0.05.

### 2.6. Principal Dynamic Modes Analysis

The intrinsic nonlinear dynamics of BP and MCAv fluctuations have been examined with the use of Laguerre-Volterra kernels based PDM analysis in recent studies (Mitsis et al., 2009; Marmarelis et al., 2012, 2013a). Briefly, the proposed methodology starts with the estimation of first- and second-order (self and cross) Volterra kernels using Laguerre expansions from the given inputs (BP and P_{ET}CO_{2}) and output (MCAv). These first- and second-order self-kernels of all subjects (of both pre- and post-blockade conditions) are combined to form a rectangular matrix that is used to compute the global PDMs via singular value decomposition. The resulting global PDMs form a filter bank where each filter generates the signal (via convolution of each input signal with the respective global PDM) which is subsequently passed through an associated nonlinear function. The intermodulatory interactions among the inputs can be included by computing the cross-kernels and pair-products of PDM outputs (Marmarelis et al., 2013a). Collectively, the polynomial transformed PDM outputs of both inputs and the cross-terms form the system output (MCAv). The mathematical relations of the above described methodology are given below.

The second-order Volterra model for the dual input (BP + P_{ET}CO_{2} − MCAv) system of cerebral hemodynamics (Mitsis et al., 2009) can be written as,

where *P*, *F*, and *C* are BP, MCAv and P_{ET}CO_{2} respectively, *k*_{0} is a zeroth-order kernel (a constant), {*k*_{P}, *k*_{C}} are first-order kernels of *P* and *C* respectively, {*k*_{PP}, *k*_{CC} and *k*_{PC}} respectively are second-order self- and cross-kernels describing the nonlinear interactions of *P* and *C* at time lags (*m*_{1},*m*_{2}), and ϵ(*n*) includes the effects of measurement or modeling errors. *M* is the order of the system memory. It is assumed that M has the same value for each series expansion term for Equation (1).

First- and second-order kernels of the BP and P_{ET}CO_{2} inputs for each subject can be estimated using orthonormal Laguerre functions {*b*_{j}, *j* = 1, …, *L*}. The kernels estimated in this way include exponential decay which is consistent with the dynamic characteristics of the physiological systems (Marmarelis, 2004) and can be written as Kang et al. (2015),

where the Laguerre functions are given by Mitsis and Marmarelis (2002),

with a parameter α which determines the asymptotic exponential decay rate. In this study the same value of α, and the same number of Laguerre functions (i.e., *L*) are employed for estimation of kernels (both BP and P_{ET}CO_{2} inputs) as defined in Equations (2–7), and β = [β_{0} β_{P} β_{PP} β_{C} β_{CC} β_{PC}] are the expansion coefficients estimated using the ordinary least-square technique (Kang et al., 2015) i.e.,

with,

where *I* is a *S* × 1 unit vector (containing 1s) and *V*_{P}, *V*_{C} (each a *S* × *L* matrix) are given by,

and each of *V*_{PP}, *V*_{CC} and *V*_{PC} is a *S* × *L*(*L* + 1)/2 matrix with columns defined by the complete set of *L*(*L* + 1)/2 unique pairs of (*b*_{j1} ∗ *P*) ⊙ (*b*_{j2} ∗ *P*), (*b*_{j1} ∗ *C*) ⊙ (*b*_{j2} ∗ *C*) and (*b*_{j1} ∗ *P*) ⊙ (*b*_{j2} ∗ *C*), respectively, with *j*_{1} = 1, …, *L* and *j*_{2} = 1, …, *j*_{1}. The notation ∗ represents convolution and ⊙ represents the Hadamard product. *S* = *T* − 1 where *T* is the data length of each input (BP or P_{ET}CO_{2}), and *F* is the output MCAv.

The minimum set of basis functions, namely Principal Dynamic Modes (PDMs) (Marmarelis et al., 2012, 2013a), that can represent the BP + P_{ET}CO_{2} to MCAv dynamics adequately is determined via the singular value decomposition of a rectangular matrix containing the estimated first- and second-order (self) kernels of all subjects of both (pre- and post-blockade) conditions of a specific group (i.e., sympathetic blockade or placebo treatment) for each input, given as,

Here,

is a rectangular matrix where ${{k}}_{{1}{,}{N}}^{{c}}$ is a first-order kernel (in the form of a column vector) and ${{k}}_{{2}{,}{N}}^{{c}}$ is a second-order self-kernel (in the form of a block matrix) multiplied by the standard deviation of the respective input [i.e., ${{\sigma}}_{{N}}^{{c}}$ for the *N*th subject for one input (BP or P_{ET}CO_{2})]. The superscripts *c* and *b* correspond to pre- and post-blockade conditions, respectively. For subsequent analyses we have made use of matrix *U* (of equation 12 instead of matrix *V*) to derive global PDMs because the columns of matrix *U* show similar characteristics across subjects. The singular vectors (the columns of matrix *U*) corresponding to the significant singular values (of diagonal matrix λ) form the global PDMs (i.e., *g*_{i,p} or *g*_{i,c} in Equation 14) for the respective input (BP or P_{ET}CO_{2}), and are used to describe the dynamics of the underlying (BP – MCAv or P_{ET}CO_{2} – MCAv) relationships. The resulting global PDMs form a filter bank where each filter takes in the input and provides an output signal which is subsequently passed through a nonlinear function (i.e., *f*_{i,p} or *f*_{i,c} in Equation 14) to form the system output. The output equation for the dual-input PDM based model is given by Marmarelis et al. (2013a, 2014),

where *g*_{i,p}, *g*_{i,c} are the *i*th global PDMs for the BP and P_{ET}CO_{2} inputs respectively, *H* is the number of the global PDMs, and *f*_{i,p}, *f*_{i,c} are the nonlinear functions associated with each global PDM *g*_{i,p}, *g*_{i,c} respectively, termed associated nonlinear functions (ANFs), generally given in a polynomial form as Kang et al. (2015),

The ANFs (defined in Equation 15) capture the intrinsic nonlinearities of the underlying (BP–MCAv and P_{ET}CO_{2}–MCAv) relations. The cross-terms (in Equation 14) are composed of pair products of the PDM outputs and indicate the interactions between BP and P_{ET}CO_{2} fluctuations. In the present study only significant cross-terms are included in the final model. A cross-term is considered significant if its correlation coefficient with the output signal exceeds 99% significance level using the *w*-statistic (derived from the Fisher transform of correlation coefficient) (Marmarelis et al., 2012). The gain of each ANF was determined by the slope of its best linear fit (Kang et al., 2015). This parameter represents the linear gain coefficient of the respective global PDM, and signifies the relative strength of a particular contributing mechanism with a larger slope value reflecting a stronger contribution (Marmarelis et al., 2012, 2013a). These gains were then used as classification parameters to discriminate between pre- and post-blockade conditions.

The optimal values for the Laguerre parameter α, the number of Laguerre functions *L*, the number of global PDMs *H*, and the order of ANFs were selected using Bayesian information criteria (BIC) (Marmarelis et al., 2016). In this study, α = 0.5, *L* = 5 and five global PDMs with cubic ANFs were found to be appropriate for the BP + P_{ET}CO_{2} – MCAv relations for all subjects across both groups (sympathetic blockade and placebo) and experimental conditions (pre- and post-blockade). The shapes of global PDMs remain similar when different (8 out of 9) subjects were selected across one group. The robustness of the proposed method was verified by a leave-one-subject-out cross-validation technique where *N* − 1 (out of *N*) subjects were used to estimate the global PDMs. Next, these estimated global PDMs were applied to the remaining subject, and the ANFs and corresponding gains were determined. This procedure was repeated for all subjects.

The first-order kernels for the BP and P_{ET}CO_{2} inputs {*k*_{P}, *k*_{C}} were also estimated in the frequency-domain. This was done by computing the magnitude of their discrete Fourier transforms and calculating their spectral power by squaring the magnitude of the respective kernel in the frequency-domain. The second-order self-kernels {*k*_{PP}, *k*_{CC}} and cross-kernel {*k*_{PC}} were also determined in the frequency domain by calculating the magnitude of their two-dimensional discrete Fourier transforms. The kernel spectral power for frequencies <0.03 Hz was calculated by averaging the spectral power of the respective kernel over the frequency range <0.03 Hz. The frequency limit of <0.03 Hz was chosen based on the frequency-domain characteristics of the pressure-flow relations under similar blockade conditions (Saleem et al., 2016) to those of the current study.

## 3. Results

### 3.1. Spontaneous Baseline Parameters

The baseline physiological parameters before and after the α_{1}-adrenergic receptor blockade and placebo treatment are shown in Table 1. As reported previously (Saleem et al., 2016), no significant condition × group interactions were found for heart rate, BP, MCAv, P_{ET}CO_{2} and breathing rate, indicating that α_{1}-adrenergic receptor blockade changed none of these baseline parameters as compared with the placebo treatment.

**Table 1. Spontaneous baseline parameters before and after the α _{1}-adrenergic receptor blockade and placebo treatment**.

### 3.2. Laguerre-Volterra Kernel Analysis

The first-order kernel for BP (Figures 1A,B) resembled a high-pass filter under pre-blockade conditions for both sympathetic blockade (Figure 1C) and placebo groups. The α_{1}-adrenergic receptor blockade increased the first-order kernel spectral power significantly (Figure 1D; *p* = 0.00052 for pre- vs. post-blockade paired-*t* test) for the frequencies <0.03 Hz, suggesting that very low frequency BP-MCAv relations become more pronounced after sympathetic blockade. This increase was significantly different from the placebo response (*p* = 0.00023 for condition × group interaction). In contrast, the first-order kernel for P_{ET}CO_{2} resembled a low-pass filter under pre-blockade conditions in both sympathetic blockade (Figure 2C) and placebo groups. The α_{1}-adrenergic blockade did not qualitatively alter the low-pass characteristics but did reduce the magnitude significantly for the frequencies <0.03 Hz (Figure 2D; *p* = 0.0042 for pre- vs. post-blockade paired-*t* test). This decrease was significantly different from the placebo response (*p* = 0.0001 for condition × group interaction). The average changes in spectral powers of both BP and P_{ET}CO_{2} first-order kernels for frequencies <0.03 Hz are summarized in Figures 3A,C for all subjects of both groups.

**Figure 1. First-order kernels of the BP {k_{P}} in the time- (upper panels) and frequency-domain (lower panels) before and after α_{1}-adrenergic receptor blockade, averaged over all subjects (means ± SD)**.

**(A)**time-domain representation before blockade,

**(B)**time-domain representation after blockade,

**(C)**frequency-domain representation before blockade,

**(D)**frequency-domain representation after blockade.

**Figure 2. First-order kernels of the P _{ET}CO_{2} {k_{C}} in the time- (upper panels) and frequency-domain (lower panels) before and after α_{1}-adrenergic receptor blockade, averaged over all subjects (means ± SD)**.

**(A)**time-domain representation before blockade,

**(B)**time-domain representation after blockade,

**(C)**frequency-domain representation before blockade,

**(D)**frequency-domain representation after blockade.

**Figure 3. Changes in spectral powers of (A)** first-order BP kernel {*k*_{P}}, **(B)** second-order BP self-kernel {*k*_{PP}}, **(C)** first-order P_{ET}CO_{2} kernel {*k*_{C}}, **(D)** second-order P_{ET}CO_{2} self-kernel {*k*_{CC}}, and **(E)** cross-kernel between BP and P_{ET}CO_{2} {*k*_{PC}} before and after α_{1}-adrenergic receptor blockade and placebo treatment for frequencies <0.03 Hz. **p* < 0.05 represents a significant two-way (condition × group) interaction between placebo and sympathetic blockade groups.

Second-order self-kernels for BP before and after α_{1}-adrenergic receptor blockade are shown in Figure S1. Average spectral power increased after blockade for frequencies <0.03 Hz but was not significant across all individuals (*p* = 0.39 for pre- vs. post-blockade paired-*t* test; Figure 3B). Second-order self-kernels for the P_{ET}CO_{2}, shown in Figure S2, demonstrated similar trends to the P_{ET}CO_{2} first-order kernels (Figure 2) but the changes between pre- and post-blockade conditions were not significant (*p* = 0.17 for pre- vs. post-blockade paired-*t* test; Figure 3D) across individuals for frequencies <0.03 Hz. No significant condition × group interactions were found for second-order BP and P_{ET}CO_{2} self-kernels between the sympathetic blockade and placebo groups (Figures 3B,D).

The cross-kernels between the BP and P_{ET}CO_{2} before and after sympathetic blockade are shown in Figure S3. The changes in average spectral power associated with sympathetic blockade were not significant (*p* = 0.89 for pre- vs. post-blockade paired-*t* test; Figure 3E) with similar magnitude responses to placebo (*p* = 0.52 for condition × group interaction). We found most individuals did not exhibit significant cross-terms with three exceptions (i.e., two subjects in the sympathetic blockade group, and one subject in the placebo cohort).

Figure 4 shows scatter-plots of the logarithmic spectral power of the first-order BP and P_{ET}CO_{2} kernels for frequencies <0.03 Hz for both placebo and sympathetic blockade groups, and achieves complete separation between pre- and post-blockade conditions of sympathetic blockade group.

**Figure 4. Scatter-plots of logarithmic spectral power of first-order BP {k_{P}} and P_{ET}CO_{2} {k_{C}} kernels for frequencies <0.03 Hz for (A)** placebo treatment,

**(B)**sympathetic blockade. The classification curve in subplot

**(B)**has been obtained using a linear discriminator.

### 3.3. Global PDM Analysis

The estimated five global PDMs from a set of eight subjects (excluding subject 1; leave-one-subject-out cross validation) of sympathetic blockade group are shown in Figure 5 for BP and P_{ET}CO_{2} in the time- and frequency-domain. The first global PDM for BP exhibited high-pass filter characteristics with a blunt spectral peak at ≈0.1 Hz for the sympathetic blockade group. A spectral peak was observed at ≈0.04 Hz for the 2nd PDM, at ≈0.03 Hz for the 3rd, and at ≈0.01 Hz for the 5th PDM. In contrast, the 4th PDM resembled a low-pass filter for BP.

**Figure 5. The five global PDMs for the BP (left panels)** and P_{ET}CO_{2} **(right panels)** of the dual-input (BP + P_{ET}CO_{2} – MCAv) models of sympathetic blockade group in the time- **(upper panels)** and frequency-domain **(lower panels)**.

Similar to BP, the first global PDM for P_{ET}CO_{2} showed high-pass characteristics with a blunt spectral peak at ≈0.07 Hz. The 2nd, 3rd, and 5th global PDMs had spectral peaks at ≈0.04 Hz, ≈0.03 Hz and ≈0.01 Hz, respectively. Also similar to BP, the 4th global PDM for P_{ET}CO_{2} showed low-pass characteristics. The frequency patterns of global PDMs for both BP and P_{ET}CO_{2} were similar between placebo and sympathetic blockade groups.

The ensemble averages of estimated cubic ANFs, using mean coefficient values, along with their standard deviation bounds for sympathetic blockade group under pre-blockade condition are shown in Figure 6 for the BP and P_{ET}CO_{2}. The ANFs of the first global PDMs for BP were found to be almost rectilinear and dominant over other ANFs having curvilinear characteristics. Considerable inter-subject heterogeneity was found across ANFs of BP, especially for the 2nd, 4th, and 5th ANFs. Positive slopes were found for the 1st and 3rd BP ANFs across all subjects while few (3, 6, and 5) subjects had negative slopes for the 2nd, 4th, and 5th BP ANFs (respectively). Considerable variability across ANFs of P_{ET}CO_{2} was found for all subjects with all ANFs having curvilinear characteristics.

**Figure 6. The ensemble averages of estimated ANFs along with their standard deviation bounds for the sympathetic blockade group under pre-blockade conditions for the BP (upper panels)** and P_{ET}CO_{2} **(lower panels)**. The solid lines represent means and dotted lines represent standard deviation bounds.

The estimated cubic ANFs along with their best linear fits for one representative subject of sympathetic blockade group are shown in Figure 7. The ANF of the first global PDM of BP had the highest slope and is almost linear. However, other ANFs had curvilinear characteristics with the 2nd and 3rd showing positive trends, and the 4th and 5th showing negative trends. In contrast to BP, all ANFs associated with P_{ET}CO_{2} had curvilinear characteristics with the 1st, 3rd and 4th having positive trends, and the 2nd and 5th showing negative trends. The ensemble averages of the estimated gains of ANFs for the sympathetic blockade group under pre- and post-blockade conditions are shown in Figure 8. We found no significant changes between pre- and post-blockade conditions across all ANFs.

**Figure 7. The estimated ANFs for subject 1 of sympathetic blockade group corresponding to the five global PDMs (shown in Figure 5) for the BP (upper panels)** and P

_{ET}CO

_{2}

**(lower panels)**, along with their best linear fits (red line).

**Figure 8. The ensemble averages of estimated gains of ANFs for the sympathetic blockade group under pre- and post-blockade conditions for the BP (upper panels)** and P_{ET}CO_{2} **(lower panels)**. Errors bars represent means ± SD.

To determine the relevance of interactions between BP and P_{ET}CO_{2} we examined pair-products for all global PDMs. In general pair-products were not significant across the majority of subjects, so were not included in the final models. However, significant pair-products were found for two subjects in sympathetic blockade group (i.e., the global PDM 1 of BP with the global PDM 2 of P_{ET}CO_{2} in subject 2 before blockade, and the global PDM 2 of BP with the global PDM 2 of P_{ET}CO_{2} in subject 5 after blockade) and across one subject for placebo group (i.e., the global PDM 1 of BP with the global PDM 4 of P_{ET}CO_{2} after placebo treatment in subject 3).

We tested different combinations of ANFs for sympathetic blockade group and found that separation of pre- and post-blockade conditions could be attained for all subjects by using the gains of ANFs associated with the 3rd and 4th global PDMs for BP (Figure 9B). However, no clear separation was evident for the pre- and post-blockade conditions of placebo treatment as shown in Figure 9A. This observation suggests that sympathetic nervous system may be linked with the BP global PDMs that have low-pass (4th global PDM) and ≈0.03 Hz resonant peak (3th global PDM) characteristics.

**Figure 9. Scatter-plots of estimated gains of ANFs corresponding to the 3th and 4th global PDMs for the BP for pre- and post-blockade conditions of (A)** placebo treatment, and **(B)** sympathetic blockade group. The classification curve in subplot **(B)** has been obtained using a linear discriminator. Leave-one-subject-out cross-validation technique was adopted to derive the global PDMs and the corresponding ANFs.

## 4. Discussion

### 4.1. Main Findings

This study sought to identify vascular biomarkers that could delineate between states of intact and impaired cerebral autoregulatory control. The proposed approach examines CBF dynamics in response to spontaneous BP and P_{ET}CO_{2} fluctuations using Laguerre-Volterra kernels based PDM analysis. This methodology makes use of global PDMs derived from the first- and second-order kernels, and the linear gains of their ANFs. We found that the mutual coupling of linear components (i.e., first-order kernels) of the BP and P_{ET}CO_{2} for frequencies <0.03 Hz could detect the cerebrovascular effects of α_{1}-adrenergic blockade. Moreover, in keeping with our hypothesis, the gains of ANFs associated with the very low frequency global PDMs (having low-pass and ≈0.03 Hz peak characteristics) for BP could clearly discriminate between the baseline and blockade states. Collectively these findings suggest that global PDMs and ANFs may have potential utility as biomarkers of sympathetic cerebrovascular dysfunction.

### 4.2. Modeling with Laguerre-Volterra Kernels and Global PDMs

Volterra formulation of kernels represents the system output in terms of linear and nonlinear interactions between different input(s), as well as between lagged epochs within a given set of input(s). The linear interactions are quantified by the first-order kernel representing the linear component (similar to the impulse response) of the system (Panerai et al., 1999; Marmarelis et al., 2012; Angarita-Jaimes et al., 2014). In contrast, the nonlinear interactions are summarized by second-order kernels, which reflect the interactions between different epochs within a respective input (i.e., self-kernels) or interactions among the different epochs of two inputs (i.e., cross-kernels) (Marmarelis et al., 2012). The spectral representation of the second-order kernels reflects the nonlinear interactions between frequency-dependent dynamics of the respective input(s) (Mitsis et al., 2002). The diagonal entries (of the frequency spectrum of the second-order kernel) reflect nonlinear physiological processes operating within the respective frequencies while the off-diagonal entries reflect nonlinear intermodulatory interactions among different mechanisms that operate at different frequencies (Mitsis and Marmarelis, 2002; Mitsis et al., 2002).

In this study a fusing procedure, via singular value decomposition, was applied to the estimated Laguerre-Volterra kernels from an ensemble of subjects to obtain a minimal set of basis functions. As these global PDMs are common to all subjects, they form a compact representation of individual's system dynamics (Marmarelis et al., 2013b; Angarita-Jaimes et al., 2014). To complete the nonlinear PDM-based model, the static nonlinearities associated with each global PDM are captured in the form of subject-specific ANFs which can be used to distinguish the dynamics of individual subjects. Recent studies (Angarita-Jaimes et al., 2014) have shown that the compact formulation of “global” PDMs reduces the inter-subject variability and the subject-specific ANFs capture the intra-subject variance across repeated recordings.

The robustness of the global PDMs (determined by the leave-one-subject-out cross-validation) shows that the apparent shape of global PDMs does not change when a different combination of participants (i.e., 8 out of 9) is chosen for one group. Our findings are consistent with previous studies (Marmarelis et al., 2013a) showing that the characteristic features of global PDMs do not vary appreciably when different subsets of reference subjects were chosen to extract global PDMs. Although global PDMs are common for all subjects of both groups, the ANFs associated with each global PDM are subject-specific and have the ability to uniquely characterize the dynamics of each subject.

### 4.3. Physiological Interpretation of Laguerre-Volterra Kernels and Global PDMs

In general the physiological origins of the observed principal dynamic modes are not well understood but several plausible explanations should be considered. The linear dynamics of BP represented by its first-order kernel have high-pass characteristics. This indicates that slow BP fluctuations are more effectively attenuated than higher frequency BP fluctuations and a plausible biophysical model that could account for these dynamics is the arterial Windkessel (Zhang et al., 2009). In contrast, the first-order kernel for P_{ET}CO_{2} has low-pass characteristics suggesting that very slow P_{ET}CO_{2} fluctuations have a pronounced effect on MCAv dynamics. These observations are consistent with studies showing that CO_{2} modulation of cerebrovasculature occurs with a slow dynamic time constant under spontaneous baseline conditions (Ainslie and Duffin, 2009; Peebles et al., 2012). Most of the spectral power of the second-order (self and cross) kernels are concentrated at frequencies <0.03 Hz suggesting that significant intermodulatory interactions between BP and P_{ET}CO_{2} fluctuations are found across very low frequencies (<0.03 Hz). In keeping with this interpretation, our ability to separate between pre- and post-blockade states of α_{1}-adrenergic group was optimal when we considered both BP and P_{ET}CO_{2} kernels for frequencies <0.03 Hz using a two-dimensional separation regime (Figure 4B).

It has previously been suggested that the distinct spectral characteristics of global PDMs may carry functional information about specific physiological mechanisms (Marmarelis et al., 2012, 2013a). Our analyses show that the first global PDM for BP has high-pass characteristics having spectral peak at ≈0.1 Hz. Similar to the first-order kernel, we speculate that this feature may relate to the Windkessel properties of the vasculature. The large positive value of linear gain coefficient of its ANF indicates that this component may be an influential determinant of CBF dynamics. In contrast, the 4th global PDM has low-pass characteristics with a cut-off frequency of ≈0.025 Hz indicating the integration of the BP input over several seconds. Speculatively, this feature may reflect the activity of neurogenic or humoral blood pressure control mechanisms that operate with a slow dynamic time constant (Ursino and Magosso, 2004). The 2nd, 3rd, and 5th global PDMs have complex spectral patterns with spectral peaks at ≈0.04, ≈0.03, and ≈0.01 Hz, respectively. The physiological origins of these components are unclear but our data suggests that sympathetic cerebrovascular control may contribute to the 3rd and 4th global PDMs since their linear gain coefficients were affected following sympathetic blockade.

The global PDMs for P_{ET}CO_{2} exhibit spectral characteristics that were similar to BP with the first global PDM also showing high-pass characteristics with a spectral peak at ≈0.07 Hz. The 4th global PDM shows low-pass characteristics with a half-max frequency of ≈0.01 Hz indicating the integration of the P_{ET}CO_{2} over several seconds. It has been suggested that this component might be sensitive under hypo- or hyper-ventilation conditions (Marmarelis et al., 2012). Consistent with one previous study (Mitsis et al., 2006), we observed that a significant portion of the P_{ET}CO_{2} power resides within the very-low frequencies of this low-pass global PDM. Similar to BP, the 2nd, 3rd, and 5th global PDMs also have complex spectral patterns with resonant peaks at ≈0.04, ≈0.03, and ≈0.01 Hz, respectively. However, in contrast to BP, α_{1}-adrenergic blockade did not affect these components suggesting that P_{ET}CO_{2} global PDMs may not be mediated by sympathetic activity.

An interesting observation was that the resonant peaks of BP and P_{ET}CO_{2} global PDMs occurred at common frequencies. This might occur as a result of BP and P_{ET}CO_{2} fluctuations engaging common vascular control mechanisms. Though the exact nature of these interactions is not well understood, the general absence of significant cross-terms between BP and P_{ET}CO_{2} shows that their effects on cerebrovasculature might be independent of each other (Marmarelis et al., 2012, 2013a). It is also important to note that significant changes after placebo treatment occurred in a few subjects. This observation suggests that there are time dependent effects in cerebral hemodynamics consistent with a non-stationary system. Such changes highlight the importance of including placebo group in the design of experimental studies involving multiple measurements across time.

### 4.4. Implications

The findings of the present study have important methodological implications. The spectral assessment of cerebral hemodynamics is often performed within predefined continuous frequency ranges that were originally based on simple transfer function analysis (Zhang et al., 1998). However, the current literature lacks the experimental verification of these bands. There is growing evidence that CBF dynamics within the VLF range is an amalgam of many redundant physiological processes that individually may only contribute to a sub-component of this predefined frequency range. Indeed, in a recent study using wavelet decomposition analysis we showed that sympathetic blockade altered the phase synchronization between BP and MCAv within the 0.02–0.03 Hz range (Saleem et al., 2016). Here we demonstrate that within this range, α_{1}-adrenergic blockade affected only global PDMs having low-pass and ≈0.03 Hz resonant peak characteristics. These findings also prompt the development and potential application of filter banks that can characterize and quantify cerebral haemodynamics along a wide array of narrow-band frequency scales.

### 4.5. Methodological Considerations

Several experimental and methodological limitations need to be considered for interpretation of the findings of this study. First, CBF dynamics were assessed using TCD, which estimates intracranial blood velocity in the middle cerebral artery instead of the volumetric flow. However, intracranial blood velocity can be used as an adequate surrogate of actual blood flow only when the diameter of the insonated vessel remains unchanged under different haemodynamic conditions. This assumption is generally met under normocapnic conditions (Serrador et al., 2000; Coverdale et al., 2014) but other complementary monitoring techniques such as near-infrared spectroscopy can be used in future studies to provide surrogates of blood volume and cortical oxygenation.

Second, we acknowledge that the duration of recorded time series can affect the model estimation. However, the 6 min recordings examined in the present study were of comparable duration to those of previous studies employing Laguerre-Volterra kernels and global PDM-based models for cerebral haemodynamic characterization (Mitsis et al., 2009; Marmarelis et al., 2012, 2013a). The adopted methodology in the present study requires 66 free parameters from a data record of 6 min, satisfying the conventional rule of thumb that the number of independent data samples should be ≈5 times as many as the estimated free parameters to avoid overfitting. It has been observed that the number of free parameters can be reduced to 41 if a Laguerre-Volterra Network having five hidden units is adopted (Marmarelis, 2004). Moreover, the selection of different values of α, and the number of Laguerre functions (i.e., *L*) for each input have been reported resulting in (1) the higher value of α (for example, 0.85 in Marmarelis et al., 2016), and (2) the less number of Laguerre functions for the P_{ET}CO_{2} input than that of the BP input. The increased α and reduced *L* value lead to longer kernel memory and less number of free parameters in the model, respectively, thereby improving the estimation accuracy of the kernels from short recorded data lengths.

Third, the proposed functional biomarkers have been developed based on an open-loop model of cerebral circulation with BP and PCO_{2} as driving inputs and CBF as an output. However, recent studies suggest that there may also be reverse information flow from CBF to BP, suggesting that cerebral haemodynamic characterization may be further improved using closed-loop models (Marmarelis et al., 2013a). In relation to this point, since CBF regulation is a multifaceted process mediated by various physiological variables, future studies should also consider the role of additional inputs to the cerebrovascular system such as cardiac output (Marmarelis et al., 2016). Finally, the present study showed that α_{1}-adrenergic blockade altered the 3rd and 4th global PDMs suggesting that sympathetic activity may underlie these global PDMs. However, our findings do not preclude the possibility that other mechanisms (e.g., cholinergic, myogenic, metabolic factors) may also be contributing to these PDM components.

## 5. Conclusion

In summary, very low frequency (<0.03 Hz) dynamics of CBF regulation in response to spontaneous BP and P_{ET}CO_{2} fluctuations are associated with the human sympathetic modulation. Moreover, the linear gains of ANFs of very low frequency (<0.03 Hz) global PDMs may be used as autoregulatory biomarkers to determine the impairment of human sympathetic control.

## Data Availability

The data can be made available to the interested readers upon request.

## Author Contributions

SS, YT, WK, and PT contributed to study design. YT contributed to data acquisition. SS, YT, WK, and PT contributed to data analysis. SS, YT, WK, and PT contributed to data interpretation. SS, YT, WK, and PT contributed to preparation of manuscript. All authors gave final approval for publication.

## Funding

This research was supported in part by funding from the New Zealand Research Health Council (Ref: 11/125; to YT).

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fphys.2016.00685/full#supplementary-material

## References

Ainslie, P. N., and Duffin, J. (2009). Integration of cerebrovascular CO_{2} reactivity and chemoreflex control of breathing: mechanisms of regulation, measurement, and interpretation. *Am. J. Physiol. Regul. Integr. Compar. Physiol.* 296, R1473–R1495. doi: 10.1152/ajpregu.91008.2008

Angarita-Jaimes, N., Kouchakpour, H., Liu, J., Panerai, R. B., and Simpson, D. M. (2014). Optimising the assessment of cerebral autoregulation from black box models. *Med. Eng. Phys.* 36, 607–612. doi: 10.1016/j.medengphy.2013.12.012

Chon, K. H., Chen, Y. M., Marmarelis, V. Z., Marsh, D. J., and Holstein-Rathlou, N. H. (1994). Detection of interactions between myogenic and TGF mechanisms using nonlinear analysis. *Am. J. Physiol. Renal Physiol.* 267, F160–F173.

Coverdale, N. S., Gati, J. S., Opalevych, O., Perrotta, A., and Shoemaker, J. K. (2014). Cerebral blood flow velocity underestimates cerebral blood flow during modest hypercapnia and hypocapnia. *J. Appl. Physiol.* 117, 1090–1096. doi: 10.1152/japplphysiol.00285.2014

Hamner, J. W., Tan, C. O., Lee, K., Cohen, M. A., and Taylor, J. A. (2010). Sympathetic control of the cerebral vasculature in humans. *Stroke* 41, 102–109. doi: 10.1161/STROKEAHA.109.557132

Jaillon, P. (1980). Clinical pharmacokinetics of prazosin. *Clin. Pharmacokinet.* 5, 365–376. doi: 10.2165/00003088-198005040-00004

Kang, Y., Escudero, J., Shin, D., Ifeachor, E., and Marmarelis, V. (2015). Principal dynamic mode analysis of EEG data for assisting the diagnosis of Alzheimer's disease. *IEEE J. Trans. Eng. Health Med.* 3, 1–10. doi: 10.1109/JTEHM.2015.2401005

Korpelainen, J. T., Sotaniemi, K. A., and Myllylä, V. V. (1999). Autonomic nervous system disorders in stroke. *Clin. Auton. Res.* 9, 325–333. doi: 10.1007/BF02318379

Latka, M., Turalska, M., Glaubic-Latka, M., Kolodziej, W., Latka, D., and West, B. J. (2005). Phase dynamics in cerebral autoregulation. *Am. J. Physiol. Heart Circul. Physiol.* 289, H2272–H2279. doi: 10.1152/ajpheart.01307.2004

Marmarelis, V. Z., Mitsis, G. D., Shin, D. C., and Zhang, R. (2016). Multiple-input nonlinear modelling of cerebral haemodynamics using spontaneous arterial blood pressure, end-tidal CO_{2} and heart rate measurements. *Philos. Trans. R. Soc. A* 374:20150180. doi: 10.1098/rsta.2015.0180

Marmarelis, V. Z., Shin, D. C., Orme, M. E., and Zhang, R. (2013a). Closed-loop dynamic modeling of cerebral hemodynamics. *Ann. Biomed. Eng.* 41, 1029–1048. doi: 10.1007/s10439-012-0736-8

Marmarelis, V. Z., Shin, D. C., Orme, M. E., and Zhang, R. (2014). Model-based physiomarkers of cerebral hemodynamics in patients with mild cognitive impairment. *Med. Eng. Phys.* 36, 628–637. doi: 10.1016/j.medengphy.2014.02.025

Marmarelis, V., Shin, D., and Zhang, R. (2012). Linear and nonlinear modeling of cerebral flow autoregulation using principal dynamic modes. *Open Biomed. Eng. J.* 6:42. doi: 10.2174/1874230001206010042

Marmarelis, V. Z. (1997). Modeling methology for nonlinear physiological systems. *Ann. Biomed. Eng.* 25, 239–251. doi: 10.1007/BF02648038

Marmarelis, V. Z. (2004). *Nonlinear Dynamic Modeling of Physiological Systems*, Vol. 10. Hoboken, NJ: John Wiley & Sons.

Marmarelis, V. Z., Chon, K. H., Chen, Y. M., Marsh, D. J., and Holstein-Rathlou, N. H. (1993). Nonlinear analysis of renal autoregulation under broadband forcing conditions. *Ann. Biomed. Eng.* 21, 591–603. doi: 10.1007/BF02368640

Marmarelis, V. Z., Shin, D. C., Song, D., Hampson, R. E., Deadwyler, S. A., and Berger, T. W. (2013b). Nonlinear modeling of dynamic interactions within neuronal ensembles using principal dynamic modes. *J. Comput. Neurosci.* 34, 73–87. doi: 10.1007/s10827-012-0407-7

Mitsis, G. D., and Marmarelis, V. Z. (2002). Modeling of nonlinear physiological systems with fast and slow dynamics. I. Methodology. *Ann. Biomed. Eng.* 30, 272–281. doi: 10.1114/1.1458591

Mitsis, G. D., Zhang, R., Levine, B. D., and Marmarelis, V. (2002). Modeling of nonlinear physiological systems with fast and slow dynamics. II. Application to cerebral autoregulation. *Ann. Biomed. Eng.* 30, 555–565. doi: 10.1114/1.1477448

Mitsis, G. D., Ainslie, P. N., Poulin, M. J., Robbins, P. A., and Marmarelis, V. Z. (2005). “Nonlinear modeling of the dynamic effects of arterial pressure and blood gas variations on cerebral blood flow in healthy humans,” in *Post-genomic Perspectives in Modeling and Control of Breathing*, eds J. Champagnat, M. Denavit-Saubié, G. Fortin, A. S. Foutz, and M., Thoby-Brisson, (New York, NY: Springer), 259–265.

Mitsis, G. D., Poulin, M. J., Robbins, P. A., and Marmarelis, V. Z. (2004). Nonlinear modeling of the dynamic effects of arterial pressure and CO_{2} variations on cerebral blood flow in healthy humans. *IEEE Trans. Biomed. Eng.* 51, 1932–1943. doi: 10.1109/TBME.2004.834272

Mitsis, G. D., Zhang, R., Levine, B. D., and Marmarelis, V. Z. (2006). Cerebral hemodynamics during orthostatic stress assessed by nonlinear modeling. *J. Appl. Physiol.* 101, 354–366. doi: 10.1152/japplphysiol.00548.2005

Mitsis, G. D., Zhang, R., Levine, B. D., Tzanalaridou, E., Katritsis, D. G., and Marmarelis, V. Z. (2009). Autonomic neural control of cerebral hemodynamics. *IEEE Eng. Med. Biol. Magaz.* 28, 54–62. doi: 10.1109/MEMB.2009.934908

Panerai, R. B., Dawson, S. L., and Potter, J. F. (1999). Linear and nonlinear analysis of human dynamic cerebral autoregulation. *Am. J. Physiol. Heart Circul. Physiol.* 277, H1089–H1099.

Peebles, K. C., Ball, O. G., MacRae, B. A., Horsman, H. M., and Tzeng, Y. C. (2012). Sympathetic regulation of the human cerebrovascular response to carbon dioxide. *J. Appl. Physiol.* 113, 700–706. doi: 10.1152/japplphysiol.00614.2012

Peng, T., Rowley, A. B., Ainslie, P. N., Poulin, M. J., and Payne, S. J. (2008). Multivariate system identification for cerebral autoregulation. *Ann. Biomed. Eng.* 36, 308–320. doi: 10.1007/s10439-007-9412-9

Saleem, S., Teal, P. D., Kleijn, W. B., Ainslie, P. N., and Tzeng, Y. C. (2016). Identification of human sympathetic neurovascular control using multivariate wavelet decomposition analysis. *Am. J. Physiol. Heart Circul. Physiol.* 311, H837–H848. doi: 10.1152/ajpheart.00254.2016

Serrador, J. M., Picot, P. A., Rutt, B. K., Shoemaker, J. K., and Bondar, R. L. (2000). MRI measures of middle cerebral artery diameter in conscious humans during simulated orthostasis. *Stroke* 31, 1672–1678. doi: 10.1161/01.STR.31.7.1672

Takahashi, C., Hinson, H., and Baguley, I. J. (2014). Autonomic dysfunction syndromes after acute brain injury. *Handb. Clin. Neurol.* 128, 539–551. doi: 10.1016/B978-0-444-63521-1.00034-0

Ursino, M., and Magosso, E. (2004). Interaction among humoral and neurogenic mechanisms in ventilation control during exercise. *Ann. Biomed. Eng.* 32, 1286–1299. doi: 10.1114/B:ABME.0000039362.26850.17

Zanos, T. P., Courellis, S. H., Berger, T. W., Hampson, R. E., Deadwyler, S. A., and Marmarelis, V. Z. (2008). Nonlinear modeling of causal interrelationships in neuronal ensembles. *IEEE Trans. Neural Syst. Rehabil. Eng.* 16, 336–352. doi: 10.1109/TNSRE.2008.926716

Zhang, R., Behbehani, K., and Levine, B. D. (2009). Dynamic pressure–flow relationship of the cerebral circulation during acute increase in arterial pressure. *J. Physiol.* 587, 2567–2577. doi: 10.1113/jphysiol.2008.168302

Keywords: biomarkers, blood flow, blood pressure, kernels, principal dynamic modes

Citation: Saleem S, Tzeng Y-C, Kleijn WB and Teal PD (2017) Detection of Impaired Sympathetic Cerebrovascular Control Using Functional Biomarkers Based on Principal Dynamic Mode Analysis. *Front. Physiol*. 7:685. doi: 10.3389/fphys.2016.00685

Received: 18 September 2016; Accepted: 23 December 2016;

Published: 09 January 2017.

Edited by:

Zbigniew R. Struzik, University of Tokyo, JapanReviewed by:

Can Ozan Tan, Harvard Medical School, USAVasilis Z. Marmarelis, University of Southern California, USA

Stephen John Payne, University of Oxford, UK

Copyright © 2017 Saleem, Tzeng, Kleijn and Teal. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Paul D. Teal, paul.teal@vuw.ac.nz