Towards robotized beating heart TECABG: Assessment of the heart dynamics using high-speed vision

Abstract

Active robotic filtering is probably the solution for beating heart Totally Endoscopic Coronary Artery Bypass Grafting (TECABG). In this work, we assess the heart motion dynamics by simultaneous use of high-speed imaging of optical markers attached to the heart, ECG signals and ventilator airflow acquisitions. Our goal is to assess the heart motions (shape, velocity, acceleration) in order to be able to make more accurate specifications for a novel, dedicated robot that could follow these motions in real time. Furthermore, using two additional inputs (ECG and airflow), we propose a novel robust prediction algorithm that could be used with a predictive control algorithm to improve the tracking accuracy.

• Heart motion
• TECABG
• arrhythmia
• motion compensation
• robotics

Introduction

The heartbeat is a major disturbance affecting medical procedures. Its effect is localized in the chest and parts of the upper abdomen, with smaller cardioballistic effects elsewhere in the body as the blood is pumped through the arteries. The heart itself undergoes two different types of quasi-periodic motion; one due to respiration, and the other due to its own beating.

Current techniques for counteracting quasi-periodic biological motion during surgery typically consist of either passive suppression or simply shutting down the source of disturbance. In cardiac surgery, the heart is often stopped and cardiopulmonary bypass used. This extracorporeal circulation has several major drawbacks, of which the most critical is a higher risk of neurological complications [1].

To avoid this problem, passive stabilizers have been developed which use pressure or suction in an attempt to immobilize a small portion of the myocardium [2]. Unfortunately, these systems have significant residual motion, especially in minimally invasive surgery [3].

Beating heart Totally Endoscopic Coronary Artery Bypass Grafting (TECABG) is probably the ultimate goal for coronary surgery: the major benefits are a shorter recovery period, less pain and fewer complications. In order to achieve this long-term goal for all patients and all types of coronary pathology, however, major technological breakthroughs are needed in order to provide more effective solutions to the problem of residual motion. Robotics is probably the key to the problem.

In general, most surgical robotic systems developed to date operate in either a teleoperative or collaborative mode, keeping the surgeon in the loop. Active physiological motion compensation can be seen as an autonomous mode, transparent to the surgeon, that uses an active mechanical device, i.e., a robot, to cancel unwanted motion while accurately following the command input given by the surgeon [4]. This active robotized cancellation device relies on an accurate predictive model of the heart motion [5]. Indeed, this model combined with a predictive control scheme has been proven to drastically reduce the cancellation error [6]. Furthermore, using additional measurements such as ECG or ventilator pressure signals allows improvement in the robustness of the motion estimation, e.g., in case of occlusions [7].

The motion of the heart is exceedingly complex: it is the result of the superposition of two non-harmonic motions: respiration and beating. This motion has been extensively studied to enhance the sharpness of medical imaging of the heart, e.g., for coronary MR angiography [9] or for the per-operative correction of X-ray fluoroscopy images of the heart [8]. It has been shown that the respiration component exhibits some hysteresis [9], that the rotation is significant, and that the heart undergoes some deformation correlated to the respiratory cycle that modulates the beating component [10].

In this work, we studied the model of myocardial motion on an anesthetized pig. The pig heart has been proven to be a good model for the human heart [11]. A high-speed camera was used to measure the motion of the myocardium with respect to 6 degrees of freedom at a refresh rate of 500 fps. ECG signals and ventilator airflow were simultaneously acquired at the same sampling rate (500 Hz).

An extensive assessment of the myocardium dynamics (3D path, rotation, velocity, acceleration) is presented. We were able to measure accelerations up to 0.6 G thanks to the high sampling rate of the camera. This assessment is a first step towards the design of a dedicated robotized system that could follow the heart surface in real time.

Since the accuracy of such a robotized system relies on the accuracy of a motion predictor, we focused our efforts on improving a frequency-based algorithm previously presented in reference [6]. In the proposed approach, the motion prediction model is computed online with an adaptive algorithm that can cope with variations in the cardiac and respiratory frequencies and modifications in the amplitude or shape of the two components. Furthermore, our model takes into account the modulation of the beating component due to respiration, thanks to a Linear Parameter Varying (LPV) technique. This is the main contribution of this work when compared to other approaches that neglect this modulation.

Experimental results show that the motion can be predicted with good accuracy several seconds in advance. In the case of arrhythmias, although the shape of the motion can no longer be predicted, we demonstrate that the occurrence of the arrhythmic behavior can be anticipated several tenths of a millisecond in advance thanks to ECG monitoring.

In the next section, we present the materials and methods used in this work. This is followed by a description of the assessment of the heart motion dynamics. Finally, a prediction algorithm based on the LPV technique is presented and validated by experimental results in vivo.

Materials and methods

Experiments were carried out on a pig that underwent a sternotomy with general anesthesia. After positioning the chest retractor, the left anterior descending (LAD) coronary artery was stabilized using an Octopus v4.3 tissue stabilizer (Medtronic, Minneapolis, MN). A rotating knob allows modification of the stiffness of the arm of this device, thus freeing or constraining local motion of the heart.

To measure heart motion, seven LEDs were attached to the stabilizer's suction fingers to serve as optical markers. A 500 Hz high-speed camera with a 256 × 256 pixel grayscale sensor (DALSA CAD6) was placed on a tripod, with its lens focused on the stabilized area (Figure 1). The 3D position of the Octopus with respect to the camera could then be computed from the marker positions in the image using the modified version of the DeMenthon algorithm for coplanar points [12]. Since the relative motion between the myocardium and the optical markers was negligible, the motion of the heart surface could be reconstructed from the motion of the markers.

Figure 1. The heart stabilizer OCTOPUS and its visual markers. [Color version available online.]

ECG signals were acquired through a classical 3-lead ECG cable and a self-made differential amplifier with a gain of 1000. We built this simple amplifier based on the AD624 instrumental amplifier (analog device) to avoid the 20–25 ms delay found in commercial ECGs due to built-in signal postprocessing.

Two AWM700 airflow sensors from Honeywell were used for real-time measurement of the ventilator flow (Figure 2). This uni-directional sensor, which is specially designed for biomedical use, has a 6-millisecond response time.

Figure 2. The two uni-directional air flowmeters (Honeywell, AWM700). [Color version available online.]

Both ECG signals and airflow measurements were then acquired at 500 Hz by a PCI acquisition board synchronized with the image acquisition. The whole acquisition software runs on RTAI, a real-time operating system, to ensure perfect synchronization and minimal jitter.

Assessment of the myocardium motion

Motion without ventilation

Stabilized residual motion

As expected, the amplitude of the heart motion was highly dampened by the stabilizer. Three-dimensional trajectories of 10 consecutive heart beats with and without stabilization were extracted from camera information (Figures 3 and 4). It appears that the area of interest (AOI) excursion range is approximately 2.3 mm in the frontal plane and 15 mm in the sagittal plane for the free motion. With a constrained motion (maximum stiffness of the octopus arm), the excursion range is smaller: 500 µm and 6 mm, respectively. These values for the excursion range are similar to those reported by Cattin et al. [13] for the residual motion of a stabilized AOI with constant hemodynamics.

Figure 3. Trajectories of the AOI in 3D space without ventilation (time parametrization, depth z coded in color intensity). [Color version available online.]

Figure 4. Trajectories of the AOI in 3D space with ventilation (time parametrization, depth z coded in color intensity). [Color version available online.]

Repeatability of the beating trajectory

Figures 3 and 5 show that the 3D path due to cardiac beating is very repetitive. The maximum deviation and standard deviation around a mean trajectory for the free and stabilized motion are reported in Table I.

Figure 5. Translational and angular motions with respect to cardiac phase χ (IS = inferior-superior, RL = right-left, PA = posterior-anterior.) [Color version available online.]

Table I.  Maximum and standard deviations with respect to the mean trajectory of 15 consecutive heartbeats.

	Frontal plane		Sagittal plane
	Max deviation	Std deviation	Max deviation	Std deviation
Free motion	146.6 µm	25.6 µm	1.6mm	0.39mm
Stabilized	63.7 µm	14.7 µm	1.5mm	0.4mm

The repeatability of heart beating trajectories in the absence of ventilation is a very interesting feature. Indeed, considering this almost perfect periodicity, a repetitive control scheme (an example of which is given in reference [6]) should give really good results if it is assumed that respiration can be stopped for a short time. The standard and maximal deviation then provide a clue to the accuracy that such a system could achieve.

Velocity and acceleration assessment

The heart local velocity is plotted in Figure 6. It is obtained by deriving the position with respect to time. The velocities of 6 consecutive heartbeats are plotted with respect to the ECG phase for two pigs. The two plots depict the magnitude of the translational and rotational instantaneous velocity vector. Maximum velocities are 140 mm/s for translation and 60 degrees/s for rotation.

Figure 6. AOI translational and rotational velocity magnitudes with respect to the ECG phase. [Color version available online.]

Since we acquire the position at a very high frequency, we are able to estimate the acceleration by twice deriving the position with respect to time without too much alteration of the signal/noise ratio, as shown in Figure 7. The maximum linear and angular accelerations are 0.57 m · s^{− 2} and 40 degrees · s^{− 2} respectively. These values will be useful when writing the specification for an active robotic stabilizer.

Figure 7. AOI translational and rotational acceleration magnitudes with respect to the ECG phase. [Color version available online.]

In pig 2, we tested the effect of injection of adrenaline. The cardiac frequency jumped to 200 bpm and we measured maximum velocities and accelerations of almost twice the values recorded at 90 bpm.

Arrhythmias

In Figures 8 and 9, we assess the possibility of predicting arrhythmic behavior of the heart. The vertical cursors in both plots indicate the time at which the error between arrhythmic ECG (respiratory motion) (solid lines) and normal ECG (respiratory motion) (dashed lines) exceeded 10% of the peak-to-peak amplitude. This figure shows that abnormal ECG could be detected approximately 90 ms before the abnormal motion (mean delay = 80 ms; standard deviation = 32 ms on 12 arrhythmias), allowing sufficient time for a robotic system to switch to failsafe mode.

Figure 8. Arrhythmia 1. [Color version available online.]

Figure 9. Arrhythmia 2. [Color version available online.]

Motion with ventilation

Influence of ventilation on heart motion

During respiration, variations in the lung volume result in motion of the heart within the chest. In Figure 4, the motion of the AOI and the corresponding lung volume are recorded over a full respiration cycle.

The low frequency component of the heart AOI appears – as expected – to be directly correlated with the lung volume (see Figure 10). Note the different paths followed by the AOI during forced inspiration and free expiration. Moreover, for low lung volume, i.e., at the end of expiration, the 3D trajectory of the heart AOI during a heartbeat is similar to that observed in the total absence of ventilation (Figure 3).

Figure 10. Motion in the frontal plane, ECG and lung volume with respect to time. [Color version available online.]

Extraction of the respiratory component

ECG signal acquisition can be used as a clock to sample the heart AOI position when, for example, the heart is “at rest”. This way, it is possible to suppress the heart's beating motion and have access to the respiration component alone. This technique called “gating” is often used in medical imaging to suppress motion blur (see, for example, reference [8]).

The detection of the QRS complex can be performed online by derivation and adaptive thresholding of the ECG signal. From the original ECG signal, a discrete clock signal QRS[k] with impulses corresponding to QRS detection time is created. In order to improve precision, the previous clock can be delayed by half a cardiac period (χ = 0.5) to sample the respiration while the heart is at rest (see Figure 5). By interpolating between the respiratory samples with a spline, the respiration motion is reconstructed and subtracted from the global motion in order to isolate the heart's beating motion.

The upper plot in Figure 11 shows the reconstructed respiratory component and the lower plot shows the isolated beating component. The resulting beating motion is almost identical to that in the absence of ventilation. However, the plot shows clearly that it is modulated by the lung volume. Therefore, we cannot assume that the beating is decoupled from respiration. QRS sampling is compared to the Fourier Linear Combiner (FLC) technique presented in reference [5]. Note that the modulation is more significant when using the QRS sampling method.

Figure 11. Extracted respiratory motion with QRS sampling and FLC. [Color version available online.]

Figure 12 corroborates the observations of Manke et al. [9] regarding hysteresis in the ventilation cycle and also those of McLeish et al. [10] concerning the significance of heart rotation due to respiration.

Figure 12. Extracted respiratory motion in the IS, RL and PA planes. [Color version available online.]

A heart motion predictor based on biological signals

Description of the model

The proposed model is based on the extraction of the respiration component presented earlier. The integrated airflow – and thus the air volume – and QRS occurrences of the ECG signal are both used. The proposed algorithm is a two-stage identification:

• 1) In a first step, the profile of the respiration component r is extracted with the approach described in the previous section. Let k be the current sample number and T_oυ be the sample number corresponding to the beginning of a new respiration cycle. It can be easily detected by a simple thresholding of the volume information. Then,where Fr is a cubic smoothing spline function which interpolates one period of the respiration component.

• 2) In a second step, the heartbeat component is extracted from the global motion by subtracting the respiration component. As shown previously, the heart beating shape depends on lung volume. In order to take this coupling into account, a Linear Parameter Varying (LPV) [14] Finite Impulse Response (FIR) model is proposed. The model input signal is an impulse at each QRS occurrence and its parameters are variable, depending on lung volume and being adaptively refreshed to match the current measurements.

Let υ[k] be the lung volume and QRS[k] the impulse signal corresponding to the QRS occurrence. Then, the heart position due to cardiac beating b[k] is given bywhere d is the number of samples in one cardiac period.

With this formulation, the coefficients a_i0 represent the mean value of the beating position, i samples after the previous QRS complex, and a_i1 is the variation of this position related to lung volume. Identification of the coefficients {(a_i0, a_i1), i = 1 … d} can be performed online via a recursive least mean square (RLS) algorithm for LPV systems, where the criterion of minimization is the error e[k] between the model output and the reference signal [14].

This identification scheme is summarized in Figure 13. The global motion h[k] of the heart AOI at sample k is given by

Figure 13. Block diagram of the motion predictor.

Validation and comparison with other approaches

The output of the proposed model is used to predict the motion of the heart several steps in advance. We compared the accuracy of these predictions to that with other methods found in the literature. We used a set of data down-sampled at 50 Hz with a ventilation frequency of f_r = 0.242 Hz and a cardiac frequency f_b = 1.23 Hz.

In Figure 14, we predict the motion Δ = T_b/2 = 20 samples in advance, i.e., half a heartbeat in advance (χ = 0.5, i.e., 0.4 s), whereas in Figure 15 the prediction is 40 samples in advance (χ = 1, i.e., 0.8 s). The curve in black gives the true measured motion of the heart.

Figure 14. Predictor output for 0.5_χ anticipation. [Color version available online.]

Figure 15. Predictor output for 1_χ anticipation. [Color version available online.]

These plots compare the prediction accuracy of a frequency-based method, i.e., the FLC technique used by Thakral et al. [5], which is very similar to the adaptive filtering technique presented in reference [6], to the QRS-sampling + LPV technique. Frequency-based methods are based on the assumption that the extracted signals, beating motion and respiratory motion, remain quasi-periodic from one period to the next. As a result, the modulation of the beating component is neglected. With our approach, the system identifies adaptively the coupling between respiration and cardiac beating and uses this relationship to enhance the accuracy of the prediction, as shown in these figures.

In Table II, we compare the mean prediction error for a prediction in the future with three different algorithms. For all tests, we waited 10 respiratory cycles to allow the algorithms to converge, then computed the average error over 4 respiration cycles.

Table II.  Comparison of the mean prediction error.

	Translations			Rotations
Method	Tx (IS)	Ty (RL)	Tz (PA)	Rx (IS)	Ry (RL)	Rz (PA)
FLC	0.1017mm	0.1238mm	0.9212mm	0.0055rad	0.0027rad	0.0037rad
QRS	0.1045mm	0.1455mm	0.5528mm	0.0052rad	0.0022rad	0.0036rad
QRS + LPV	0.0453mm	0.0810mm	0.5155mm	0.0030rad	0.0017rad	0.0018rad

In Table II, “FLC” is the method used in reference [5], being first introduced by Riviere et al. [15]. “QRS” is a method that uses QRS sampling to extract the respiratory motion, but the modulation effect of the beating motion is neglected by the predictor. In “QRS + LPV”, QRS sampling is used to extract the respiratory component and LPV is used to take into account the modulation of the beating component due to respiration. These data show clearly the improvement achieved with this last method.

Discussion

This study shows that there is a coupling between the two components of heart motion: the shape of the beating is modulated by the state of the respiration. We propose a new approach to predict the motion of the heart that takes this coupling into account. Experimental results show improved accuracy with this motion predictor, as compared to other approaches found in the literature.

Furthermore, an assessment of the heart dynamics and their relationship with ECG signals shows that the cardiac beating component exhibits very sharp transients that yield important accelerations, and that arrhythmic motions can be predicted approximately 80 ms in advance by monitoring the ECG.

These results are the starting point for the design of an active robotized filtering system that could synchronously follow the surface of the heart and allow minimally invasive TECABG for almost every patient and coronary pathology.