Acta Polytechnica doi:10.14311/AP.2013.53.0854 Acta Polytechnica 53(6):854–861, 2013 © Czech Technical University in Prague, 2013 available online at http://ojs.cvut.cz/ojs/index.php/ap QUANTIFICATION OF RESPIRATORY SINUS ARRHYTHMIA WITH HIGH-FRAMERATE ELECTRICAL IMPEDANCE TOMOGRAPHY Christoph Hoog Antink∗, Steffen Leonhardt Philips Chair for Medical Information Technology, Helmholtz-Institute for Biomedical Engineering, RWTH Aachen University, Pauwelsstr. 20, 52074 Aachen, Germany ∗ corresponding author: hoog.antink@hia.rwth-aachen.de Abstract. Respiratory Sinus Arrhythmia, the variation in the heart rate synchronized with the breathing cycle, forms an interconnection between cardiac-related and respiratory-related signals. It can be used by itself for diagnostic purposes, or by exploiting the redundancies it creates, for example by extracting respiratory rate from an electrocardiogram (ECG). To perform quantitative analysis and patient specific modeling, however, simultaneous information about ventilation as well as cardiac activity needs to be recorded and analyzed. The recent advent of medically approved Electrical Impedance Tomography (EIT) devices capable of recording up to 50 frames per second facilitates the application of this technology. This paper presents the automated selection of a cardiac-related signal from EIT data and quantitative analysis of this signal. It is demonstrated that beat-to-beat intervals can be extracted with a median absolute error below 20 ms. A comparison between ECG and EIT data shows a variation in peak delay time that requires further analysis. Finally, the known coupling of heart rate variability and tidal volume can be shown and quantified using global impedance as a surrogate for tidal volume. Keywords: Electrical impedance tomography, respiratory sinus arrhythmia, cardio-respiratory cou- pling, heart-rate variability. 1. Introduction Electrical Impedance Tomography (EIT) is a powerful imaging tool. It seeks to reconstruct the impedance distribution inside a patient from measurements on the boundary. These measurements are non-invasive, painless and have no known side effects. This makes EIT an ideal tool for long-term monitoring of patients. Since the electrical impedance of lung tissue varies greatly with air content, the most common use for med- ical EIT is in pulmonary monitoring. Here it serves to visualize and analyze the regional distribution of ventilation, which in turn can be used for example to automatically optimize the respirator settings for mechanically ventilated patients [1]. Although the most prominent changes in conductivity in the thorax originate from respiration, a signal that is roughly an order of magnitude smaller can be attributed to cardiac activity. A question still unanswered is the op- timal electrode configuration to maximize the quality of the cardiac related signal [2]. In mechanically ventilated patients, the distribution of the ventilation is a very important parameter. If water accumulates in the dorsal lung regions, these regions may collapse and thus may not be ventilated at all. At the same time, the ventral regions may be over-ventilated, leading to potentially lethal lung damage. However, distribution of ventilation is not the only parameter of importance, perfusion is as well. If the pressure is set very high to optimize the ventilation, this may simultaneously hinder cardiac activity, thus leading to a reduction in gas exchange. Monitoring cardiac functionality is therefore a task of equal importance. In addition to heart rate (HR), heart rate variability (HRV), i.e. the change in timespan between two con- secutive heartbeats, has received increasing attention recently. It has been examined as an indicator for stress [3], for sleep stages [4] and even as a predictor for septic shock [5]. Respiratory sinus arrhythmia (RSA) is an oscillation of the heart rate in sync with respiration. In the inspiration phase, an increase in heart rate can be observed, whereas a decrease in heart rate occurs in the expiration phase. This fact can be exploited to extract the respiratory rate (RR) from an electrocardiogram (ECG) [6]. At the same time, if an individual model for the cardio-respiratory coupling can be calculated or learned, the breathing signal could be used to improve the estimation of the cardiac related signal in a multi-sensor data fusion setting. EIT is widely used in medical research and its intro- duction into standard clinical practice is now immi- nent. The world’s first commercial medical EIT device was released in 2011 by Dräger Medical GmbH, under the name PulmoVista® 500 [7]. This device works with a 16-electrode belt and is capable of recording up to 50 frames per second (FPS). With such a high frame rate, the noninvasive nature of the measurement, and the possibility to measure lung- and cardiac-related signals simultaneously, the device seems an ideal tool for analyzing cardio-respiratory coupling. This paper presents a proof-of-concept study con- 854 http://dx.doi.org/10.14311/AP.2013.53.0854 http://ojs.cvut.cz/ojs/index.php/ap vol. 53 no. 6/2013 Quantification of RSA with High-Framerate EIT Figure 1. Electrode position in the experiment: A EIT electrode belt, B EIT reference electrode, C1−3 ECG electrodes. sisting of simultaneous EIT and ECG recordings. In Section 2, technical details are presented. Section 3 presents the results and a discussion, followed by the conclusions in Section 4. 2. Materials and methods 2.1. Hardware EIT measurements were obtained using the Pulmo- Vista® 500 EIT device, operating at a rate of 50 FPS. Additionally, a single channel ECG was recorded us- ing the SOMNOlab 2 PSG Standard by Weinmann Medical Technology. Designed for polysomnography studies, this device features a variety of sensors. For this study, only ECG functionality sampling at 256 Hz was used. 2.2. Trial setup The study was conducted as a self-experiment. No recent history of pulmonary or cardiac related diseases is known. Electrodes were positioned according to Figure 1. Ten measurements were conducted in sitting position: In the first 8 runs, the respiratory frequency was controlled to a specific value and EIT data was recorded for two minutes, see also Table 1. In run 9, breath was held after expiration for as long as possible (25 seconds) and in run 10 the same was done after inspiration 2.3. Data preprocessing Data was recorded and exported using the manufac- turers’ respective software tools, and was imported into MATLAB. In each trial, a single EIT file was HF/LF Ratio [A .U .] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Figure 2. Ratio of HF to LF components in the EIT signal. The heart region can be clearly distinguished from the background. recorded, ECG was recorded continuously. The ”for- ward problem” in EIT is governed by the equations ∇· ~J(~x) = −∇·σel(~x)∇φ(~x) = 0 (1) for ~x ∈ Ω \∂Ω and j(~x) = −~J(~x) ·~n(~x) = σel(~x)∇φ(~x) ·~n(~x) (2) for ~x ∈ ∂Ω, see also [8]. Here, ~x is a coordinate inside a body Ω or on its boundary ∂Ω with normal vector ~n, ~J is current density, φ is electric potential, j is scalar current density on the boundary, and σel is conductivity distribution inside the body. If j(~x) and σel(~x) are known, the partial differential equation 1 has a unique solution φ(~x). The ”inverse problem” in EIT is where the underlying conductivity distribu- tion σel(~x) is unknown and needs to be inferred from current injections and voltage measurements on the boundary. This is a much harder problem to solve due to its ill-posed nature, which means that at least one of the criteria • existence of a solution, • uniqueness of the solution or • continuity of the solution is violated. To be still able to solve the problem, a-priori knowledge in the form of regularization is in- corporated. This usually results in a spatially smooth reconstructed distribution σel. The EIT system used here works with an electrode belt with N = 16 electrodes, which leads to N(N − 3) = 208 volt- age measurements per frame, of which due to reci- procity only 104 are linearly independent. The de- vice’s internal linear method of reconstruction was used, leading to a sequence of impedance images Z with 32 by 32 pixels at 50 FPS. To extract the cardiac signal, a single pixel Z was selected. In a first proof of concept, this selection process was 855 C. Hoog Antink, S. Leonhardt Acta Polytechnica 0 20 40 60 80 100 120 −0.5 0 0.5 1 1.5 G lo b a l Im p e d a n c e [A .U .] Run # 1 0 20 40 60 80 100 120 −5 0 5 C a rd ia c P ix e l [A .U .] Time [seconds] Figure 3. Time course of the global impedance (top) and of a single pixel in the cardiac region (bottom). performed manually and yielded reasonable results [9]. As a subsequent development, the selection pro- cess was optimized: For this, the time course of each pixel of the first trial was converted into fre- quency domain using the Fast Fourier Transforma- tion (FFT). Next, the ratio of high frequency (HF, 0.99 to 1.12 Hz, cardiac activity) to low frequency components (LF, 0.08 to 0.16 Hz, respiratory activ- ity) was calculated. To avoid division by (nearly) zero, a small offset was added to the low frequency component. In Figure 2, the HF/LF ratio is spa- tially visualized, clearly distinguishing the heart re- gion from the background. The coordinates of the maximum value were used in the subsequent exper- iments. No spatial averaging was performed, since the reconstructed EIT impedance image data is in- trinsically ”spatially low-pass filtered”, as argued above. The difference between the time course of the global impedance (GI) and the selected pixel is visualized in Figure 3. While the average impedance shows changes mainly due to respiration, a strong cardiac related signal (CRS) can be deduced in the time course of a single pixel. This signal was re- sampled to 256 Hz to match the ECG sampling rate and was high-pass filtered to remove the residual res- piratory component, resulting in CRS∗. The same high-pass was applied to the ECG signal ECG′, re- sulting in ECG. CRS∗ and ECG were then cross- correlated; a very distinct peak could be observed that was used to calculate tsync, which could then be used to synchronize the EIT and ECG data by selecting the appropriate part in the ECG stream. To remove high-frequency artifacts, CRS∗ was low- pass filtered before peak detection. No further pro- cessing was applied to ECG. The algorithm is out- lined in Figure 4. The resulting arrays tP eak,EIT and tP eak,ECG were used in the subsequent analy- sis. Pixel Selection EIT-Device Resampling High-Pass Z 50Hz [32 x 32] Z 50Hz Z‘ 256Hz CRS* ECG‘ 256Hz ECG Cross Correlation Low-Pass Peak Detection CRS ECG-Device High-Pass Peak Detection tPeak, EIT tsync Max tPeak, ECG Figure 4. Outline of the algorithm to extract heart- beats from EIT and ECG as well as synchronizing the raw data. 2.4. Data analysis First, beat-to-beat intervals (BBI) were calculated from the arrays containing the peaks. Since these time intervals are by definition unevenly sampled and thus located on an irregular grid, the BBI de- rived from CRS were linearly interpolated to the lo- cations of the BBI obtained from ECG. Then, the root mean square error (RMSE), the median absolute error (MAE) and correlation coefficient r were cal- culated. To validate the RR, the global impedance signal was evaluated in the frequency domain, using the FFT. To evaluate RSA, the Lomb periodogram 856 vol. 53 no. 6/2013 Quantification of RSA with High-Framerate EIT 0 20 40 60 80 100 120 −0.5 0 0.5 1 G lo b a l Im p e d a n c e [A .U .] Run # 5, RR = 20 BPM 0 20 40 60 80 100 120 400 600 800 1000 B e a t− to − B e a t In te rv a l [m s ] RMSE = 23.0382 ms, r = 0.91413 CRS ECG 0 20 40 60 80 100 120 150 200 250 300 P e a k D e la y [m s ] Time [seconds] t delay,mean = 203.683 ms ECG to EIT Figure 5. Time course for run # 5. 0 20 40 60 80 100 120 −0.5 0 0.5 1 1.5 G lo b a l Im p e d a n c e [A .U .] Run # 7, RR = 12 BPM 0 20 40 60 80 100 120 400 600 800 1000 1200 B e a t− to − B e a t In te rv a l [m s ] RMSE = 33.4943 ms, r = 0.9393 CRS ECG 0 20 40 60 80 100 120 150 200 250 300 350 P e a k D e la y [m s ] Time [seconds] t delay,mean = 218.0469 ms ECG to EIT Figure 6. Time course for run # 7. was calculated. This allows the use of BBI directly without interpolating them to the regular grid. The respiratory rate was identified as the peak in the peri- odogram above 1 BPM. The heart rate was calculated as the inverse of the median BBI. Additionally, the delay between peaks derived via ECG and EIT was calculated. Here, values below 150 and above 400 milliseconds were considered to be outliers, and were discharged. It is known that the tidal volume (TV) influences the degree of RSA: It was found in [10], that the change in BBI shows an almost linear dependence on TV. Additionally, the coefficient was found to be frequency dependent, being almost constant up to a corner frequency, which was found to be different for each subject, with a mean value of 7.1±1.5 BPM. For higher breathing frequencies, it showed an exponen- tial roll-off of around 20.4 ± 2.4 dB per decade ([10], Table 2). In this study, only the breathing frequency was actively controlled with the help of an acoustic and visual signal. No such device was employed to keep TV constant. Since the global impedance correlates very well with TV, it was used as a surrogate. Visual 857 C. Hoog Antink, S. Leonhardt Acta Polytechnica 0 5 10 15 20 25 30 35 40 −1 0 1 2 G lo b a l Im p e d a n c e [A .U .] Run # 10, Breath Hold (Inspiration) 0 5 10 15 20 25 30 35 40 600 800 1000 1200 B e a t− to − B e a t In te rv a l [m s ] RMSE = 18.3771 ms, r = 0.98304 CRS ECG 0 5 10 15 20 25 30 35 40 150 200 250 300 P e a k D e la y [m s ] Time [seconds] t delay,mean = 223.0569 ms ECG to EIT Figure 7. Time course for run # 10. Run RR RREIT RRCRS RRECG # [BPM] [BPM] [BPM] [BPM] 1 6 6.00 5.99 5.95 2 10 10.00 10.00 10.00 3 14 14.00 14.01 14.05 4 18 18.00 18.09 18.09 5 20 20.00 20.11 20.10 6 16 16.00 15.99 16.00 7 12 12.00 12.05 12.05 8 8 8.00 7.99 7.99 Table 1. Respiratory rate according to the protocol (RR), calculated from EIT (RREIT) and calculated via BBI periodogram analysis from CRS (RRCRS) and ECG (RRECG). Run 9 and 10 were breath-hold experiments. inspection of the time course of the global impedance showed a small variation in TV over time, see Figure 3. For exact analysis, all breaths were identified using a simple peak detector. Next, the difference of the maximum and the minimum in global impedance as well as the difference between maximum and mini- mum BBI in a small window around each peak were determined. 3. Results and discussion Tables 1 and 2 present the values described above. Additionally, Figures 5 to 7 present the time course of the global impedance change, BBI derived from CRS and ECG as well as the ECG to CRS peak delay. One can observe that high frame rate EIT is capable of analyzing cardiac function with great precision. Table 1 shows that the effect of RSA can clearly be inferred from the BBI signal derived via EIT, since the calcu- lated respiratory rates are almost identical. Figures 5 to 7 also show very good agreement of BBI calculated via EIT and ECG. A careful analysis of the raw data shows that relatively high RMSE stems from singular artifacts, especially from the misclassification of beats in CRS. These errors greatly influence RMSE, whereas MAE is influenced only minimally. This is especially prominent in run 9, where the correlation coefficient is below 0.7 and RMSE is above 90 ms, while the MAE shows a very low value below 10 ms, proving that singular outliers influence the result. Combin- ing all experiments, RMSE is 34.31 ms and MAE is 15.91 ms. This is a promising observation, considering that the EIT device used here has a frame rate of 50 FPS, resulting in a sampling time of 20 ms. An interesting observation can be made from Fig- ures 6 and 7. In general, the ECG and CRS peaks are not synchronous, since ECG stems from the electrical activity of the heart, whereas a peak in CRS stems from a maximum in contraction, i.e. a mechanical activity. The delay of roughly 200 milliseconds be- tween the ECG peak (R-wave) and the CRS peak is consistent with values in the literature, relating ECG and blood volume inside the heart [11]. It is very interesting to observe a slow drift in Figure 7, peak- ing just before expiration, i.e. just before maximum discomfort was reached. Figure 6, however, shows a sharp jump. While at first glance an artifact seems the most reasonable explanation, the consistency of the BBI derived from CRS and from ECG does not indicate an obvious source of error. 858 vol. 53 no. 6/2013 Quantification of RSA with High-Framerate EIT Run HRCRS HRECG RMSE MAE r # [BPM] [BPM] [ms] [ms] 1 69.50 69.03 28.92 12.60 0.97 2 74.56 75.29 30.14 16.98 0.96 3 77.58 79.59 37.46 23.13 0.91 4 87.77 88.79 31.39 18.13 0.86 5 85.81 85.81 23.04 15.06 0.91 6 78.37 80.21 29.68 17.20 0.94 7 75.29 76.42 33.49 21.24 0.94 8 74.38 76.04 24.89 11.72 0.96 9 65.08 67.22 98.08 7.32 0.68 10 61.69 61.69 18.38 6.14 0.98 Table 2. Heart Rate calculated from CRS (HRCRS) and from ECG (HRECG). Additionally, Root-mean-square error (RMSE), median absolute error (MAE) as well as correlation coefficient (r) of BBI calculated from ECG and CRS. Run RR std(∆GI) r(∆GI, ∆BBI) p(∆GI, ∆BBI) a b # [BPM] [A.U.] [ms/A.U.] [ms] 1 6 0.13 0.91 0.0002 871 -596 2 10 0.05 0.19 0.4445 200 -83 3 14 0.08 0.01 0.9741 3 20 4 18 0.05 0.10 0.6757 63 -29 5 20 0.04 0.22 0.2893 184 -112 6 16 0.04 -0.29 0.1593 -398 366 7 12 0.06 0.09 0.6938 58 18 8 8 0.08 0.80 0.0011 341 -223 Table 3. Respiratory rate according to protocol (RR), standard deviation of global impedance peaks (std(∆GI)) as well as the correlation coefficient (r(∆GI, ∆BBI)) of the change in GI and the change in BBI and the corresponding p-Value. a and b represent the parameters of the linear fit ∆BBI = a · ∆GI + b. An analysis of the tidal volume dependence of RSA is presented in Table 3. It can be seen that although the change in TV is relatively small, a correlation greater than 0.9 could be measured in run 1, which is the one with the lowest RR. The run with the second lowest RR (run 8) also showed a relatively high correlation. Beyond the RR of 8, no correlation could be determined. This is probably due to the low variation in TV in combination with the exponential decay of the slope parameter a discovered in [10] — see also its determined values for RR of 6 and 8 in Table 3. Figures 8 and 9 plot the time course as well as the direct comparison of the change in global impedance and the change in BBI, clearly showing the correlation. 4. Conclusions and outlook This paper has proved the general feasibility of EIT as a tool for analyzing cardio-respiratory coupling. HRV could be measured with MAE below 20 ms, allowing the calculation of RR from BBI analysis. This is a redundant task here, since global impedance can be used with less effort. However, it serves to prove the accuracy of the CRS extracted from EIT. Moreover, it shows the redundancy in the measured data that could be used in a model for multi sensor data fusion, exploiting the cardio-respiratory coupling. Interesting observations can be made on the basis of an analysis of the relationship between CRS and ECG. Slow drifts could be observed in the time interval between ECG- peak (R-wave) and CRS-peak (maximum contraction of the heart) in the breath-hold experiment. This ob- servation needs further investigation and might allow non-invasive analysis of traditional parameters such as blood pressure, or the derivation of novel surrogates for cardiac health. Finally, the known coupling of heart rate variability and tidal volume could be shown using a change in global impedance as a surrogate for 859 C. Hoog Antink, S. Leonhardt Acta Polytechnica 10 20 30 40 50 60 70 80 90 100 110 0.6 0.8 1 1.2 1.4 Run # 1, RR = 6 BPM, r = 0.91337 Time [seconds] C h a n g e i n G lo b a l Im p e d a n c e [A .U .] 10 20 30 40 50 60 70 80 90 100 110 0 100 200 300 400 C h a n g e i n B B I [m .s .] Figure 8. Time course of the change in GI and the change in BBI for run # 1. 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 0 50 100 150 200 250 300 350 400 450 Change in Global Impedance [A.U.] C h a n g e i n B B I [m .s .] Run # 1, RR = 6 BPM, r = 0.91337 Measurement Linear Fit Figure 9. The change in BBI over the change in GI for run # 1. the latter. Since this study presents only a proof of concept, some limitations do exist: Only a single subject par- ticipated in the experiment, a larger cohort needs to be examined to test robustness. Moreover, only the respiratory rate was controlled. Additionally, a sim- ple peak detector is likely to be only a suboptimal approach to identifying the true maximum of CRS. Firstly, this may lead to a systematic offset if the shape of the peak changes, which could be an explana- tion for the jump observed in Figure 6. Secondly, this leads to some misclassifications which are responsible for outliers. While small in number, their high ampli- tude leads to serious degrading in non-robust metrics like RMSE and the correlation coefficient. Finally, other sources of information, e.g. the CRS waveform or its slope should be examined. This could lead the way to new, non-invasive tools for cardio-respiratory model verification, and potentially to new methods of diagnosis. Acknowledgements The author would like to thank Dipl.-Ing. Boudewijn Venema for helping with the SOMNOlab equipment. References [1] T. Meier, H. Luepschen, J. Karsten, et al. Assessment of regional lung recruitment and derecruitment during a PEEP trial based on electrical impedance tomography. Intensive Care Medicine 34(3):543–50, 2008. [2] J. Nasehi Tehrani, C. Jin, A. L. McEwan. Modelling of an Oesophageal Electrode for Cardiac Function Tomography. Computational and Mathematical Methods in Medicine 2012:585786, 2012. [3] M. Pagani, G. Mazzuero, A. Ferrari, et al. Sympathovagal interaction during mental stress. A study using spectral analysis of heart rate variability in healthy control subjects and patients with a prior myocardial infarction. Circulation 83(4):II43–51, 1991. [4] M. H. Bonnet, D. L. Arand. Heart rate variability: sleep stage, time of night, and arousal influences. Electroencephalography and Clinical Neurophysiology 102(5):390–6, 1997. 860 vol. 53 no. 6/2013 Quantification of RSA with High-Framerate EIT [5] W.-L. Chen, C.-D. Kuo. Characteristics of heart rate variability can predict impending septic shock in emergency department patients with sepsis. Academic Emergency Medicine 14(5):392–7, 2007. [6] G. Moody, R. Mark. Derivation of respiratory signals from multi-lead ECGs. Computers in Cardiology 12:113–116, 1985. [7] E. Teschner, M. Imhoff. Elektrische Impedanztomographie: Von der Idee zur Anwendung des regionalen Beatmungsmonitorings. Dräger Medical GmbH 2010. [8] W. Lionheart, N. Polydorides, A. Borsic. The Reconstruction Problem. In D. Holder (ed.), Electrical Impedance Tomography: Methods, History and Applications, vol. 0750309520, pp. 3–64. Institute of Physics Publishing, 2004. [9] C. Hoog Antink. Analyzing Cardio-Respiratory Coupling with High-Framerate EIT : A Proof of Concept. In L. Husník (ed.), Proceedings of the 17th International Scientific Student Conference POSTER 2013. Czech Technical University in Prague, 2013. [10] J. A. Hirsch, B. Bishop. Respiratory sinus arrhythmia in humans: how breathing pattern modulates heart rate. The American Journal of Physiology 241(4):H620–9, 1981. [11] R. Klinke, H.-C. Pape, S. Silbernagl. Physiologie. Thieme, 5th edn., 2005. 861 Acta Polytechnica 53(6):854–861, 2013 1 Introduction 2 Materials and methods 2.1 Hardware 2.2 Trial setup 2.3 Data preprocessing 2.4 Data analysis 3 Results and discussion 4 Conclusions and outlook Acknowledgements References