## Abstract

The linear dynamic relationship between systemic arterial blood pressure (ABP) and cerebral blood flow velocity (CBFV) was studied by time- and frequency-domain analysis methods. A nonlinear moving-average approach was also implemented using Volterra-Wiener kernels. In 47 normal subjects, ABP was measured with Finapres and CBFV was recorded with Doppler ultrasound in both middle cerebral arteries at rest in the supine position and also during ABP drops induced by the sudden deflation of thigh cuffs. Impulse response functions estimated by Fourier transfer function analysis, a second-order mathematical model proposed by Tiecks, and the linear kernel of the Volterra-Wiener moving-average representation provided reconstructed velocity model responses, for the same segment of data, with significant correlations to CBFV recordings corresponding to*r* = 0.52 ± 0.19, 0.53 ± 0.16, and 0.67 ± 0.12 (mean ± SD), respectively. The correlation coefficient for the linear plus quadratic kernels was 0.82 ± 0.08, significantly superior to that for the linear models (*P* < 10^{−6}). The supine linear impulse responses were also used to predict the velocity transient of a different baseline segment of data and of the thigh cuff velocity response with significant correlations. In both cases, the three linear methods provided equivalent model performances, but the correlation coefficient for the nonlinear model dropped to 0.26 ± 0.25 for the baseline test set of data and to 0.21 ± 0.42 for the thigh cuff data. Whereas it is possible to model dynamic cerebral autoregulation in humans with different linear methods, in the supine position a second-order nonlinear component contributes significantly to improve model accuracy for the same segment of data used to estimate model parameters, but it cannot be automatically extended to represent the nonlinear component of velocity responses of different segments of data or transient changes induced by the thigh cuff test.

- cerebral blood flow
- mathematical model
- thigh cuff test
- Volterra-Wiener kernels
- cerebral hemodynamics

the temporal relationship between beat-to-beat changes in mean arterial blood pressure (MABP) and mean cerebral blood flow (MCBF) is modulated by the mechanism of cerebral autoregulation, which tends to maintain MCBF relatively constant despite changes in systemic MABP. In animals, sudden changes in ABP are transmitted to the cerebral circulation, inducing similar changes in CBF, but under normal conditions the CBF tends to return to its original value with a time constant of a few seconds (13, 26, 28, 32). A similar transient response of CBF to ABP disturbances has been observed in humans, and it is usually referred to as dynamic cerebral autoregulation (1, 12, 19). The recent advent of transcranial Doppler ultrasound (TCD) has been instrumental in the study of the dynamic relation between changes in systemic ABP and CBF because it has the temporal resolution required to detect fast changes in CBF, unlike more traditional methods of measuring CBF. Although the Doppler technique measures blood velocity rather than absolute flow, these two quantities will change in parallel unless the diameter of the insonated vessel changes significantly, which has not been shown to occur (11, 18).

Studies of dynamic cerebral autoregulation have used different techniques to induce rapid changes in ABP. These include the sudden deflation of thigh cuffs (1, 29), Valsalva maneuvers (12, 30), forced breathing (8), periodic squatting (4), or tilting (2). Other investigators have relied on spontaneous fluctuations in ABP to observe the corresponding transient changes in CBF velocity (CBFV) (5, 9, 10,14, 17, 21-24, 33).

Dynamic cerebral autoregulation has been shown to be disturbed by a number of pathophysiological conditions (1, 2, 4, 8, 10, 21, 23, 24,29) and to correlate well with assessments of static autoregulation (21, 29). Both physiological and clinical studies can benefit from the dynamic approach compared with the more classic view of cerebral autoregulation as a static phenomenon (25, 29). The time course of the CBF transient response can shed light on the structure and efficiency of the feedback mechanisms responsible for CBF control (9, 22, 23, 33). With the use of noninvasive methods to estimate ABP and CBF, it is possible to conduct repeated tests, thus facilitating studies of the influence of other variables or physiological interventions (e.g., posture, blood gases, drugs) on CBF regulation. Modeling the dynamic pressure-flow (or velocity) relationship represents an important element of this research as a conceptual and quantitative tool to refine our understanding and assessment of dynamic cerebral autoregulation for physiological studies and clinical applications.

For step changes in ABP, as obtained by the sudden deflation of thigh cuffs, the CBFV response has been modelled by a second-order differential equation with a set of parameters that can be used to grade the performance of autoregulation (29). The typical CBFV transient that follows the negative step change in ABP and the phase changes observed during forced breathing have also prompted analogies between the dynamic autoregulatory response and the behavior of a simple high-pass filter (8, 9, 33). Empirical modeling has also been performed by means of frequency-domain analysis through calculation of the coherence, transfer function, and phase-frequency response between CBFV and ABP (5, 9, 10, 14, 17, 21-23, 33). Some studies have also obtained estimates of the CBFV impulse response, which reflects the dynamic CBFV response to an impulselike disturbance in ABP (22, 23,33). This approach is quite attractive because the time integral of the impulse response represents the CBFV step response, which is approximately the same transient generated by the sudden deflation of thigh cuffs (1, 29). We have recently demonstrated that the model proposed by Tiecks et al. (29) can also be used to estimate the impulse response of dynamic autoregulation during baseline spontaneous fluctuations in ABP (24).

All of these models rely on the assumption that the dynamic autoregulatory mechanism can be approximated by a linear system. This assumption can be questioned on theoretical grounds, because changes in pressure lead to changes in cerebrovascular resistance. Because resistance is in the denominator of Pouseuille’s law, and because it depends on the fourth power of arteriolar diameter, which may be influenced by myogenic, neurogenic, endothelial, or metabolic mechanisms (25), it becomes clear that there are several nonlinearities that cannot be ignored (31). On the other hand, the assumption of linear systems has been justified for the case of spontaneous fluctuations in ABP, for example, because of the relatively small changes in ABP compared with the larger pressure drops produced by the thigh cuff technique. The objective of the present study is to examine the limitations of the linear assumption and to compare the relative performance of different modeling options. For this purpose we have adopted the approach pioneered by Marmarelis and co-workers (6, 7) in studies of the dynamic renal autoregulation of rats, showing that the renal regulatory mechanism cannot be represented by a linear system. Their analysis was based on the Volterra-Wiener representation of dynamic nonlinear systems (15, 16), which we have also used to obtain linear and quadratic models for the dynamic relationship between beat-to-beat changes of ABP and CBFV in a large group of normal subjects.

## METHODS

#### Subjects and measurements.

Subjects recruited were free from cardiovascular disease as determined from history, full physical examination, and a 12-lead surface electrocardiogram (ECG). None was taking any medication known to affect the autonomic nervous or cardiovascular system. All studies were conducted in the morning, at least 2 h after a light breakfast, and the subjects were asked to refrain from alcohol-, nicotine-, or caffeine-containing products for a minimum of 12 h. Measurements were performed in a quiet, dimly lit room at constant ambient temperature (23°C) with subjects lying supine on a couch with the head supported by two pillows and the right arm supported at atrial height. Arterial blood pressure was measured with a noninvasive finger blood pressure monitor (Finapres 2300, Ohmeda) that has been extensively validated against intravascular measurements (20, 27) and has also been used in previous studies of dynamic autoregulation (5, 14, 24, 33). A finger cuff of the appropriate size was attached to the middle finger of the right hand. During the recording the servo-adjust mechanism of the Finapres was disabled, and an ABP calibration signal was recorded before each measurement. An electrocardiogram (ECG) was obtained using three standard surface chest leads, and a transcutaneous gas monitor (TINA, Radiometer) was attached to record carbon dioxide partial pressure. A dual-channel TCD (SciMed QVL 120) was used to measure CBFV in the right and left middle cerebral arteries (MCA) using 2-MHz probes held immobile by a purpose-built head frame.

The ABP, dual-channel CBFV Doppler-shift signals, and ECG were continuously recorded on digital tape for posterior analysis (DAT, Sony PC 108M). After a 30-min rest period, signals were recorded for two 5-min periods with subjects breathing normally at rest. To perform the thigh cuff test, wide blood pressure cuffs were attached to both legs and inflated 40 mmHg above systolic pressure for 45 s. The cuffs were rapidly deflated, and signals were recorded for another 2 min. Only one thigh cuff test was performed on each subject. Tests were rejected if the ABP drop was <10 mmHg or if it was not accompanied by a simultaneous drop in CBFV in both channels.

#### Data analysis.

The DAT recording was downloaded to a microcomputer in real time. A fast Fourier transform (FFT) was used to extract the maximum frequency velocity envelope, using a time window of 6.25 ms. The ABP and ECG signals were sampled at 200 samples/s per channel.

The ABP signal was calibrated at the start of each recording, and all signals were visually inspected for artifacts or noise. Narrow spikes on the CBFV signals were removed by linear interpolation, and the four signals were low-pass filtered with a zero-phase, eighth-order, Butterworth digital filter with a cutoff frequency of 20 Hz. The beginning and end of each cardiac cycle were detected from the upstroke of the arterial pulse-pressure wave, and the MABP and mean CBFV (MCBFV) were calculated for each cycle. Ectopics occurring during the thigh cuff test led to the rejection of the data as did more than one ectopic per 30 s during the baseline recording. Occasional ectopics could be marked and removed by linear interpolation. Spline interpolation was used to resample the data at 0.2 s to create a uniform time base. For baseline recordings, the MABP and MCBFV time series were normalized by the mean value during the 5-min period and subtracted by 1. The resulting zero-mean signals, reflecting relative changes in ABP and CBFV, are represented by P(*t*) and*V*(*t*), respectively, where *t* is time. A similar procedure was adopted for the thigh cuff, but in this case the reference values were the mean values of MABP and MCBFV preceding the cuff deflation.

Four linear models and one nonlinear (quadratic) model were analyzed by assuming P(*t*) as the input and*V*(*t*) as the output. The mathematical formulation for the different models is given in the
. The simplest linear model, labeled “zero-order model,” corresponds to the approximation*V*(*t*) = P(*t*). Because both signals are normalized in relative units, it is assumed that all velocity changes are purely reflecting the pressure transients in the absence of autoregulation. The corresponding impulse response for this model is*i*
_{Z}(*t*) = 1 for *t* = 0 and*i*
_{Z}(*t*) = 0 for all other values of *t*. Although the zero-order model can be expected to perform poorly, it fits the purpose of providing a “background” condition against which all other models can be compared. The model proposed by Tiecks et al. (29) was implemented using exactly the same set of equations and parameters proposed by the authors (see
). Different combinations of the model parameters (*T*,*D*, and*K*, defined in
) were selected corresponding to 10 grades of autoregulation expressed by a dynamic autoregulation index (ARI) ranging from 0 (absence of autoregulation) to 9 (best autoregulation). The Tiecks model was fitted to baseline data by selecting the value of ARI leading to the minimum quadratic error between the model-generated velocity [*V*
_{T}(*t*)] and*V*(*t*). The corresponding impulse response [*i*
_{T}(*t*)] was calculated as the first derivative of*V*
_{T}(*t*).

Another impulse response estimate was obtained using the FFT approach as described previously (23). The P(*t*) and*V*(*t*) time series were divided into 4 segments with 256 samples (51.2 s) each, multiplied by a cosine-tapered window, and transformed with the FFT algorithm using 30% superposition of segments (Welch method). The corresponding frequency resolution is 0.019 Hz per harmonic. The auto- and cross-spectra, the coherence function, and the transfer function (see
) were estimated using the average of the 4 segments of data and smoothed with a 3-point triangular window, leading to spectral estimates with 18 degrees of freedom per harmonic (3). The phase-frequency response was estimated without unwrapping (23). The inverse FFT was applied to obtain the time-domain impulse response [*i*
_{F}(*t*)] after low-pass filtering the transfer function with a cutoff frequency of 1 Hz.

The Volterra-Wiener approach to estimate the linear and nonlinear representation of the dynamic P(*t*)-*V*(*t*) relationship was implemented as proposed by Marmarelis (15) with the use of Fortran software specially written for this purpose. We have taken advantage of the orthogonal property of Laguerre polynomials to expand the first- and second-order (quadratic) kernels, thus obtaining a more efficient and accurate estimation of these models (15). This method will be referred to as the Wiener-Laguerre representation for nonlinear systems. The software Lysis 6.2, distributed by the Biomedical Simulations Resource, University of Southern California, was used as a benchmark to test and validate our own software. As described in the
, the Volterra series contains a term for each order considered. Because both P(*t*) and*V*(*t*) have zero mean, the constant term was not estimated. Only the second term,*k*
_{1}(*t*), was used to obtain the impulse response corresponding to the linear assumption. To obtain a nonlinear, quadratic model, both the first- and second-order terms were considered (see
). The Wiener formulation assumes that the input signal is random, but, as demonstrated by Marmarelis (15, 16), this condition can be relaxed when Laguerre polynomials are used to expand the different kernels. To increase the whiteness of P(*t*), we decimated both input and output signals by a factor of five, taking every fifth sample to increase the sampling interval to 1 s. Decimation reduces the frequency bandwidth of the signals to 0.5 Hz (3), which is appropriate to describe the fastest components of cerebral autoregulation (5, 14,23, 33). To estimate the linear and quadratic kernels, we used segments of data with 160-s duration, which is approximately the same signal duration used for obtaining the FFT impulse response,*i*
_{F}(*t*). The Laguerre expansion used 10 terms with a value of α = 0.2 (see
). The number of lags used for both the linear and nonlinear kernels was 32.

All of the analyses were conducted separately for the CBFV signals from the right and left MCA. Separate average curves of*i*
_{T}(*t*),*i*
_{F}(*t*), and*k*
_{1}(*t*) were also obtained for each MCA by averaging the individual impulse responses of all subjects studied. Similarly, the average second-order kernel,*k*
_{2}(*t, t*), was obtained from the mean of all subjects.

The performance of each linear model was assessed by comparing the predicted model response [*V*
_{M}(*t*)] with the real data,*V*(*t*). A standard procedure was adopted for all models using the convolution operation to calculate the model response from the impulse response and the input signal, P(*t*), as
Equation 1where*i*(*m*) can be any of the impulse responses*i*
_{Z}(*t*),*i*
_{T}(*t*),*i*
_{F}(*t*), or*k*
_{1}(*t*), and *L* is the number of lags (see
). Model-estimated velocities were also obtained for the quadratic or second-order kernels by adding the corresponding two-dimensional convolution as described in the
. The same convolution integrals were used to predict the model step responses by substituting P(*t*) with a step function in*Eq. 1
*.

The performance of the different models studied was assessed in three distinct conditions. First, a training set of data with 160-s duration was used to estimate the impulse response of the different models, and these were used to test the predicted velocity response (*Eq. 1
*) and compare it with the observed velocity in the same segment of data. The second comparison used the same impulse response but a different test set of data to obtain the velocity predictions. The test set was chosen as the 160-s segment of data following the training segment. The third analysis attempted to investigate the extent to which baseline-derived models can be used to represent the pressure-velocity relationship recorded after the sudden deflation of thigh cuffs. The large pressure and velocity changes produced by the thigh cuff test present a serious limitation to model identification because of the nonstationarity of the data and the limited duration (∼30 s) of the characteristic velocity response. As suggested by Zhang et al. (33), we used the impulse responses estimated during baseline with the normalized pressure drop recorded during the thigh cuff test to estimate the corresponding velocity response, using *Eq.1
*, and compared the model prediction with the velocity data as described in *Statistical analysis*.

#### Statistical analysis.

Model-predicted velocity responses (*Eq.1
*) were compared with velocity data by means of Pearson’s correlation coefficient (*r*) and the mean square error (MSE). The latter represents a more rigorous test because it takes into account the amplitude similarity between model and data in addition to the temporal pattern that is reflected by the correlation coefficient. The whiteness of model residuals was assessed with the Durbin-Watson test. For each model, means ± SD of*r* and MSE were calculated for all subjects for the right and left MCA responses. Results for the two MCAs were grouped if no significant differences were observed between the right and left sides. Comparisons between models were performed for each subject by testing the change in*r* with the*z* test and the MSE with the*F* test. The sign test was then used to assess the significance of the differences between models for all subjects. A value of *P* < 0.05 was considered significant.

## RESULTS

We studied 47 subjects (27 male) with mean ± SD age of 66.9 ± 9.2 yr (range 44–80 yr), body mass index of 26.5 ± 3.4 kg/m^{2}, heart rate of 75 ± 9 beats/min, clinical MABP of 101.6 ± 11.5 mmHg, systolic BP of 138.3 ± 15.2 mmHg, and diastolic BP of 83.3 ± 10.6 mmHg. Only 33 subjects had thigh cuff responses that were acceptable according to the criteria described previously. There were no major changes in transcutaneous CO_{2} during measurements or between the baseline recordings and the thigh cuff test.

The transfer function analysis results are presented in Fig.1, and the corresponding impulse response functions, estimated using the FFT method, are shown in Fig.2
*A* as mean values for the population studied. At frequencies <0.1 Hz, the mean coherence function was relatively low, suggesting uncoupling between changes in CBFV and ABP. As the frequency increased to ≥0.1 Hz, the squared coherence reached a plateau around 0.6 for both right and left channels. The mean gain, or amplitude-frequency response (Fig.1
*B*), also increased with frequency, in this case rising continuously up to a maximum at 0.2 Hz, leveling off until reaching ∼0.4 Hz, and then rising again to peak at 0.5 Hz. The mean phase-frequency response (Fig.1
*C*) presented an inverse trend; it was more positive below 0.1 Hz and then decreases more or less continuously until crossing zero at ∼0.35 Hz for both channels.

Figure 2
*B* presents the group-average Tiecks impulse response, which was almost identical for both channels, with the smallest SE for all linear models studied. The corresponding mean ± SD of the ARI for the combined right- and left-side responses was 4.5 ± 2.9.

The use of the Wiener-Laguerre method to treat the ABP-CBFV relationship as linear led to the impulse response functions*k*
_{1}(*t*) depicted in Fig. 2
*C*, also showing excellent agreement between the two different CBFV channels. The nonlinear analysis, with the formulation described in the
, produced a first-order term with a temporal pattern that resembles*k*
_{1}(*t*) in Fig. 2
*C* and the mean second-order kernel*k*
_{2}(*t, t*), shown in Fig. 3. The differences between the linear and quadratic implementations of the Wiener-Laguerre method are exemplified by a representative segment of data and the corresponding model-predicted responses in Fig.4. Figure 4 also shows the phase lead of CBFV in relation to ABP during spontaneous oscillations in pressure.

For the training set of data, the superiority of the nonlinear representation, exemplified in Fig. 4, has been confirmed for the whole group of subjects as indicated by the results in Table1. As expected, the zero-order model gave the worst results for both *r* and MSE. The FFT and Tiecks methods yielded similar results, but for the linear assumption the best results were obtained with the linear Wiener-Laguerre approach. When applied to the change in correlation and MSE, the sign test indicated that*1*) the nonlinear method was superior to the linear Wiener-Laguerre method (*P* < 10^{−6});*2*) the linear Wiener-Laguerre method was superior to both the FFT and Tiecks methods (*P* < 10^{−6});*3*) the FFT and Tiecks methods were not significantly different from each other; and*4*) both the FFT and Tiecks methods were superior to the zero-order model (*P* < 10^{−6}). Durbin-Watson tests were inferior to the critical limit (*d*= 1.7, where *d* is the Durbin-Watson statistic), indicating that residuals were serially correlated in all models (Table 1). However, for the nonlinear Wiener-Laguerre model there was a very significant (*P* < 10^{−8}) change of this statistic, showing an improvement in the direction of less correlated residuals.

Despite the fact that the nonlinear Wiener-Laguerre model was much superior to all linear versions (Table 1), the estimated step responses obtained with the quadratic and linear implementations of the Wiener-Laguerre method show only a small difference in the temporal pattern of the CBFV response to a step change in ABP, as indicated in Fig. 5.

Table 2 presents the model performance results when impulse responses estimated from the training set were used to predict the velocity responses of the test set of data. The worst model in this case was the quadratic Wiener-Laguerre model (*P* < 0.0002), followed by the zero-order model (*P* < 10^{−6}). The Tiecks, FFT, and linear Wiener-Laguerre models had comparable performance, with nonsignificant differences by the sign test between these three models, except for the correlation coefficient, which was superior for the linear Wiener-Laguerre model in relation to the FFT model (*P* < 0.05). The Durbin-Watson test indicated that residuals were serially correlated in every case, although the quadratic Wiener-Laguerre model provided significantly higher values in relation to the other models examined.

The results of the estimation of CBFV responses during the thigh cuff test are given in Table 3. In this case the zero-order model performed even worse than in the baseline situations, but neither the mean results for the FFT method nor those for the Tiecks method are very different from those in Table 1, although the SD is larger. The linear Wiener-Laguerre method was only slightly superior to the other linear methods, but the nonlinear model was clearly inadequate, performing even worse than the zero-order model with respect to MSE (Table 3). The sign test for these data indicated that*1*) the nonlinear model was significantly worse than the linear Wiener-Laguerre representation (*P* = 0.0003);*2*) the linear Wiener-Laguerre model was not significantly different from the Tiecks method;*3*) both the Tiecks and linear Wiener-Laguerre methods were superior to the FFT method (*P* = 0.008,*P* = 0.041); and*4*) the FFT method was superior to the zero-order model (*P* = 0.003). For all models the Durbin-Watson statistic indicated that residuals were serially correlated.

## DISCUSSION

Analyses of dynamic cerebral autoregulation based on spontaneous fluctuations of ABP have often adopted transfer function analysis to quantify the relationship between beat-to-beat changes in MCBFV and MABP (5, 9, 10, 14, 17, 23, 33). Several investigators have reported coherence-, gain-, and phase-response curves that allow cerebral autoregulation to be interpreted as a frequency-dependent mechanism (5,9, 14, 17, 22, 23, 33). The common characteristics of these findings are a very significant phase lead of MCBFV in relation to ABP and a coherence function that is relatively high for frequencies >0.2 Hz [say, γ^{2}(*f*) > 0.6] but drops to values of ≤0.50 for frequencies <0.1 Hz (14, 23, 33). This reduction in coherence has been interpreted as the result of active autoregulation operating in this frequency range, whereby changes in ABP are buffered and the corresponding fluctuations in MCBFV are attenuated by feedback adjustments in cerebrovascular resistance. As the frequency of the input ABP changes increases, the autoregulatory mechanism cannot respond fast enough and both the coherence and gain, or amplitude-frequency response, tend to increase. These previous findings have been confirmed by the present study, incorporating data from both MCA arteries of a large number of healthy subjects. We have obtained frequency-domain functions of mean coherence, gain, and phase that are very similar to those presented by Zhang et al. (33). The relatively high values of squared coherence (γ^{2}) obtained by ourselves and others (5, 14, 23, 33) indicate that a large fraction of the beat-to-beat variability of Doppler velocity measurements in the MCA can be explained by fluctuations in MABP. This is an important result because of the ongoing discussion about the possibility of diameter changes of the MCA affecting the flow-velocity relationship. Although Newell et al. (18) have provided evidence that the relationship between flow and velocity remains relatively constant during thigh cuff tests, we cannot rule out the possibility of small spontaneous changes in diameter, leading to CBFV fluctuations of amplitude comparable to what we have observed during baseline recordings. The high values of coherence, however, indicate that if such fluctuations in diameter occur, they have only a minor influence and cannot obscure the relationship between MABP and CBFV. In fact, it is possible to assume that these hypothetical changes in diameter are even less influential because other factors, such as nonlinearity and noise, also have a bearing on the variance unaccounted for (i.e., 1 − γ^{2}). When our results on transfer function analysis are compared with those reported by other investigators, two aspects should be noted. First, we studied an older population, ranging from 44 to 80 yr of age (mean 66.9 yr), but without apparent qualitative differences in relation to studies involving younger subjects. Second, we used dual-channel Doppler recordings, obtaining an excellent agreement between the right- and left-side measurements, not only for the transfer function analysis (Fig. 1) but also for the other modeling approaches investigated (Fig. 2, Table 1). These results are very encouraging regarding the confidence that can be placed on the Doppler technique and the reproducibility of the analytic methods.

Within the linear assumption for the MABP-MCBFV relationship at rest, we have also studied the performance of the model of Tiecks et al. (29) and the Wiener-Laguerre representation proposed by Marmarelis et al. (6, 7, 15, 16). Contrary to our initial expectations, the Tiecks model gave results comparable to those of the transfer function approach, although it was originally proposed and validated for use with the thigh cuff test (29). These results should not be interpreted, however, as an indication that the Tiecks model is the best mathematical representation of the pressure-velocity relationship in the presence of an active autoregulation (31). As an empirical mathematical model, the formulation proposed by Tiecks et al. (29) was aimed solely at grading dynamic autoregulation for step changes in ABP and, consequently, does not provide any insights into the underlying physiology of cerebral autoregulation. On the other hand, the finding that the Tiecks model can provide a reasonable representation of dynamic autoregulation at rest has some relevant implications for the assessment of autoregulation for clinical purposes. Panerai et al. (24) have recently demonstrated that the Tiecks model used in conjunction with baseline measurements can detect the occurrence of carotid artery disease. Similar results could possibly be obtained with the transfer function method, but the Tiecks model has the potential advantage of better performance for short segments of data. As a drawback, this model cannot reflect the same temporal evolution of the negative phase of the impulse responses given in Fig. 2 as those of the transfer function and Wiener-Laguerre impulse responses.

For the linear methods studied, the best results were obtained with the linear Wiener-Laguerre representation (Table 1). This approach has been applied to study renal autoregulation in rats, but we are not aware of other applications involving the human cerebral autoregulation. As described in the
, the method requires the choice of three basic parameters, namely, the lag duration*L*, the number of Laguerre functions (*J*), and the value of α for the Laguerre functions. We have obtained satisfactory results with lags in the range of 16–32 s and values of α in the range of 0.05–0.2. We have adopted *J* = 10 but have not explored the sensitivity of the method to this parameter. Impulse responses obtained with the transfer function and Tiecks methods indicate that the response is fairly completed for lag times of ≤8 s. This feature of dynamic autoregulation makes the Wiener-Laguerre method a convenient approach because of the limited time lag and, hence, computational burden required (15). Despite the apparent practical advantage of the Tiecks model in providing a simple index (ARI) for grading autoregulation, it should be possible to obtain a similar measure for the linear Wiener-Laguerre approach, as indicated by the step response depicted in Fig. 5. Similarly to the Tiecks model, the linear Wiener-Laguerre can also provide robust estimates for much shorter segments of data than we have adopted in our study.

Application of impulse responses derived with the three linear models in the training set to other two segments of data for each subject, the test set and the thigh cuff response, respectively, have shown similar performances for the three methods and also considerable robustness of the linear approach. Although extrapolation of dynamic responses to different segments of data can provide an indication of model robustness, in the case of cerebral autoregulation this test needs to be interpreted with caution because there is no evidence that the dynamic relationship between MABP and CBFV should be immutable. The opposite is more likely, in fact, because this relationship can be influenced by other variables, such as and mental activation (1, 25), that may show small oscillations with time even during baseline recordings. The superiority of the linear Wiener-Laguerre model for the training set, which disappears for the other two segments of data, suggests that this method may be more sensitive to the time-varying influences of other variables that can affect the regulation of CBF.

We and others (5, 22, 23, 33) have used the small-signal argument to justify the utilization of linear models for baseline measurements. The finding that all three linear models produced satisfactory results when baseline impulse responses were applied to thigh cuff recordings suggests that these impulse responses can be regarded as a surrogate for situations with much larger changes in ABP such as the thigh cuff test, tilting, and Valsalva or lower body negative-pressure maneuvers. Zhang et al. (33) have obtained similar results showing that the linear impulse response derived for baseline measurements can predict the MCBFV for the thigh cuff test reasonably well, quoting a mean difference of 0.22 ± 1.45 cm/s between the model prediction and the data. However, according to their Fig. 5, it is likely that the error quoted includes the balance of positive and negative deviations, and therefore it is not comparable to our results using the mean quadratic error (Table 3). Interestingly enough, the Tiecks model presented a performance comparable to that of the linear Wiener-Laguerre representation and was significantly, albeit only slightly, superior to the transfer function approach (Table 3). Despite the fact that the thigh cuff test can provide an almost immediate visualization of the efficiency of the autoregulatory mechanism, our results indicate that similar information can be derived from measurements at rest. The feasibility of substituting the thigh cuff analysis with the modeling of baseline recordings has been recently demonstrated (24), but the present set of results provides more general and encouraging evidence that should facilitate the diffusion of cerebral autoregulation assessment as a clinical tool.

When the linear Wiener-Laguerre representation was expanded to include a quadratic term, excellent results were obtained in the training conditions (Table 1), but extension of the corresponding kernels to the test set (Table 2) and the thigh cuff data (Tables 3) led to disappointing results. It could be argued that the superior results obtained with the nonlinear model for the training set simply reflect the fact that the additional quadratic term can more easily fit the noise rather than reflect a true characteristic of the autoregulatory mechanism. Figure 4 shows that this is not the case and is representative of observations in the majority of subjects. The predicted velocity response for the nonlinear model does seem to represent important features of the MCBFV time series rather than simply following the noise. This conclusion is reinforced by the analysis of residuals as expressed by the Durbin-Watson test. Although none of the models provided uncorrelated residuals (Table 1), the nonlinear Wiener-Laguerre model led to a very significant change in the Durbin-Watson statistic, thus indicating that the quadratic term helps to explain nonrandom features of the CBFV time series. Nevertheless, the nonwhiteness of residuals for all the models studied is of concern and may be due to the fact that the output noise is not itself random because it may be the result of low-frequency oscillations induced by respiration, intracranial pressure, or small changes in MCA diameter, as mentioned previously. The two-dimensional mean kernel for the quadratic term*k*
_{2}(*t, t*), represented in Fig. 3, shows relatively large secondary peaks at lag values of several seconds. This feature is not observed in any of the linear impulse responses in Fig. 2. The presence of these peaks suggests that the second-order kernel may account for some delayed effects of pressure on MCBFV transients. This behavior is exemplified in Fig. 5, which shows that the main difference between the linear and nonlinear estimated responses occurs several seconds after the sudden ABP transition. Although the difference between the linear and nonlinear predicted responses is relatively small, this result cannot be generalized to other input functions because of the presence of the quadratic term.

Unfortunately, the finding that the quadratic term produces significantly better models does not immediately shed much light on the mechanism of autoregulation. Moreover, it does not make sense to attempt to find a physiological interpretation for the second-order term. Whatever the nature of the nonlinearities involved (31), the quadratic term is only the first of their equivalent series expansion (see ), and that it is how it should be interpreted. By applying the nonlinear Wiener-Laguerre method to study the renal autoregulation of rats, Chon et al. (6) have proposed that the third-order kernel is more relevant than the second-order term to represent the nonlinear behavior of the system. The same authors previously analyzed the second-order component in the frequency domain, using the two-dimensional Fourier transform, extracting information about possible interactions between the myogenic and tubuloglomerular feedback mechanisms for control of kidney blood flow (7). Similar investigations could be performed on the cerebral autoregulation, for example, to explore the possibility of interaction between the myogenic and metabolic control pathways (25, 31). One limitation for human studies, however, is that Chon et al. (6, 7) have imposed randomlike changes in ABP to increase whiteness and improve the convergence of estimation for higher order kernels. With spontaneous fluctuations we have no control of the nature and amplitude of the ABP transients, but we have been able to fit the second-order component using 160-s segments of data. However, for shorter segments of data the Wiener-Laguerre method failed because of poor conditioning of the single-value decomposition solution for the second-order kernel (see ). It is an open question as to whether it would be possible to obtain the third-order component by using baseline data of similar duration in humans.

The most appropriate explanation for the poor results of the nonlinear Wiener-Laguerre model, when extended to other segments of data, lies with the temporal pattern of MABP fluctuations. Although the Wiener method requires strict randomness of the input variable, Marmarelis (15, 16) claims that this condition can be relaxed in the case of the Laguerre transformation. Nevertheless, Chon et al. (6, 7) have imposed disturbances in MABP to increase its randomness, and we are not aware that a separate test set of data has been used to validate their results. In our study, only spontaneous fluctuations in MABP were used, and decimation to a sampling rate of 1 s was adopted to increase the randomness of the MABP signal. Because of the larger number of coefficients in the quadratic kernel compared with that in the linear Wiener-Laguerre component (see ), we hypothesize that in the absence of strict randomness the nonlinear Wiener-Laguerre model becomes too sensitive to the particular MABP temporal pattern of the training set. Consequently, when kernels are applied to different sets of data, nonrandom changes in MABP, such as the thigh cuff response, produce dominant quadratic components with aberrant responses, as shown by the results in Tables 2 and 3. Further work is necessary to determine whether increases in randomness of MABP will lead to parallel increases in robustness of the nonlinear Wiener-Laguerre model. With spontaneous fluctuations in ABP, it may be possible to increase randomness by performing much longer baseline recordings (≥20 min) than we have adopted in the present study.

Notwithstanding these limitations, the very significant improvement in model performance obtained with the nonlinear Wiener-Laguerre method in the training segment of data deserves some attention because it may have important implications for clinical investigations and future research into dynamic autoregulation. First, investigators need to be aware of the general limitation of linear models when searching for mathematical models of cerebral dynamic autoregulation. This might be the case of attempts to improve empirical models (29) or obtain models that can provide greater insight into the physiology of the cerebral circulation. Both the FFT and linear Wiener-Laguerre techniques (Fig.2, *A* and*C*) can be regarded as general linear models whose flexibility of coefficients (or impulse response) allows the representation of relatively complex mathematical models formed by a large set of linear differential equations. The differences in performance between the FFT and linear Wiener-Laguerre approaches and the Tiecks model (Table 1) provide an indication of how much can be gained by sticking to the linear formulation. On the other hand, nonlinear models, as proposed by Ursino and Lodi (31), stand a much better chance of explaining the dynamic relationship between MABP and CBFV beat-to-beat changes. Second, the clinical relevance of nonlinear models can only be established by future research. Linear models are much simpler to use in routine applications and have been shown to possess diagnostic value in certain conditions (23, 24, 29). What remains to be established, however, is whether nonlinear models can lead to significant improvements in diagnostic accuracy in different clinical conditions.

## Acknowledgments

S. Dawson was funded by The Stroke Association of the United Kingdom. The support of the Engineering and Physical Sciences Research Council (UK) is also gratefully acknowledged.

## Appendix

#### Tiecks model.

For an ABP drop P(*t*), normalized by the baseline values preceding the deflation of thigh cuffs, the relative pressure change dP(*t*) is given by (29)
Equation A1where CCP_{r} is a parameter introduced by Tiecks et al. (29) to represent the critical closing pressure, which in*Eq. EA1
* is assumed to be normalized by the baseline pressure. The relative velocity change predicted by the model is given by
Equation A2where*K* is a parameter reflecting system autoregulatory gain (29), and*x*
_{2}(*t*) is an intermediate variable obtained from the following pair of equations
Equation A3
Equation A4where*f*, *D*, and *T* represent the sampling frequency, damping factor, and time constant parameters, respectively.

*Equations EA1-EA4
* were implemented with the same combination of parameter values proposed by Tiecks et al. (29) in their Table 3. For each segment of data P(*t*), the model generates a corresponding predicted velocity*V*
_{M}(*t*) that can be compared with the original data*V*(*t*). Ten possible combinations of parameter values (*K, D, T*; corresponding to ARI ranging from 0 to 9) are tested, and the one leading to the minimum square error corresponds to the solution adopted. For this combination of parameters, a step function of pressure is entered as P(*t*), and the*i*
_{T}(*t*) impulse response is obtained as the time derivative of*V*
_{M}(*t*), given by *Eq. EA2
*.

#### Transfer function analysis.

From the temporal sequences P(*t*) and*V*(*t*), the frequency-domain transforms*P*(*f*) and*V*(*f*) are computed with an FFT algorithm. The power spectrum of P(*t*) [G_{PP}(*f*)] is calculated as
Equation A5where the expected value of the complex product E[P(*f*)*P(*f*)] was obtained by smoothing the spectra with a triangular window and averaging four segments of data. Similarly, the cross-spectrum G_{PV}(*f*) was computed as
Equation A6The squared coherence function [γ^{2}(*f*)] was estimated by (3)
Equation A7where G_{VV}(*f*) is the power spectrum of*V*(*t*). The squared coherence reflects the fraction of output power that can be linearly related to the input power at each frequency. Similarly to a correlation coefficient, it varies between 0 and 1.

The complex transfer function H(*f*) was given by
Equation A8With amplitude (gain) ‖H(*f*)‖ and phase spectrum Φ(*f*) obtained from the real part H_{R}(*f*) and imaginary part H_{I}(*f*) of the complex transfer function, this becomes
Equation A9
Equation A10Finally, the impulse response*i*
_{F}(*t*) is computed as the inverse FFT transform of the transfer function
Equation A11where Ϝ^{−1} represents the inverse Fourier transform (3).

#### Marmarelis method.

The Volterra series expansion for a nonlinear dynamic system can be written as
Equation A12where P(*t*) is the ABP input and*V*
_{M}(*t*) is the model velocity output;*k*
_{1}(*t*) and*k*
_{2}(*t, t*) are the first- and second-order kernels, respectively; and*L* is the number of lag time intervals used to represent the duration of the kernels. For the linear implementation,*k*
_{2}(*t, t*) and all other higher order terms are assumed to be zero.

As suggested by Wiener, this series can be rewritten using an instrumental set of orthogonal functions
Equation A13where *J* is the number of functions used in the decomposition and*c* values are coefficients. The orthogonal basis functions ϑ_{j }(*t*) are obtained from the convolution integral

Equation A14where*L _{j}
* (

*t*) corresponds to the

*j*th-order Laguerre function given by Equation A15The parameter α controls the amount of damping of the Laguerre function and needs to be carefully selected for different applications. One advantage of using the Laguerre series is the fact that, in practice,

*Eq. EA14*does not need to be computed because the orthogonal set ϑ

_{j}(

*t*) can be obtained with the recursive relation Equation A16which can be initialized with Equation A17where

*T*is the sampling interval.

In summary, the basis functions ϑ_{j }(*t*) are calculated with *Eqs. EA16
* and *
EA17
*. We assume*c*
_{0} =*k*
_{0} = 0, because P(*t*) and*V*(*t*) are zero mean. From *Eq. EA13
*, the set of coefficients*c*
_{1}( *j*) and*c*
_{2}( *j, j*) can be calculated by single-value decomposition (16). Finally, the first- and second-order kernels can then be obtained from
Equation A18
Equation A19With the use of these estimates, it is then possible to obtain model-predicted responses with *Eq.EA12
*.

## Footnotes

Address for reprint requests and other correspondence: R. B. Panerai, Dept. of Medical Physics, Leicester Royal Infirmary, Leicester LE1 5WW, United Kingdom (E-mail: rp9{at}le.ac.uk).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

- Copyright © 1999 the American Physiological Society