Medical Engineering & Physics 36 (2014) 1636–1643

Contents lists available at ScienceDirect

Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy

Rapid pressure-to-flow dynamics of cerebral autoregulation induced by instantaneous changes of arterial CO2 Jia Liu a,b,∗ , David M. Simpson c , Hesam Kouchakpour c , Ronney B. Panerai d , Jie Chen e , Shan Gao e , Pandeng Zhang a,b , Xinyu Wu a,b a Guangdong Provincial Key Laboratory of Robotics and Intelligent System, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, No. 1068 Xueyuan Avenue, Shenzhen University Town, Shenzhen 518055, PR China b Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, Shatin, Hong Kong, PR China c ISVR, University of Southampton, University Rd., Southampton SO17 1BJ, UK d Department of Cardiovascular Sciences and National Institutes for Health Research (NIHR), Biomedical Research Unit, University of Leicester, Leicester Royal Infirmary, Level 1, Sandringham Building, Leicester LE1 5WW, UK e Department of Neurology, Peking Union Medical College Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, No. 1 Shuaifuyuan, Dongcheng District, Beijing 100730, PR China

a r t i c l e

i n f o

Article history: Received 28 November 2013 Received in revised form 12 August 2014 Accepted 7 September 2014 Keywords: Cerebral hemodynamics Multivariate model Adaptive filters Hilbert transform

a b s t r a c t Continuous assessment of CA is desirable in a number of clinical conditions, where cerebral hemodynamics may change within relatively short periods. In this work, we propose a novel method that can improve temporal resolution when assessing the pressure-to-flow dynamics in the presence of rapid changes in arterial CO2 . A time-varying multivariate model is proposed to adaptively suppress the instantaneous effect of CO2 on CBFV by the recursive least square (RLS) method. Autoregulation is then quantified from the phase difference (PD) between arterial blood pressure (ABP) and CBFV by calculating the instantaneous PD between the signals using the Hilbert transform (HT). A Gaussian filter is used prior to HT in order to optimize the temporal and frequency resolution and show the rapid dynamics of cerebral autoregulation. In 13 healthy adult volunteers, rapid changes of arterial CO2 were induced by rebreathing expired air, while simultaneously and continuously recording ABP, CBFV and end-tidal CO2 (ETCO2 ). Both simulation and physiological studies show that the proposed method can reduce the transient distortion of the instantaneous phase dynamics caused by the effect of CO2 and is faster than our previous method in tracking time-varying autoregulation. The normalized mean square error (NMSE) of the predicted CBFV can be reduced significantly by 38.7% and 37.7% (p < 0.001) without and with the Gaussian filter applied, respectively, when compared with the previous univariate model. These findings suggest that the proposed method is suitable to track rapid dynamics of cerebral autoregulation despite the influence of confounding covariates. © 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

1. Introduction Patient specific continuous assessment of CA can help in determining the optimal blood pressure to achieve favorable outcomes in a number of clinical conditions, where cerebral hemodynamics may change within relatively short periods. This may occur for example in patients in a neurosurgical intensive care unit, for

∗ Corresponding author at: Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, No. 1068 Xueyuan Avenue, Shenzhen University Town, Shenzhen 518055, PR China. Tel.: +86 755 86392138; fax: +86 755 86392299. E-mail address: [email protected] (J. Liu). http://dx.doi.org/10.1016/j.medengphy.2014.09.005 1350-4533/© 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

example following subarachnoid hemorrhage or head injury, or during cardiovascular surgery [4,10,35,39]. Univariate linear models of dynamic CA, such as the moving window based autoregulation index, transfer function analysis, and linear regression, using beat-to-beat data of CBFV and ABP (or cerebral perfusion pressure—CPP) have been proposed to track this vital physiological function continuously [6,10,26]. However, the sensitivity and specificity of these methods might be constrained by a number of physiological confounding covariates as well as the limitations of modeling techniques that are applied [6,26,28]. Arterial CO2 , usually measured as the partial pressure of end-tidal CO2 (ETCO2 ), is known as a powerful cerebral vasodilator that can change blood flow greatly within seconds [33]. In

J. Liu et al. / Medical Engineering & Physics 36 (2014) 1636–1643

addition, hyper- or hypocapnia induced by increased or decreased ETCO2 can result in significant changes in the effectiveness of autoregulation [1,3,23], and a number of studies have used modeling to understand the underlying variations in pressure-to-flow dynamics [2,9,14,17,18,20,21,24,31,34,37]. Panerai et al. showed that a multivariate model with ETCO2 as an additional input can improve the prediction of CBFV [24]. Peng et al. [31] also reported that the linearity of the system can be significantly increased when using a multiple input model, especially in the frequencies below 0.1 Hz. In a recent study of phase dynamics, they also found that the low synchronization of phase difference (PD) between ABP and CBFV can be attributed to CO2 , which can be corrected by a CO2 term derived from a multivariate model [32]. There thus appears to be a growing consensus that a multipleinput-and-single-output (MISO) system is more suitable to model the cerebral hemodynamics than a single-input-and-single-output (SISO) system [14,17,18,32]. In addition to the multivariate characteristics, time-varying assessments are also of growing interest [8,14,18,25,29,30]. Dineen et al. proposed a continuous autoregulation index and showed that the response of dynamic autoregulation to a change of ETCO2 is delayed with respect to the change in CBFV [8]. This finding agrees with the results reported by Liu et al. using adaptive filtering techniques to track time-varying dynamics of autoregulation in response to step-wise changes of CO2 [16]. These studies indicate that the pressure-to-flow dynamics might not be changing at the same pace as the dilation of arterioles. However, univariate models were employed by these studies and the tracking speed might be relatively slow, as the updating rate of the continuous estimate is either constrained by the relatively long length of the moving window or the slow forgetting factor for the adaptive filter. We thus aim now to elucidate if the reported delay in autoregulatory response is due to physiological phenomena or reflects the limitations of the signal processing methods. Autoregulation can be quantified by many parameters. Instantaneous phase dynamics has been reported as a standalone autoregulatory parameter in a number of recent papers, using Hilbert transform (HT) and wavelet transform techniques [11,12,15,22,32]. The instantaneous phase estimated by these techniques is intrinsically a continuous parameter, as it is calculated sample by sample from the recordings. This approach therefore readily lends itself to the investigation of the time-varying property of autoregulation. Previous publications have shown that dynamics of autoregulation can be assessed within a relatively narrow band around 0.1 Hz [1,3,7,15,16,22]. The choice of bandwidth (and thus the band-pass filter that needs to be applied to the signals prior to using the Hilbert transforms [5]) requires a compromise between temporal and frequency resolution. In the present study, we thus combine multivariate modeling and instantaneous signal processing methods to shed further light on the continuous evolution of CA with improved temporal resolution in the presence of rapid changes in arterial CO2 . The objective is to provide a method that may track rapid changes of dynamic CA despite the influence of confounding covariates.

2. Methods 2.1. Outline The approach may be summarized as follows: we model the cerebral hemodynamics as a multivariate system. The recursive least square (RLS) method is applied to adaptively remove the contribution of changes in ETCO2 to CBFV. Gaussian filters with minimal spread in the time-frequency domain are used prior to applying the RLS filter to limit fluctuations to the frequency band

1637

where autoregulation is usually considered as most evident. The instantaneous PD is then estimated from the HT of the filtered residual CBFV and ABP. In the following section a simulation is described to validate our method, which was then applied to the recorded data to show the underlying time-varying pressure-toflow dynamics during changes of ETCO2 induced by rebreathing. 2.2. Mathematical methods Cerebral autoregulation is modeled as a multivariate system [24], where ABP (p[n]), and ETCO2 (c[n]) are the inputs and CBFV (v[n]) is the output,

v[n] = vp [n] + vc [n] + e[n] Lpv −1

=





Lc v −1

hpv [i]p[n − i] +

i=0

hcv [i]c[n − j] + e[n],

(1)

j=0

and n denotes sample number at 1 Hz sampling rate. e[n] is the residual CBFV unexplained by the model. vp [n] and vc [n] are parts of CBFV that can be attributed to ABP and ETCO2 , respectively; hpv [i] and hcv [i] are their respective causal FIR filter impulse responses; Lpv and Lcv denote the orders of the filters. Based on a compromise between the known time for a physiological response to occur, the number of free parameters, and the length of data available, these orders are chosen to be equal 10 for both filters. Based on (1), we extended our previous time-varying univariate model identified by a SISO RLS adaptive filter to a time-varying multivariate model identified by a MISO RLS adaptive filter. We rearrange (1) by vectors, where the vector of the filter impulse responses can be written as:



H [n] =

H pv [n]



,

H cv [n]



(2)



t

H pv [n] = hpv [0] . . .hpv Lpv − 1

,

(3)

t

H cv [n] = [hcv [0] . . .hcv [Lcv − 1]] .

(4)

The inputs can be defined as:



x [n] =

xp [n] xc [n]



,



(5)



t

xp [n] = p [n] , p [n − 1] , . . ., p n − Lpv + 1

,

xc [n] = [c [n] , c [n − 1] , . . ., c [n − Lcv + 1]]t .

(6) (7)

According to (1), the error can be expressed as: e [n] = v [n] − H[n]t x [n] .

(8)

We can then follow the algorithm of the RLS adaptive filter to update H[n]. We first estimate the Kalman gain vector, k[n], as: k [n] =

−1 P [n − 1] x [n] 1 + −1 xt [n] P [n − 1] x [n]

(0 ≤  ≤ 1) ,

(9)

where P[n] is the inverse autocorrelation matrix of the input signals. This is updated by: P [n] = −1 P [n − 1] − −1 k [n] x[n]t P [n − 1] .

(10)

The vector of the filter impulse responses, H[n], is then updated by: H [n] = H [n − 1] + en [n] k [n] .

(11)

1638

J. Liu et al. / Medical Engineering & Physics 36 (2014) 1636–1643

To this end, the residual CBFV, vr [n], with vc [n], the component of CBFV attributed to ETCO2 removed, can be estimated sample by sample as,

vr [n] = v [n] − H tcv [n] xc [n] .

(12)

Following HT, the instantaneous PD can then be estimated between the analytic signals of ABP (pa [n])and the residual CBFV (var [n]), which are defined as: pa [n] = p [n] + jp˜ [n] = Ap [n] ejϕp [n] ,

(13)

var

(14)

[n] = vr [n] + j˜vr [n] = Avr

[n] ejϕvr [n] ,

where p˜ [n] and v˜ r [n] are the corresponding HTs. The instantaneous PD can then be calculated as: ϕvr [n] − ϕp [n] = arg

 va

r [n] pa [n]



.

(15)

The right hand side in (15) is preferred to simply subtracting the phases of the signals p and vr , as this reduces discontinuities in results due to phase unwrapping. In order to obtain meaningful results for the phase of the analytic signal and to focus the analysis on a clinically relevant frequency band, band pass filtering is required prior to HT to extract the signal component that is of interest. In the current work, this filter was applied to all the recorded signals prior to RLS system identification, to ensure that the minimization of the mean-squared error is focused on the relevant frequency band. The specification of this band-pass filter requires a compromise between temporal and frequency resolution of the instantaneous phase [5]. Given that the Gaussian function gives minimal spread in both time and frequency domains [13], we defined a Gaussian filter, with an impulse response, GF[n], as follows: GF [n] =



1

2 / ) gf

2gf

sin (2fc n) e(−n

,

(16)

where fc is the center frequency and  gf is the bandwidth parameter. We set the center frequency, fc , at 0.1 Hz for ABP and CBFV, and at 0.05 Hz for ETCO2 , as these are the most commonly used frequencies to show the phase dynamics [1,5,8,15,32,38]. In order to compare the performance with alternative band-pass filters, the temporal and frequency resolution are calculated as  t and  f , respectively [19]:

N/2

t2

=

n=−N/2

N/2

n2 h2 [n]

n=−N/2

M f2 =

m=0

M

h2 [n]

m2 |H|2 [m]

m=0

|H|2 [m]

,

(17)

,

(18)

where n and m are the sample of time and frequency, respectively; N and M stand for the length of the impulse response and frequency response; h and H refer to the impulse and frequency responses, respectively. 2.3. Simulation study In order to show the performance of the proposed method in an ideal condition, a simulation test was designed to compare four different methods of estimating time-varying parameters. Step changes of PD were simulated by generating sinusoidal ABP changes at 0.1 Hz, that is: ps [n] = sin(2 0.1 n), Corresponding fluctuations in CBFV were generated with Tiecks filters [36] in a sequence of autoregulation indexes (ARI) of 7, 3, and, 5 and represented as vsp [n]. The simulated ETCO2 was provided by another sinusoidal wave at 0.05 Hz cs [n] = sin(2 0.05 n) and a sequence of phase changes (changing at the same time as the ARI) with values

of ␲/2, ␲, and ␲/2 were used to simulate, vsc [n] = sin(2 0.05 n + ϕ), ϕ = (/2), , (/2). The resultant simulated CBFV is therefore given by vs [n] = vsp [n] + vsc [n]. We then compared performance of the following techniques: 1. PD estimated by HT with the effect of CO2 removed by a multivariate model: We estimated the residual CBFV, vsr [n], by (14) and then calculated the instantaneous PD between the estimated residual CBFV and the simulated ABP, ps [n], using HT. This is the technique as described in methods above. 2. PD estimated by HT with the effect of CO2 removed by a univariate model: Instead of using the proposed multivariate model, we also tried to extract the effect of CO2 on CBFV using a univariate model, considering ETCO2 as the input and CBFV as the output. A SISO RLS adaptive filter was applied to estimate the simulated CBFV that can be attributed to CO2 , vsc [n]. The residual CBFV was then calculated by subtracting this estimate from the estimated CBFV, vs [n] The instantaneous PD was then calculated from the residual CBFV, vsr [n], and the simulated ABP, ps [n], using HTs. 3. PD estimated by HT without removing the effect of CO2 : We applied HT directly to both the simulated CBFV, vs [n], and simulated ABP, ps [n]. 4. PD estimated by RLS adaptive filter only: Following the approach used in our previous work [16], the PD was also calculated directly from the SISO model with ps [n] and vs [n]as the input and the output, respectively. The RLS algorithm was again used for system identification and hence to obtain sample-by-sample estimates of PD. In order to remove the CO2 effect for the step changes of PD, we chose a fast adaptation with a forgetting factor at 0.95 (below which erratic results were observed), except for the last case, where we compared with our previous RLS method with a forgetting factor of 0.98. 2.4. Experimental methods Data collection was approved by the Local Research Ethics Committee of the Peking Union Medical College Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College and informed consent was signed by each participant. Thirteen healthy individuals (24–52 years old) were consecutively recruited in the Department of Neurology, Peking Union Medical College Hospital (Beijing, China), who were free from known cardiovascular or cerebrovascular diseases or specific risk factors (including hypertension, diabetes mellitus, dyslipidemia, smoking). CBFV in middle cerebral artery (MCA) was measured by transcranial Doppler (TCD, mutidop-X, DWL, Germany) with a 2 MHz transducer fixed to the head with the Marc 500 TCD probe fixation headframe (Spencer Technologies, USA). Continuous beat-to-beat ABP was recorded by a servo-controlled finger plethysmograph (Nexfin-cc, BMEYE, Netherlands) in the right hand positioned at heart level. ETCO2 was measured in mmHg with an infrared capnometer (OEM-Module, Envitec, Germany) during nasal expiration. All the subjects breathed spontaneously in a supine position with their heads slightly elevated. After a resting time of 10 min, the CBFV, ABP, and ETCO2 were recorded for at least 2 min. Rebreathing was then provoked by connecting a half blocked long tube (1 m, 600 mL) to a face mask to increase the inhalation of CO2 for 2 min before a return to normal breathing for another 2 min. All the data were recorded with a sampling rate of 100 Hz and stored on a hard disk for further analysis. Recordings of ABP, CBFV and ETCO2 were filtered by an anti-alias low-pass filter with a cutoff frequency at 0.5 Hz and then down-sampled to 1 Hz. In order to obtain the narrow band signals with minimal spread in both time and frequency domain, all signals were then filtered by the

J. Liu et al. / Medical Engineering & Physics 36 (2014) 1636–1643

2.5. Statistical analysis

80

Phase (degree)

proposed Gaussian band-pass filter with center frequency fc = 0.1 Hz for ABP and CBFV and fc = 0.05 Hz for ETCO2 , respectively. The first and last minute of data were discarded due to possible noise caused by movements of the subjects at the beginning and end of data recording. We then applied the same methods described in the simulation to the recorded data to track the rapid phase dynamics.

1639

HT with CO2 effect removed by a multivariate model HT with CO2 effect removed by a univariate model HT SISO adaptive filter Liu et al. 2010 Ground Truth

60

40

20

• Time-course of variables: In order to show the effect of rebreathing on ETCO2 , CBFV, and the instantaneous phase dynamics across the volunteers, the distributions of the samples at each time point for the 13 subjects were calculated. Mean ± standard deviation (std) was used for ETCO2 and CBFV across 13 subjects, and median ± semi interquartile deviation (SID) was computed for PD in order to avoid possible outliers introduced by the discontinuity from phase wraparound occasionally noted. • Test of the effect of rebreathing on ETCO2 and autoregulation estimates: The blocks of steady states were chosen and average ETCO2 and PD calculated. One-way ANOVA was applied to test if rebreathing can significantly change ETCO2 in the volunteers and the Friedman test was used with pairwise post hoc to test if the described methods can differentiate between the dynamic autoregulation at these three steady states. • Comparison of the spectra of the residues from different models: The NMSEs from these methods with and without the proposed Gaussian filter applied were calculated in time domain. Student’s t-tests were used to evaluate if the proposed method can better fit the hemodynamic system output. 3. Results 3.1. Simulation study Fig. 1 shows that, when applying HT directly on ABP and CBFV without removing the effect of CO2 on CBFV, the estimated PD (red line) is markedly contaminated with fluctuations of the simulated ETCO2 at 0.05 Hz. With the effect of CO2 removed by the multivariate model (black line), the noise from ETCO2 is almost completely suppressed after a few seconds. The transient errors can be explained by filter adaptation to step changes in phase shift due to changing levels of ARI. However, with the effect of CO2 removed by a univariate model (blue line), though the effect is also suppressed, the PD is deviates from the ground truth (dashed cyan line), indicating a biased estimate. When applying a univariate adaptive filter to calculate PD (green line) in accordance with previous work [16], the phase difference changes more slowly than with estimates from HT and an evident contamination by the CO2 component is observed. The biased estimate (blue line) in Fig. 1 appears to arise from the relatively fast adaptation of the RLS filter, with changing gains over a single cycle of the lower-frequency ETCO2 signal. This leads to a distortion in the predicted lower frequency component, i.e. it is no longer a pure sine-wave but includes additional frequencies,

0

0

100

200

300

400

500

600

Time (s) Fig. 1. Step changes of phase dynamics estimated from simulated data. (1) The black line is the PD estimated by HT with CO2 effect removed by a multivariate model; (2) the blue line is estimated by HT with CO2 effect removed by a univariate model; (3) the red line is estimated by HT only; (4) the green line is the time-varying phase estimated by our previous method [15], using an adaptive filter and no HT; (5) the dashed cyan line is the ground truth calculated from Tiecks’ model (ARI = 7, 3, and 5). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

6 Simulated Vs [n] c s Distorted Esitmated V [n] c by a univariate model

4

%/%

Continuous samples of all recordings and PD estimated from the models were aligned by the time of the start and end of rebreathing, respectively. The samples were divided into 5 blocks (one minute of continuous samples for each block), which consist of three steady states (Normocapnia—before rebreathing, Hypercapnia, and Normocapnia—after rebreathing) and two transients (CO2 onset and offset from the start and end of rebreathing for 1 min, respectively). The following steps of analysis were carried out and statistical tests were applied with a significance level of 0.05.

2

0

−2

0

100

200

300

400

500

600

Time (s) s

Fig. 2. The simulated (black line) and the predicted vˆ c [n] (gray line) signals. When a relatively fast adaptation of the RLS filter is chosen (forgetting factor = 0.95), the adaptive filter gains change over a single cycle of the input signal and thus distort the predicted signal, generating higher frequency components that interfere with s the simulated vˆ c [n] component.

as illustrated in Fig. 2. Subtracting this predicted CBFV signal from the original left-shifts the oscillation at the higher frequency, and biases the estimated PD, as observed in Fig. 1. Slower adaptation of the filter removes this effect—but of course also leads to slower adaptation to changing autoregulation. 3.2. Physiological study The effect of rebreathing on ETCO2 is illustrated in Fig. 3. The left half of the plot is the mean ± std of recorded ETCO2 aligned by the start of rebreathing, whereas the right half is aligned by the end of rebreathing, given that the duration of rebreathing was not identical in all subjects. Table 1 shows that ETCO2 was significantly elevated (by approximately 16 mmHg) after the start of rebreathing and drops back to the steady state when rebreathing stops. Similarly, in Fig. 4, the mean ± std of recorded CBFV to the changes of ETCO2 are plotted. It is evident that CBFV responds to the changes of ETCO2 very quickly, but as expected it takes longer for CBFV to increase than to drop back to baseline at the end of the manoeuvre [8,16]. Table 1 also shows that all methods can detect significant change of PD induced by rebreathing.

39.8 ± 33.9

55 Aligned by CO onset 2

ETCO (mmHg)

45

2

40.6 ± 25.0

M3

39.7 ± 18.9



M2

20 0

38.4 ± 27.0

M1 M1

50

100

150 Time (s)

200

250

300

Fig. 3. The mean (thick solid line) ± std (thin dashed lines) ETCO2 recorded from the 13 subjects during the start and end of rebreathing. The ETCO2 is elevated by approximately 16 mmHg (see Table 1) after the start of rebreathing and drops back to the steady state at the end of the manoeuvre.

8.4 ± 23.0 15.6 ± 17.0 14.4 ± 24.0 54.1 ± 12.1

Fig. 4. The mean (thick solid line) ± std (thin dashed lines) of CBFV change recorded from the 13 subjects during the start and end of rebreathing. As expected, the decrease of CBFV after the CO2 offset is steeper than the increase of CBFV induced by CO2 onset.

51.1 ± 12.3 50.9 ± 19.5 57.6 ± 23.5 Median ± SID of PD (degree)

M1-4 stands for Methods 1–4 as defined in Fig. 5. * p < 0.05 N1 vs H. † p < 0.05 H vs N2.  p < 0.05 N1 vs N2.

M1 Model

N1 Mean ± std of ETCO2 (mmHg)

30.7 ± 4.1*

*

*

M3 M2

*

M4

*

16.8 ± 27.0



M2



M3



M4

35

25



31.3 ± 4.8 46.2 ± 5.3†

40

30



N2 H

Aligned by CO2 offset

50



M4

J. Liu et al. / Medical Engineering & Physics 36 (2014) 1636–1643 Table 1 Test of the efficacy of rebreathing for inducing hypercapnia and the models for differentiating autoregulation at three ETCO2 levels. One-way ANOVA was applied to test if ETCO2 levels are different at 3 steady states and Friedman test was used with pairwise post hoc to test if the models can detect the difference of dynamic autoregulation at the steady states, where N1: normocapnia (before CO2 onset), H: hypercapnia, and, N2: normocapnia (after CO2 offset) denote three steady states, respectively.

1640

Fig. 5 shows the phase dynamics estimated from the recorded signals. The colored lines are PD estimated by different models. It is evident that without removing the effect of ETCO2 , the instantaneous PD (red lines) overshoots markedly when subjects started and stopped rebreathing. The overshoot of the phase dynamics during the step-wise CO2 changes is suppressed when employing the proposed multivariate model (black line) to remove the partial CBFV introduced by ETCO2 . The univariate model to remove the effect of CO2 is not shown here, given the results of our simulation tests. There is no overshoot observed (green line) when a univariate adaptive filter with a relatively slower tracking speed (forgetting factor = 0.98) was applied [16]. However, when we reduced the forgetting factor to 0.95, the overshoot becomes evident. Fig. 6 shows the averaged spectrum of the measured CBFV compared with the spectra of the residue estimated from different models (with and without the Gaussian filter applied). The NMSE of these methods are listed in Table 2. Using the timevarying approach, shown in Fig. 6(a) and (c), the multivariate model produces much smaller errors than the univariate model, especially below 0.1 Hz. Table 2 indicates that reductions of 38.7% and 37.7% of NMSE were achieved by the multivariate models without and with the Gaussian filter applied, respectively, which is highly significant when compared with the univariate cases (p < 0.001).

J. Liu et al. / Medical Engineering & Physics 36 (2014) 1636–1643

120 100

140

80

120

60

100

40

80

20

60

0

40

−20

20

−40 0

50

100

150 Time (s)

200

250

SID of Phase (degree)

Median of Phase (degree)

Table 2 The NMSE estimated from 4 different models.

160 Aligned by CO2 offset

Aligned by CO2 onset

0 300

4. Discussion 4.1. Main findings The simulation study validated the proposed approach for detecting instantaneous changes of PD with a reduced effect of CO2 from neighbouring frequency bands (0.05 Hz). This effect can also be reduced by a univariate adaptive filter, but biased estimates were

Time-varying univariate model (Liu et al., 2010 [16]) Time-varying multivariate model Time-invariant univariate model Time-invariant multivariate model

58.6*

55.6*

19.9**

17.9**

74.3

70.0

40

CBFV

30

Residue from a time−varying univariate model

2

Residue from a time−varying multivariate model

20 10

10 0.1

0.2 0.3 Frequency (Hz)

0.4

0 0

0.5

c) 20

67.3†

Residue from a time−varying multivariate model

0.1

0.2 0.3 Frequency (Hz)

0.4

0.5

d) 20 CBFV

15 Residue from a time−varying univariate model

10

Residue from a time−varying multivariate model

5

0.1

0.2 0.3 Frequency (Hz)

0.4

0.5

(cm s−1)2Hz−1

CBFV −1 2 −1 (cm s ) Hz

71.0

50

Residue from a time−varying univariate model

20

0 0



observed (Fig. 1). When both ABP and CO2 drive CBFV, but only ETCO2 is used in the univariate model, the adaptive filter appears to falsely ascribe some of the effect of ABP to ETCO2 leading to errors in correcting CBFV for the effect of CO2 and consequently errors in phase difference. This problem was largely overcome by the multivariate filter. When comparing with our previous method of a univariate adaptive filter, it is evident that the tracking speed of the current methods is faster as the delay of our previous method is observed in both simulated and real data (Figs. 1 and 5). This finding suggests that changes of autoregulation, as expressed by PD, that were slower than CBFV in response to changes in CO2 in previous studies were due to the slower tracking speed of the models adopted. It is known that hypercapnia can increase CBFV and reduce the effectiveness of CA [1,3,23]. However, the underlying dynamics of how CA changes during the transient changes of CO2 has not been fully described. Without removing the effect of CO2 on CBFV, the PD estimated by HT overshoots markedly during the transient changes of CO2 (Fig. 5). When the proposed method of removing the effect

(cm s−1) Hz−1

−1 2

(cm s−1) Hz

NMSE with the Gaussian filter (%)

CBFV

30

0 0

NMSE (%)

b) 60

a) 60

40

Model

* p < 0.05 t-test for comparing univariate adaptive filter with univariate Wiener filter. ** p < 0.01 t-test for comparing multivariate adaptive filter with univariate adaptive filter and uni-/multivariate Wiener filter. † p < 0.05 t-test for comparing multivariate Wiener filter with univariate Wiener filter.

Fig. 5. The colored lines are the median (solid lines, left y-axis) and SID (dashed lines, right y-axis) of the continuous PD estimated from 4 different methods. Method 1 (red): HT applied directly to CBFV and ABP. Method 2 (black): HT with CO2 effects removed by an MISO adaptive filter. Method 3 (green): an SISO adaptive filter with forgetting factor = 0.98. Method 4 (blue): an SISO adaptive filter with forgetting factor = 0.95. Without removing the effect of ETCO2 , marked overshoot of the PD (red line) is observed. The proposed method has successfully suppressed the overshoot (black line). Compared without previous method (green line), the overshoot becomes evident (blue line) when adaptation is increased with a reduced forgetting factor of 0.95. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

50

1641

15 Residue from a time−varying univariate model 10

Residue from a time−varying multivariate model

5

0 0

0.1

0.2 0.3 Frequency (Hz)

0.4

0.5

Fig. 6. Spectra of the measured CBFV averaged across the 13 subjects and the averaged spectrum of the error (residue) estimated from different models. It is evident that the error in (a) from the time-varying multivariate model in is markedly lower than from other methods, especially below 0.1 Hz. In order to show the improvement that the multivariate and time-varying models can bring in terms of the prediction of CBFV, we also applied univariate and multivariate Wiener filters (time invariant) in (b), to make the results comparable. The results in (c) and (d) are the cases with the Gaussian filter applied. The corresponding statistical analysis is given in Table 2.

1642

J. Liu et al. / Medical Engineering & Physics 36 (2014) 1636–1643

of CO2 was applied to recorded data, it appears that only a fraction of the phase change is attributed to the pressure-to-flow dynamics, suggesting that the influence of the transient changes of CO2 on CBFV is prominent and confounding. Therefore, when assessing CA during the transient changes of CO2 , the effect of CO2 should be removed. Though there are no overshoots when using our previous method (univariate adaptive filter), it becomes evident if the adaptation is increased by reducing the forgetting factor to 0.95, suggesting that the absence of the overshoot using our previous method is due to a relatively slow adaptation or the absent inclusion of CO2 . Previous studies have shown that models with an additional input of CO2 can improve the prediction of CBFV and coherence [24,31], as CBFV responds to changes in PaCO2 . We found that the time-varying method combined with multivariate model can significantly reduce (by more than 50%) the prediction error (Table 2) below 0.1 Hz, indicating that the time-varying characteristic of the multiple-input system is prominent at low and mid-frequency ranges. This result agrees with the evident nonstationary behaviour of cerebral hemodynamics in the same frequency range shown by Mitsis et al. who used continuous estimates from moving windows over time [20]. Both works found that including ETCO2 in an adaptive model can greatly improve data fitting below 0.1 Hz. All methods have a number of parameters that can be set, including forgetting factor and model order. These can significantly affect the results, and can increase or reduce the benefit of one method over others. During extensive experiments, we found that the improvement in the combination of time-varying and multivariate RLS approach was consistent when compared to the other methods. 4.2. State of the art in cerebral autoregulation modelling Hybrid models combing multivariate and time-varying features have been proposed in recent studies [14,18], also allowing for tracking time-varying autoregulation caused by the fluctuations of arterial CO2 . The authors observed time-varying autoregulation, however, with relatively slow changes. This is different from what we speculated. These models use either adaptive filters or moving windows to update the time-varying autoregulatory parameters. Hence, one possible explanation is that the tracking speed of these models is constrained by the forgetting factors or the length of the moving windows. Comparing with these models, the advantage of using instantaneous phase dynamics estimated from the Hilbert or wavelet transforms is that it is a direct estimation between two signals without the constraint of pre-defined models and it is intrinsically a continuous parameter [11,12,15,22,32], which allows for tracking rapid phase dynamics. To the best of our knowledge, this is the first attempt of using instantaneous phase to track timevarying dynamics of cerebral autoregulation within a multivariate context, removing the powerful and potentially confounding effect of changing CO2 levels. The use of the proposed band-pass filter ensures that the PD is calculated from two signals centered at the same frequency and within the same frequency range. In related previous work [22] using empirical mode decomposition (EMD) this cannot be guaranteed as the Hilbert-Huang transform applied to different signals may not lead to identical bands, with the consequent difficulties in interpreting the physical meaning of PD. The wavelet method provides an alternative set of filters [15,32], but when only one ‘scale’ is required, wavelet analysis is unnecessarily complex (most of the transform data is discarded) and restrictive (usually doubling bandwidth between scales and thus not allowing independent choice of centre-frequency and bandwidth). While the proposed method has clear advantages, its performance should be rigorously compared (and on independent datasets) against alternatives before final recommendations can be made.

4.3. Clinical perspectives Clinical applications of cerebral autoregulation indices have been limited by poor reproducibility and uncertainty about which methods provide greater sensitivity and specificity for diagnostic and prognosis in patients with different cerebrovascular conditions [6,11,28,30,35]. Our results suggest that taking into account the contribution and often confounding effects of fluctuations in PaCO2 could lead to more reliable indices with the potential to improve assessment of dynamic CA in clinical settings. Further work is needed to test this hypothesis, but to do so, it is essential that future studies always include continuous recordings of PaCO2 , either as capnographic recordings of EtCO2 or other methods such as mass-spectroscopy [33]. Although our approach, based on a multivariate RLS model has shown superior results when compared to previous implementations, alternative mathemarical models could also be considered to remove the separate influences of PaCO2 on CBFV and dynamic CA [14]. Moreover, multivariate modelling can readily be extended with more additional inputs, if these can be measured continuously. Obvious candidates are intracranial pressure (ICP) and arterial O2 , or perhaps critical closing pressure or resistance-area product [27] The relevance of these additional co-variates would be pathologydependent, with ICP being more relevant in severe head injury and PaO2 in patients at the risk of hypoxia. 4.4. Limitations In this study, all four methods could detect significant changes of PD during step-wise changes of ETCO2 . However, our recordings were not suitable for detecting the changes of PD in each individual as we only had one recording per subject. Repeated tests for each subject would therefore be needed to test the sensitivity of these methods in the future. A face mask was worn by the healthy volunteers to measure ETCO2 . This experimental setting might not be well tolerated by many patients. Alternative arrangements such as nasal sampling should be considered and evaluated to check its feasibility in clinical studies. 5. Conclusion In summary, we showed that a multivariate RLS approach can be more sensitive and accurate in tracking time-varying cerebral autoregulation than a univariate adaptive filter to take into account the simultaneous influences of arterial BP and PaCO2 on CBFV. When applied to experimental data, the proposed method may effectively suppress an apparently artefactual transient effect (overshoot) of CO2 changes on phase dynamics and significantly improve the prediction of CBFV below 0.1 Hz. The method thus potentially provides greater insight into the physiological processes and could lead to improved diagnostic tools in assessing autoregulation continuously even when CBFV is influenced by multiple physiological variables. Competing interest No competing interests. Ethical approval Data collection was approved by the Local Research Ethics Committee of Peking Union Medical College Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College and an

J. Liu et al. / Medical Engineering & Physics 36 (2014) 1636–1643

informed consent was signed by each participant before the study with reference number S-501. Acknowledgments The authors would like to acknowledge support from the National Natural Science Foundation of China under grant number 81000644 and the Science and Technology Innovation Commission of Shenzhen under grant number JSE201109160051A and ZDS Y20120617113312191. Conflict of interest statement No conflict of interest. References [1] Aaslid R, Lindegaard KF, Sorteberg W, Nornes H. Cerebral autoregulation dynamics in humans. Stroke 1989;20:45. [2] Ainslie PN, Celi L, McGrattan K, Peebles K, Ogoh S. Dynamic cerebral autoregulation and baroreflex sensitivity during modest and severe step changes in arterial PCO2 . Brain Res 2008;1230:115–24. [3] Birch AA, Dirnhuber MJ, Hartley-Davies R, Iannotti F, Neil-Dwyer G. Assessment of autoregulation by means of periodic changes in blood pressure. Stroke 1995;26:834. [4] Brady K, Joshi B, Zweifel C, Smielewski P, Czosnyka M, Easley RB, et al. Realtime continuous monitoring of cerebral blood flow autoregulation using nearinfrared spectroscopy in patients undergoing cardiopulmonary bypass. Stroke 2010;41:1951–6. [5] Bruns A. Fourier–Hilbert- and wavelet-based signal analysis: are they really different approaches? J Neurosci Methods 2004;137:321–32. [6] Czosnyka M, Brady K, Reinhard M, Smielewski P, Steiner LA. Monitoring of cerebrovascular autoregulation: facts, myths, and missing links. Neurocrit Care 2009;10:373–86. [7] Diehl RR, Linden D, Lücke D, Berlit P. Phase relationship between cerebral blood flow velocity and blood pressure a clinical test of autoregulation. Stroke 1995;26:1801–4. [8] Dineen NE, Brodie FG, Robinson TG, Panerai RB. Continuous estimates of dynamic cerebral autoregulation during transient hypocapnia and hypercapnia. J Appl Physiol 2010;108:604–13. [9] Edwards MR, Lin DC, Hughson RL. Modeling the interaction between perfusion pressure and CO2 on cerebral blood flow. Front Model Control Breath 2001;499:285–90. [10] Enevoldsen EM, Jensen FT. Autoregulation and CO2 responses of cerebral blood flow in patients with acute severe head injury. J Neurosurg 1978;48:689– 703. [11] Gommer ED, Shijaku E, Mess WH, Reulen JP. Dynamic cerebral autoregulation: different signal processing methods without influence on results and reproducibility. Med Biol Eng Comput 2010;48:1243–50. [12] Hu K, Lo MT, Peng CK, Liu Y, Novak V. A nonlinear dynamic approach reveals a long-term stroke effect on cerebral blood flow regulation at multiple time scales. PLoS Comput Biol 2012;8:e1002601. [13] Laback B, Balazs P, Necciari T, Savel S, Ystad S, Meunier S, et al. Additivity of nonsimultaneous masking for short Gaussian-shaped sinusoids. J Acoust Soc Am 2011;129:888. [14] Kostoglou K, Debert CT, Poulin MJ, Mitsis GD. Nonstationary multivariate modeling of cerebral autoregulation during hypercapnia. Med Eng Phys 2013;36(5):592–600.

1643

[15] Latka M, Turalska M, Glaubic-Latka M, Kolodziej W, Latka D, West B. Phase dynamics in cerebral autoregulation. Am J Physiol Heart Circ Physiol 2005;289:H2272–9. [16] Liu J, Simpson MD, Yan J, Allen R. Tracking time-varying cerebral autoregulation in response to changes in respiratory PaCO2 . Physiol Meas 2010;31:1291. [17] Marmarelis VZ, Shin DC, Orme ME, Zhang R. Closed-loop dynamic modeling of cerebral hemodynamics. Ann Biomed Eng 2013;41(5):1029–48. [18] Marmarelis VZ, Shin DC, Orme ME, Zhang R. Time-varying modeling of cerebral hemodynamics. IEEE Trans Biomed Eng 2014;61(3):694–704. [19] Marple Jr SL. Digital spectral analysis with applications, 1. Englewood Cliffs, NJ: Prentice-Hall; 1987. p. 512. [20] Mitsis GD, Zhang R, Levine BD, Marmarelis VZ. Modeling of nonlinear physiological systems with fast and slow dynamics II Application to cerebral autoregulation. Ann Biomed Eng 2002;30:555–65. [21] Mitsis GD, Zhang R, Levine BD, Marmarelis VZ. Cerebral hemodynamics during orthostatic stress assessed by nonlinear modeling. J Appl Physiol 2006;101:354–66. [22] Novak V, Yang AC, Lepicovsky L, Goldberger AL, Lipsitz LA, Peng CK. Multimodal pressure-flow method to assess dynamics of cerebral autoregulation in stroke and hypertension. Biomed Eng Online 2004;3:39. [23] Panerai RB, Kelsall AWR, Rennie JM, Evans DH. Analysis of cerebral blood flow autoregulation in neonates. IEEE Trans Biomed Eng 1996;43:779–88. [24] Panerai RB, Simpson DM, Deverson ST, Mahony P, Hayes P, Evans D. Multivariate dynamic analysis of cerebral blood flow regulation in humans. IEEE Trans Biomed Eng 2000;47:419–23. [25] Panerai RB, Carey BJ, Potter JF. Short-term variability of cerebral blood flow velocity responses to arterial blood pressure transients. Ultrasound Med Biol 2003;29:31–8. [26] Panerai RB, Eames P, Potter JF. Variability of time-domain indices of dynamic cerebral autoregulation. Physiol Meas 2003;24:367. [27] Panerai RB, Eames PJ, Potter JF. Multiple coherence of cerebral blood flow velocity in humans. Am J Physiol Heart Circ Physiol 2006;291:H251–9. [28] Panerai RB. Cerebral autoregulation: from models to clinical applications. Cardiovasc Eng 2008;8:42–59. [29] Panerai RB, Dineen NE, Brodie FG, Robinson TG. Spontaneous fluctuations in cerebral blood flow regulation: contribution of PaCO2 . J Appl Physiol 2010;109:1860–8. [30] Panerai RB. Nonstationarity of dynamic cerebral autoregulation. Med Eng Phys 2013;36(5):576–84. [31] Peng T, Rowley AB, Ainslie PN, Poulin MJ, Payne SJ. Multivariate system identification for cerebral autoregulation. Ann Biomed Eng 2008;36:308–20. [32] Peng T, Rowley AB, Ainslie PN, Poulin MJ, Payne SJ. Wavelet phase synchronization analysis of cerebral blood flow autoregulation. IEEE Trans Biomed Eng 2010;57:960–8. [33] Poulin MJ, Liang PJ, Robbins PA. Dynamics of the cerebral blood flow response to step changes in end-tidal PCO2 and PO2 in humans. J Appl Physiol 1996;81:1084. [34] Simpson DM, Panerai RB, Evans DH, Garnham J, Naylor AR, Bell PRF. Estimating normal and pathological dynamic responses in cerebral blood flow velocity to step changes in end-tidal pCO2 . Med Biol Eng Comput 2000;38:535–9. [35] Steiner LA, Czosnyka M, Piechnik SK, Smielewski P, Chatfield D, Menon DK, et al. Continuous monitoring of cerebrovascular pressure reactivity allows determination of optimal cerebral perfusion pressure in patients with traumatic brain injury. Crit Care Med 2002;30:733–8. [36] Tiecks FP, Lam AM, Aaslid R, Newell DW. Comparison of static and dynamic cerebral autoregulation measurements. Stroke 1995;26:1014–9. [37] Ursino M, TerMinassian A, Lodi CA, Beydon L. Cerebral hemodynamics during arterial and CO2 pressure changes: in vivo prediction by a mathematical model. Am J Physiol Heart Circ Physiol 2000;279:H2439–55. [38] Zhang R, Zuckerman JH, Giller CA, Levine BD. Transfer function analysis of dynamic cerebral autoregulation in humans. Am J Physiol Heart Circ Physiol 1998;274:H233–41. [39] Zweifel C, Castellani G, Czosnyka M, Carrera E, Brady KM, Kirkpatrick PJ, et al. Continuous assessment of cerebral autoregulation with near-infrared spectroscopy in adults after subarachnoid hemorrhage. Stroke 2010;41:1963–8.

Rapid pressure-to-flow dynamics of cerebral autoregulation induced by instantaneous changes of arterial CO2.

Continuous assessment of CA is desirable in a number of clinical conditions, where cerebral hemodynamics may change within relatively short periods. I...
773KB Sizes 1 Downloads 4 Views