Sage Journals: Discover world-class research

Abstract

Off-pump coronary artery bypass graft surgery outperforms the traditional on-pump surgery because the assisted robotic tools can cancel the relative motion between the beating heart and the robotic tools, which reduces post-surgery complications for patients. The challenge for the robot assisted tool when tracking the beating heart is the abrupt change caused by the nonlinear nature of heart motion and high precision surgery requirements. A characteristic analysis of 3D heart motion data through bi-spectral analysis demonstrates the quadratic nonlinearity in heart motion. Therefore, it is necessary to introduce nonlinear heart motion prediction into the motion tracking control procedures. In this paper, the heart motion tracking problem is transformed into a heart motion model following problem by including the adaptive heart motion model into the controller. Moreover, the model following algorithm with the nonlinear heart motion model embedded inside provides more accurate future reference by the quadratic term of sinusoid series, which could enhance the tracking accuracy of sharp change point and approximate the motion with sufficient detail. The experiment results indicate that the proposed algorithm outperforms the linear prediction-based model following controller in terms of tracking accuracy (root mean square).

Keywords

Heart Motion Nonlinear Prediction Quadratic Nonlinearity Heart Motion Tracking Model Following Control Robot-Assisted Beating Heart Surgery

1. Introduction

Coronary Artery Bypass Graft (CABG) is a common open-heart surgery for people who suffer severe coronary heart disease [1]. Essentially, there are two types of operation techniques: the traditional on-pump CABG and the off-pump CABG. In the traditional surgery, surgeons have to stop and replace the beating heart with a cardio-pulmonary bypass machine, in which the patient suffers a longer cognitive loss [2]. Consequently, the post-surgery complications are increased in terms of hospital time and cost [3]. The more preferred off-pump CABG has the advantage of a significant decrease in post-surgery complications. However, off-pump CABG is more demanding than the traditional way because it keeps the heart beating during the surgery. Generally, the beating heart motion is too quick for the surgeon to operate manually [4]. The advancement of robotic technology provides an effective alternative to resolve this challenge by eliminating the relative motion between the beating heart and the robotic tools to achieve a relatively stationary surgical environment.

1.1 Beating heart-assisted robotics system - introduction and requirements

As an assisted surgical robotic tool, a beating heart-assisted robotics system is designed to allow surgeons perform the operations as if the heart is stationary. Technically, the relatively stationary motion is realized through synchronizing the beating heart with surgical equipment, thus cancelling the relative motion between Point of Interest(POI) on the beating heart surface and the end-effecter of the robot. Figure 1 shows a conceptual model of the robotic tele-operational beating heart surgery system. On the right hand side, the robot follows the beating heart motion which is measured by an ultrasound sensor mounted on the POI. On the left hand side, the surgeon watches the still heart view on the screen and manipulates the joysticks to conduct surgical operations. With the assistance of the robotic tool, surgeons can possibly operate on the heart surface with high-frequency motions.

Figure 1.

Conceptual model of the robotic tele-operational system [5]

The challenges to figure out the solution to the beating heart motion tracking problem come from the requirement of safe surgery corresponding to sufficient accuracy and the irregular nature of heart motion corresponding to the fast tracking needs for the surgery. In off-pump CABG surgery, the surgeons operate on blood vessels with diameters ranging from 0.5-2mm, and therefore, for safety reasons, the dynamic position error has to be controlled to within 1% of the blood vessel diameter or in the order of 100–250μm root mean square [6]. More specifically, the tracking error of the whole robot system, including measurement error, reference signal prediction error and control error, should be accurately limited within this range. In addition, fast tracking is necessary because the motion of individual POI on the heart surface is in a relatively high bandwidth quasi-period and associated with non-stationary characteristics, therefore the robot tools need to follow the fast change POI with minimum delay. A traditional causal feedback controller alone will be replaced by the predictive controller with a prediction of heart motion in the feedforward path [6], since it cannot follow the moderate to high frequency components of the heart motion. Furthermore, the improved heart motion prediction could enhance the heart motion tracking [7].

1.2 Related Work

The research on the solution of the heart motion tracking problem extended the methods used in respiration motion cancellation and integrated many other techniques into the control strategy in order to enhance the tracking performance.

Patel et al. [8] introduced a behaviour selection-based navigation and obstacle avoidance algorithm allowing the robot to adapt itself to the surgery environment. Nakamura et al. [9] performed feedback control experiments to track the heart motion with a 4-DOF robot using a vision system to measure heart motion. The tracking error was relatively too large to perform beating heart surgery. Ortmaier [10] proposed a global robust prediction-based predictive control scheme in which the algorithm reconstructed the underlying dynamics of the heart motion through the past measurements. The prediction was based on taken theory which can predict the heart motion behaviour, even though the visual sensor was obscured for some period of time. Jin et al. [11] proposed the time delay estimation-based tracking control method to cancel nonlinearity in robot dynamics to increase the tracking accuracy. Bebek and Cavusoglu [5] presented an intelligent model predictive method using the previous circle of heart motion to predict the current circle of heart motion with the help of electrocardiogram (ECG) signals to correct the heart motion duration. A new adaptive prediction technique with the same control architecture was developed in [5] based on the autoregressive (AR) model [12,13] which was less susceptible to the variation of the heart motion period and robust to noise. However, this AR model was still a linear model that could not describe some of the nonlinear dynamics of the heart motion.

Thakral et al [14]separated he heart motion into the heartbeat and breathing components. The adaptive Fourier linear combiner (AFLC) method was introduced to model the heart motion. Ginhoux et al. [15] proposed a similar Fourier linear combiner model and moreover added disturbance modelling. In [16], the quasi-periodic beating heart motion was modelled as a time-varying dual Fourier series and an Extended Kalman Filter (EKF) was used to estimate the parameters. Yuen et al. [17,18] used a one degree-of-freedom motion compensation system to synchronize with tissue motions that may be approximated by 1D motion. A predictive filter before the PID controller was modelled by a time-varying truncated Fourier series with an offset to a few harmonics and the EKF was used to estimate the characteristics' parameters recursively.

Bachta et al. [19] reviewed the models of the heart motion and also treated the heart motion as the coupling of heartbeat and respiration. He modelled the coupling as amplitude modulation between two period signals. This model represented the coupling nature of the cardiac motion between the heartbeat and respiration motion faithfully. We can still expect further physical explanations of the coupling nature and a mathematical model that can further explore the nonlinearity of the cardiac dynamics [20].

Liang et al. [21] proposed a model following controller with linear prediction of the heart motion. The results showed that the model following controller outperformed any feedback-based controller because of the feedforward part. However, a linear prediction has limitations in representing the heart motion precisely.

In this paper, we outlined in particular the novel nonlinear prediction-based model following algorithm. The model following controller could separate the control law into the feedforward component and the feedback component, respectively, to track different frequency band signals in heart motion. The relative high frequency cardiac motion is predicted one step ahead by the inner heart motion model. Such that the model following controller, embedding the adaptive heart motion model, would allow the control effort, especially in the feedforward path, to adjust according to the predicted future reference. In addition, we also established the novel quadratic nonlinear embedded heart motion model, which extended the linear Fourier series superposition model in [14,15] by adding the nonlinear quadratic term. This nonlinear prediction method is naturally motivated by investigating the nonlinear dynamics of the heart motion. With the help of the bi-spectrum analysis, we found the cardio-respiration interaction is the main nonlinear factor, indicating that the respiration modulates the heartbeat motion, which holds the square law form.

This paper is organized as follows: in section 2 the fast nonlinear adaptive quadratic sinusoid series prediction algorithm is introduced after the analysis of the heart motion by using the bi-spectrum method. Next, the proposed nonlinear prediction-based model following algorithm is formulated. Finally the proposed algorithm is evaluated through comparative experiments using two sets of pre-recorded heart motion data.

2. Heart Motion Nonlinear Prediction

In practice, the statistics of heart motion are likely to change during surgery. Such a change stemming from the nonlinear phase-coupling [22] nature between the respiration component and the heartbeat component will result in variations in the underlying dynamics of the POI's motion. The greater the degree of phase-coupling, the more irregular the heart rate in the beating heart motion. In order to explore the effects of these slow variations on tracking performance and investigate how the adaptive algorithms will adjust to these changes, the heart motion characteristic model should be formed and the model-based prediction is made for the controller.

In this study, two distinct types of experimentally collected heart data [23] are used. The first type includes constant heart rate data, whereas the second type includes varying heart rate data.

2.1 Nonlinear dynamics and bi-spectral analysis

Prior to the bi-spectral analysis, the physical nature of the heart motion should be examined. The dynamics of the beating heart can be described as the coupled oscillators mainly between the cardiac and respiratory oscillations [20] which is formulated as

\frac{d X_{i}}{d t} = F (X_{i}) + δ F_{i} (X_{i}) + \sum_{j = 1}^{N} V_{i j} (X_{i}, X_{j}), (i = 1, 2, \dots, N)

(1)

where F(X) is the common structure for the N oscillators (N=2 for this study) and δF(X_i) is the derivative of the F(X_i), V_ij (X_i, X_j) which is the interacting term. This is the general nonlinear model for the beating heart which will represent the N dominant oscillators and their mutual modulation. The beating heart motion can be viewed as a process driven by coupled respiration motion and cardiac motion with each individual characteristic frequency. Coupled oscillators interact with each other through their amplitudes and phases. The coupled phase phenomenon is termed Phase-coupling (PC). In phase-coupling a new frequency component will be produced and its phase angles will be dependent on the phase angles of the original components [24]. Some research [25 –27] on the beating heart synchronization properties of the phase dynamics reveals the phase-coupling nature. The quadratic phase-coupling (coupling at sum and difference frequencies) occurs when a signal is passed through a square law device [24], which motivated us to model the heart motion including a quadratic term.

The assumption above can be proved by observing Figure 2 in [5] and Fig. 2 and Fig. 3 in [19] where the inter-modulation new frequency components other than dominate mode frequency components can be found. Also as notations shown in Figure 2, the newly produced frequency components are in a certain format instead of in random distribution. In order to further reveal the nonlinear phase-coupling nature of the beating heart motion, we use a more suitable higher order statistics bi-spectrum analysis tool.

Figure 2.

PSD of 3D heart motion data with frequency notations

Figure 3.

Heart motion bi-spectrum plot

The interactions, both amplitude and phase, can be detected by analysing the pre-recorded heart motion data using the bi-spectrum tool which quantifies relationships among the underlying oscillatory components of the observed signals [28,29]. Specifically, the bi-spectral analysis examines the relationships between the oscillations at two basic frequencies, fr(respiration mode frequency), fh(heartbeat mode frequency) and 2fr, 2fh and a modulation component at the frequency fr ± fh and etc. [30]

Figure 3 presents a typical bi-spectrum for the whole frequency domain for heart motion signal. The high peaks located at the bi-frequencies 0.16 Hz, 0.16 Hz correspond to the respiratory self-coupling and 1.82Hz, 1.82Hz correspond to the heartbeat self-coupling. The two sets of peaks located around 2Hz, 0.33 Hz and 0.33Hz, 2Hz are attributable to cardio-respiratory coupling between frequency component fr and fh. Several peaks around 2.2Hz, 2.2Hz are assumed to be coupling between the twice respiratory component 2fr and the difference fr+/−fh. The distinguished peaks and their bi-frequencies are collated in Table 1. Due to the symmetry of the bi-spectrum, only half of the main nonlinear coupling bi-frequency components are shown:

Table 1.

Bi-frequencies of beating heart motion

Peak	Nonlinear coupling bi-frequencies
1	(fr,fr)
2	(fr,fh)
3	(fr,fh+fr)
4	(fr,fh+2fr)
5	(2fr,fh)
6	(2fr,fh+fr)
7	(2fr,fh+2fr)
8	(fh,fh)
9	(fh,fh+4fr)
10	(fh+fr,fh+fr)
11	(fh+2fr,fh+2fr)
12	(fh+2fr,fh+4fr)
13	(fh+4fr,fh+4fr)

Two results from the bi-spectrum analysis of the heart motion data are derived:

The adaptive scheme is necessary. From the bi-spectrum in Figure 3 above, we can observe that the cardiac frequency span is 1.75Hz to 1.85Hz. Consequently, the coupling component has a wide frequency range resulting from the variation of the heart beating frequency. The width of the peaks in Figure 3 (b) also indicates that the instantaneous dominant mode frequencies of the heart motion are different from period to period.

The quadratic term is needed for the estimation model. In the presence of peaks at the nonlinear coupling bi-frequencies, we found that the typical nonlinear quadratic coupling regulates the dynamics of the beating heart. The square of sinusoid could appropriately approximate quadratic nonlinearity, which is the cardiac-respiration interaction. This square of sinusoid can be viewed as the modulated cardiac motion by respiration. The cardiac motion and respiration motion are modelled by truncated Fourier series respectively. The quadratic term is the product of respiration and heartbeat sinusoids, demonstrating the quadratic phase-coupling phenomena.

2.2 Nonlinear heart motion future reference prediction algorithm

The robot system requires a precise prediction of the immediate future of the heart motion in order to track the beating heart accurately. Figure 4 shows us the prediction scheme of immediate future heart motion using an adaptive modulated sinusoid model. The d[n], d̂[n]and x[n] denote the desired response, the predicted response and the input signal, respectively. The prediction problem in nature is to approximate the beating heart motion by the modulated sinusoid model with the historical heart motion measurements. The quadratic sinusoid model is updated recursively whenever a new piece of data comes in. The parameterization of the coefficient vector is achieved through EKF in a nonlinear way.

Figure 4.

Prediction scheme of the beating heart motion

2.2.1 Quadratic sinusoid nonlinear model

The heart motion can be modelled as the superposition of three parts: cardiac motion Xc, respiration motion Xr and interaction motion part Xm. The model can be represented as:

X (t) = X_{c} (t) + X_{r} (t) + X_{m} (t)

(2)

In which

X_{c} (t) = \sum_{k = 1}^{N_{c}} a_{k}^{c} \cos (2 π f_{k}^{c} t + ϕ_{k}^{c})

(3)

X_{r} (t) = \sum_{l = 1}^{N_{r}} a_{l}^{r} \cos (2 π f_{l}^{r} t + ϕ_{l}^{r})

(4)

X_{m} (t) = X_{c} \cdot X_{r} = \sum_{k = 1}^{N_{c 1}} \sum_{l = 1}^{N_{r 1}} a_{k l}^{m} \cos (2 π f_{k l}^{m} t + ϕ_{k l}^{m})

(5)

In (3–5) generally a denotes the magnitude, f denotes the frequency component and φ denotes the phase. (3) and (4) are the linear combinations of the significant cardiac harmonics and respiration harmonics respectively. (5) is the quadratic term meaning modulating the cardiac signal by respiration signal. Specifically, in (3) and (4), N and N are the number of harmonics in cardiac motion and respiration motion. respectively. In (5), $f_{k l}^{m} = f_{k}^{c} + f_{l}^{r}$ and $φ_{k l}^{m} = φ_{k}^{c} + φ_{l}^{r}$ describe the coupling frequency and coupling phase in beating heart motion. The coupling number of harmonics N_c1 and N_r₁are chosen according to selected frequencies in Figure 2 and Table 1.

Considering time-varying property of the heart motion, the time invariant model combined of three part Fourier series with a constant offset could be modified to

X_{n} (t) = c (t) + X (a (t), f (t), θ (t))

(6)

in which all the magnitude and frequency in (2) are time-varying and $θ_{i} (t) = i 2 π \int_{0}^{t} f (τ) d τ + φ_{i} (t)$

The state space model for this system is

\begin{array}{l} x (t + Δ t) = A (Δ t) x (t) + μ (t) \\ z (t) = h (x (t)) + v (t) \end{array}

(7)

with state vector defined as

\begin{array}{l} x (t) = [c (t), a^{c} ​_{k} (t), a^{r} ​_{l} (t), a^{m} ​_{k l} (t), f^{c} (t), f^{r} (t), \\ f^{m} (t), θ_{k}^{c} (t), θ_{l}^{r} (t), θ_{k l}^{m} (t)]^{T} \end{array}

Here the transfer matrixA could be derived if the variable c(t) and all a(t),f(t),φ(t) could evolve through a random walk.

A (Δ t) = [\begin{matrix} 1 & ​ & ​ & ​ \\ ​ & I_{N_{c} + N_{r} + N_{r 1}} & ​ & ​ \\ ​ & ​ & A_{f} (Δ t) & ​ \\ ​ & ​ & \begin{matrix} k Δ t \\ ⋮ \\ l Δ t \\ ⋮ \\ N_{r 1} Δ t \end{matrix} & A_{θ} (Δ t) \end{matrix}]

(8)

where A_f(Δt) and A_θ(Δt)hold the same structure to represent the coupling feature of the heart motion.

A_{f} (Δ t) = [\begin{matrix} I_{N c} & ​ & ​ \\ ​ & I_{N r} & ​ \\ \begin{matrix} 1 & ​ \end{matrix} & \begin{matrix} ​ & 1 & ​ & ​ \end{matrix} & ​ \\ \begin{matrix} 1 & ​ \end{matrix} & \begin{matrix} ​ & ​ & ​ & 1 \end{matrix} & ​ \\ \begin{matrix} ​ & 1 \end{matrix} & \begin{matrix} ​ & 1 & ​ & ​ \end{matrix} & 0_{N_{c 1}} \\ \begin{matrix} ​ & 1 \end{matrix} & \begin{matrix} ​ & ​ & ​ & 1 \end{matrix} & ​ \end{matrix}]

(9)

Additionally, the h(x(t)) could be formulated as

\begin{array}{l} H^{T} = {(\frac{δ h}{δ x})}^{T} |_{\hat{x} (t + Δ t) = A \hat{x} (t)} \\ = [\begin{matrix} 1 \\ \cos {\hat{θ}}_{1}^{c} (t + Δ t) \\ ⋮ \\ \cos {\hat{θ}}_{k l}^{m} (t + Δ t) \\ \begin{array}{l} ⋮ \\ ⋮ \end{array} \\ - {\hat{a}}_{k}^{c} (t + Δ t | t) \sin {\hat{θ}}_{1} (t + Δ t) \\ ⋮ \\ - {\hat{a}}_{k l}^{m} (t + Δ t | t) \sin {\hat{θ}}_{m} (t + Δ t) \end{matrix}] \end{array}

(10)

2.2.2 EKF parameters estimation

z(t)in the measurement equation in (7) is corrupted by measurement uncertainty variable v(t). The prediction task requires the estimation of the state variables. The EKF is used for parameter estimation since the measurement equation is nonlinear. The detailed derivation of the EKF for this model is formulized as [17]:

\begin{array}{l} P (t + Δ t | t) = A P (t | t) F^{T} + Q_{h} \\ S = σ_{R}^{2} + H P (t + Δ t | t) H^{T} \\ K = P (t + Δ t) H^{T} S^{- 1} \\ \hat{x} (t + Δ t | t + Δ t) = A \hat{x} (t | t) + K (z (t + Δ t) - h (A \hat{x} (t | t))) \\ P (t + Δ t | t + Δ t) = (I - K H) P (t + Δ t | t) \end{array}

(11)

3. Model Following Control Algorithm

Figure 5 illustrates the control mechanism diagram. With the adaptive nonlinear heart motion model in hand, the tracking control problem could be simply transformed into the heart motion model following problem. The proposed control algorithm has the advantage of prediction of future motion by using the adaptive nonlinear heart motion model. The control effort could be divided into the feedforward portion and the feedback portion to handle the high frequency band heart motion signal and lower frequency band respiration signal, respectively. Furthermore, the nonlinear heart motion model will improve the heart motion tracking performance by improving the motion prediction accuracy.

Figure 5.

Model following control architecture based on the heart motion model

The block diagram of the model following control algorithm is shown in Figure 6. Considering the nonlinear nature of heart motion, EKF is used to estimate the states of nonlinear heart motion. Note that with the help of the robot experimental model, another EKF is used to provide the integrated robot state estimation feedback to the optimal tracking controller to form the closed-loop, which is different from the one used in parameter estimation in the heart motion model.

Figure 6.

block diagram of model-following control algorithm Optimal tracking controller phase:

We stress that the feedforward controller does include the anticipatory characteristics for the heart motion reference. This is a time-varying dynamic feedforward controller because it embedded the adaptive model of the heart motion. The nonlinear model poses the time-varying property to give the one step ahead future reference to the controller for the model following task. The controller could adjust the feedforward gain to obtain the improved transit response of the robotics. The controller design procedure could be divided into two phases: optimal tracking controller design and state estimator design.

The model following method combines the heart motion states with the robot state to form a whole new state space model. The solution will be obtained by using the Hamilton-Jacobi method [31]. The formulation is as follows [21]:

The heart motion model is rewritten in terms of the state:

\begin{array}{l} z [k + 1] = f (z [k]) + w_{h} \\ \tilde{y} = z_{k} + v_{h} \end{array}

(12)

The robotics dynamics model is given:

\begin{array}{l} x [k + 1] = F x [k] + G u [k] + w_{r} \\ y [k] = H x [k] + v_{r} \end{array}

(13)

where F, G, H are system matrix, input matrix and output matrix, respectively. x[k] is the robot state vector. u[k] is the control effort. The whole new system state space model is given:

\begin{array}{l} \bar{X} [k + 1] = \hat{F} \bar{X} [k] + \hat{G} u [k] \\ \bar{y} [k] = \hat{H} \bar{x} [k] + v \end{array}

(14)

where:

\bar{x} = [\begin{array}{l} x \\ z \end{array}]

(15)

\hat{F} = (\begin{matrix} F & 0 \\ 0 & \partial f ∕ \partial t \end{matrix})

(16)

\hat{G} = [\begin{array}{l} G \\ 0 \end{array}]

(17)

\hat{H} = [\begin{array}{l} H \\ 1 \end{array}]

(18)

The cost function of the optimal tracking represents as:

\begin{array}{l} J = \sum_{k = k_{0}}^{k_{0} + T} ({(\bar{x} [k] - \tilde{x} [k])}^{T} Q (\bar{x} [k] - \tilde{x} [k]) \\ + u^{T} [k] R u [k]) \end{array}

(19)

where Q is the non-negative definite symmetric matrix andR is positive symmetric definite matrix. x˜ is the ideal state. y˜ is the ideal output. Their relationship is

\tilde{x} = L \tilde{y}

(20)

L = {\hat{H}}^{T} {(\hat{H} {\hat{H}}^{T})}^{- 1}

(21)

Q = {(I - L \hat{H})}^{T} Q_{1} (I - L \hat{H}) + {\hat{H}}^{T} Q_{2} \hat{H}

(22)

Q₁ and Q₂ are both negative definite symmetric matrixes. Thus the cost function could be reduced to:

J = \sum_{k = k_{0}}^{k_{0} + T} ({\bar{x}}^{T} [k] \hat{Q} \bar{x} [k] + u^{T} [k] R u [k])

(23)

where

\hat{Q} = (\begin{matrix} Q & - Q L {\hat{H}}^{T} \\ - \hat{H} L^{T} Q & \hat{H} L^{T} Q L {\hat{H}}^{T} \end{matrix})

(24)

Control law:

\begin{aligned} u [k] = - ({\hat{G}}^{T} S [k + 1] \hat{G} + R)^{- 1} {\hat{G}}^{T} (S [k + 1] \hat{F} \bar{x} [k] \\ M [k + 1]) \end{aligned}

(25)

where S and M could be given in a recursive way.

\begin{array}{l} S [k] = \hat{F} (S [k + 1] - S [k + 1] \hat{G} ({\hat{G}}^{T} S [k + 1] \hat{G} \\ + R)^{- 1} {\hat{G}}^{T} S [k + 1] \hat{F}) + \hat{Q}) \end{array}

(26)

M [k] = (\hat{F} + K^{T} [k] {\hat{G}}^{T}) M [k + 1] - \hat{Q} L \tilde{y}

(27)

K is the time-varying feedback gain:

K [k] = - {({\hat{G}}^{T} S [k + 1] \hat{G} + R)}^{- 1} {\hat{G}}^{T} S [k + 1] \hat{F}

(28)

The resulting control law is composed of feedback and feedforward parts, which are:

u [k] = u_{f} [k] + u_{b} [k]

(29)

u_{b} [k] = - {({\hat{G}}^{T} S [k + 1] \hat{G} + R)}^{- 1} {\hat{G}}^{T} S [k + 1] \hat{F} \bar{x} [k]

(30)

u_{f} [k] = - {({\hat{G}}^{T} S [k + 1] \hat{G} + R)}^{- 1} \hat{G} M [k + 1]

(31)

State estimator phase design:

The state space model for this system comes from equations (14,15). The process error and measurement error could be modelled as zero mean and corresponding covariance:

\begin{array}{l} w ~ (0, Q_{k}) \\ v ~ (0, R_{k}) \end{array}

(32)

The EKF is used to state estimation since the state equation is nonlinear. The time update of estimated error and estimation error covariance could be denoted as:

\begin{array}{l} P_{k}^{-} = {\hat{F}}_{k - 1} P_{k - 1}^{+} {\hat{F}}_{k - 1}^{T} + Q_{k} \\ {\hat{\bar{x}}}_{k}^{-} = {[F {\hat{\bar{x}}}_{k - 1}, f_{k - 1}]}^{T} \end{array}

(33)

The measurement update of the state estimate and estimation error covariance could be denoted as:

\begin{array}{l} K_{k} = P_{k}^{-} {\hat{H}}_{k} ​^{T} {(R_{k} + {\hat{H}}_{k} P_{k}^{-} {\hat{H}}_{k} ​^{T})}^{- 1} \\ {\hat{\bar{x}}}_{k}^{+} = \hat{\bar{x}} ​_{-}^{k} + K_{k} [\bar{y} - {\hat{H}}_{k} \hat{\bar{x}} +_{k}^{-} +_{]} \\ P_{k}^{+} = (I - K_{k} {\hat{H}}_{k}) P_{k}^{-} \end{array}

(34)

4. Experiment and Result Discussion

4.1 Experiment setup

The control algorithm mentioned above is implemented on the 3D robot test-bed named PHANToM Premium 1.5A, which has good characteristics to be able to perform the medical surgery. The general dynamics of PHANToM could be denoted as:

M (θ) \dot{θ} + C (\dot{θ}, θ) \dot{θ} + N (θ) = τ

(35)

where ⊝ is angular vector for three axes, M is rotation inertia matrix, C is Coriolis matrix and N is other forces including gravity, friction and so on. The Coulomb friction model describes the friction in the robotics. After the compensation of gravity and friction in the robotics system, C is calculated small enough such that the nonlinear of the system is overcome by the linear one. That is why we used the linear experimental model for the control algorithm. In the experimental setup, the control algorithms were executed on an Intel E5-2609 2.4GHz quad cores PC running MATLAB xPC Target v2.8 real-time kernel with a sampling time of 0.5ms. The experiment utilized pre-recorded real heart motion data rather than an on-line measurement of the target motion.

4.2 Experiments and results

In this section, we evaluated three types of tracking controller to be tested by two different pre-recorded datasets presented in Section 2: model following controller with exact reference (MF_R), model following controller with linear sinusoid series prediction (MF_LSS) (see Appendix A) and model following controller with quadratic nonlinear prediction (MF_QN). For comparison purposes, the MF_R is a non-causal method which gives the baseline performance of the model following control architecture in order to provide the standard for algorithm comparison.

The tracking results are derived through fine-tuning the parameters to minimize the tracking errors for the algorithms. The parameters for each approach on the two datasets are summarized in Table 2. The duration of the experiment is 60 seconds for both datasets. The root mean square of the tracking errors and control efforts are summarized in Table 3. Figure 7 shows the tracking results from the first dataset and Figure 8 shows the tracking results from the second dataset.

Table 2

Parameter comparison

Control Algorithm	Parameters	1^st dataset	2^nd dataset
MF_LSS	Number of harmonies Number of states	6 14	10 22
MF_QN	Number of harmonies Number of states	10 22	16 34

Table 3.

Experiment results

Tracking results	RMS position error[mm] (Max position error[mm])		RMS control effort[Nmm]
Tracking results	1^st dataset	2^nd dataset	1^st dataset	2^nd dataset
MF_R	0.251 (1.355)	0.218 (1.223)	51.34	48.15
MF_LSS	0.602 (1.497)	0.487 (1.427)	69.31	63.34
MF_QN	0.298 (1.386)	0.275 (1.263)	97.57	86.66

We found that the MF_QN algorithm outperforms the MF_LSS algorithm in terms of tracking accuracy. The results prove that by using the nonlinear heart motion prediction, the controller is more responsive to the sudden changes in heart motion. The RMS error and the worst case error are largely decreased. Comparing the results of MF_QN with that of MF_R reveals that the nonlinear prediction algorithm is very close to the exact reference. When the nonlinearity of the phase-coupling in the heart motion is relatively high, which is the situation in dataset 2, the tracking performance by MF_QN is not only improved, but also more robust. We will now discuss the experiment results through the two main properties of the model following controller:

Feedforward control effort. MF_QN poses the feedforward control path to compensate the time lag of the feedback controller only. The time-varying feedforward gain derives from the future reference predicted from the adaptive time-varying nonlinear model. MF_LSS is the same as MF_QN that has the feedforward control path and time-varying gain. This is the property which is a necessity in robot tracking of beating heart surgery.

Embedded model.

Time-varying heart rate representation. Both methods have the ability to compensate for the heart motion changes through different ways in some extent. The linear sinusoid series model varies its filter parameters to fit new blocks of data if the heart motion frequency varies. The nonlinear model changes both its linear part parameters and it quadratic part parameters to take into account a change in heart motion frequency.

Phase-coupling representation. The MF_QN is better than MF_LSS because of the coupling term between heartbeats and respiration, which are the two components of the heart motion. Linear sinusoid series model is just a linear combination of the characteristics' frequency element in heartbeats and respiration, and is not able to adapt itself to the newly produced frequency components, whereas the nonlinear model has the quadratic term to represent the frequency and phase-coupling. Figure 7 and Figure 8 (b)(c) show this results clearly. However, when the heart motion nonlinearity is not great in a certain axis, as shown in Figure8(a), the MF_QN over fits a little and is not as good as MF_LSS.

Figure 7.(a).

Heart motion model following algorithm results in X axis using the first dataset

Figure 7.(b).

Heart motion model following algorithm results in Y axis using the first dataset

Figure 7.(c).

Heart motion model following algorithm results in Z axis using the first dataset

Figure 8.(a).

Heart motion model following algorithm results in X axis using the second dataset

Figure 8.(b).

Heart motion model following algorithm results in Y axis using the second dataset

Figure 8.(c).

Heart motion model following algorithm results in Z axis using the second dataset

5. Conclusions

In this paper, we aim to improve the tracking performance by a nonlinear prediction-based model following controller. In this study, a bi-spectral method is used to analyse the phase-coupling nature between the respiration component and the heartbeat component to show physically the nonlinearity phenomenon dominating the beating heart motion. A novel nonlinear adaptive prediction model is incorporated into the model following controller. In this method, the quadratic sinusoid representing the cardiac-respiration modulation is taken into account in the model, such that the heart motion dynamics is described by more detail in the frequency domain. A comparison of the experiment results is provided to evaluate the proposed algorithm using two different vivo datasets. The MF_QN algorithm shows convincing results and properties. With regard to feedforward control, MF_QN has the time-varying control gain to cancel the fast movement caused by the heartbeat. Looking at nonlinear heart motion characteristics, MF_QN has the quadratic part in its structure to model the phase-coupling effect and changes the coefficients of the quadratic term to adapt the heart motion. To conclude, the MF_QN has achieved a significant tracking accuracy improvement because the embedded heart motion model could adaptively describe the irregular heart motion precisely and quickly.

In the future, the effect of different POI positions on the prediction algorithm will be investigated. Additionally, the multistep receding horizon prediction via the nonlinear heart motion model will be given more attention because it will provide more information and be more adjustable for the controller to output better control efforts in order to achieve better tracking performance.

6. Appendix A

Time-varying Fourier series with EKF

The time-varying model is

y (t) = c + \sum_{i = 1}^{m} r_{i} (t) \sin θ_{i} (t)

(A.1)

Where $θ_{i} (t) = i \int_{0}^{t} ω (τ) d τ + φ_{i} (t)$ , ω is the heart rate, r_i and φ_i are the amplitude and phase of the ith harmony. The state space model for this system is

\begin{array}{l} x (t + Δ t) = A (Δ t) x (t) + μ (t) \\ z (t) = h (x (t)) + v (t) \end{array}

(A.2)

with state vector defined as $x (t) = [c (t), r_{i} (t), w (t), θ_{i} (t)]^{T}$ . Here the transfer matrix A could be derived as (A.3), if the variable c(t),r_i(t),w(t),φ_i(t) could evolve through a random walk. Additionally, the linearized h(x(t)) could be formulated (A.4) as the EKF is used to parameter estimation since the measurement equation is nonlinear. The detailed derivation of the EKF for this model could be formulized as [17,32].

A (Δ t) = [\begin{matrix} I_{m + 1} & ​ & ​ & ​ & 0 \\ ​ & 1 & ​ & ​ & ​ \\ 0 & Δ t & ⋱ & ​ & ​ \\ ​ & ⋮ & ​ & 1 & ​ \\ ​ & m Δ t & ​ & ​ & 1 \end{matrix}]

(A.3)

\begin{array}{l} H^{T} = (\frac{δ h}{δ x})^{T} |_{\hat{x} (t + Δ t) = A \hat{x} (t)} \\ = [\begin{matrix} 1 \\ \sin {\hat{θ}}_{1} (t + Δ t) \\ ⋮ \\ \sin {\hat{θ}}_{1} (t + Δ t) \\ 0 \\ {\hat{r}}_{1} (t + Δ t | t) \cos {\hat{θ}}_{1} (t + Δ t) \\ ⋮ \\ {\hat{r}}_{m} (t + Δ t | t) \cos {\hat{θ}}_{m} (t + Δ t) \end{matrix}] \end{array}

(A.4)

\begin{array}{l} P (t + Δ t | t) = A P (t | t) F^{T} + Q \\ S = σ_{R}^{2} + H P (t + Δ t | t) H^{T} \\ K = P (t + Δ t) H^{T} S^{- 1} \\ \hat{x} (t + Δ t | t + Δ t) = A \hat{x} (t | t) + K (z (t + Δ t) - h (A \hat{x} (t | t))) \\ P (t + Δ t | t + Δ t) = (I - K H) P (t + Δ t | t) \end{array}

(A.5)

Footnotes

7. Acknowledgments

The study is supported by the China National Science Foundation, Group Creative Project (number 61121003) issued to the Science Technology on Inertia Laboratory, Beihang University, Beijing, China.

Thanks to Dr. Cavosoglu at Case Western Reserve University, Cleveland, OH, U.S., for the help and advice related to the research in this paper given during Fan Liang's visiting period in the U.S.

References

Park

Cha

and Park

(2010). Development of Biomedical Mcrorobot for Intravascular Therapy, International Journal of Advanced Robotic Systems, Kordic

Vedran

Lazinica

Aleksandar

Merdan

Munir

(Ed.), ISBN: 1729-8806, InTech, Available from: http://www.intechopen.com/journals/international_journal_of_advanced_robotic_systems/development-of-biomedical-microrobot-for-intravascular-therapy

Newman

M. F.

Kirchner

J. L.

Phillips-Bute

Gaver

Grocott

Jones

R. H.

(2001)., Longitudinal assessment of neurocognitive function after coronary artery bypass surgery New England Journal of Medicine, vol. 344, no. 6, pp. 395–402

Puskas

J. D.

Wright

C. E.

Ronson

R. S.

Brown

W. M.

Gott

J. P.

and Guyton

R. A.

(1998), Off-pump multi-vessel coronary bypass via sternotomy is safe and effective Annals of Thoracic Surgery, vol. 66, no. 3, pp. 1068–1072

Trejos

Salcudean

Sassani

and Lichtenstein

(1999). On the feasibility of a moving support for surgery on the beating heart. In Medical Image Computing and Computer-Assisted Intervention–MICCAI'99, 1088–1097

Bebek

, and Çavusogavusoglu

M.C.

(2007). Intelligent control algorithms for robotic-assisted beating heart surgery. Robotics, IEEE Transactions on 23 (3): 468–480.

Cavusoglu

M.C.

Rotella

Newman

W.S.

Choi

Ustin

and Sastry

S.S.

(2005). Control algorithms for active relative motion cancelling for robotic assisted off-pump coronary artery bypass graft surgery. In Advanced Robotics, 2005. ICAR'05. Proceedings., 12th International Conference on, 431–436.

Rotella

(2005). Predictive tracking of quasi periodic signals for active relative motion cancellation in robotic assisted coronary artery bypass graft surgery. M.S. Thesis, Case Western Reserve University, Cleveland, OH, USA, August 2004. [Online]. Available:http://robotics.case.edu/publications.html#thesis

Patel

H.R.H.

Amodeo

and Joseph

J.V.

(2009). Robotic Oncological Surgery: Technology That's Here to Stay?, International Journal of Advanced Robotic Systems, Kordic

Vedran

Lazinica

Aleksandar

Merdan

Munir

(Ed.), ISBN: 1729-8806, InTech, Available from: http://www.intechopen.com/journals/international_journal_of_advanced_robotic_systems/robotic-oncological-surgery-technology-that-s-here-to-stay-

Nakamura

Kishi

and Kawakami

(2001) Heartbeat synchronization for robotic cardiac surgery, in Proc. of IEEE International Conference on Robotics and Automation (ICRA), vol. 2, Seoul, Korea, May 2001, pp. 2014–2019.

10.

Ortmaier

(2003). Motion compensation in minimally invasive robotic surgery. Technische Universität München, Universitätsbibliothek

11.

Jin

Chang

P. Hun

and Choi

(2011). High-Accuracy Tracking Control of Robot Manipulators Using Time Delay Estimation and Terminal Sliding Mode, International Journal of Advanced Robotic Systems, Kordic

Vedran

Lazinica

Aleksandar

Merdan

Munir

(Ed.), ISBN: 1729-8806, InTech, Available from: http://www.intechopen.com/journals/international_journal_of_advanced_robotic_systems/high_accuracy_tracking_control_of_robot_manipulators_using_time_delay_estimation_and_terminal_slidin

12.

Franke

T.J.

Bebek

and Cavusoglu

M.C.

(2007). Improved prediction of heart motion using an adaptive filter for robot assisted beating heart surgery. In Intelligent Robots and Systems, 2007. IROS 2007. IEEE/RSJ International Conference on, 509–515.

13.

Franke

T.J.

Bebek

and Cavusoglu

M.C.

(2008). Prediction of heartbeat motion with a generalized adaptive filter. In Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on, 2916–2921.

14.

Thakral

Wallace

Tomlin

Seth

and Thakor

(2001). Surgical motion adaptive robotic technology (SMART): Taking the motion out of physiological motion. In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2001, 317–325.

15.

Ginhoux

Gangloff

De Mathelin

Soler

Sanchez

M.M.A.

and Marescaux

(2005). Active filtering of physiological motion in robotized surgery using predictive control. Robotics, IEEE Transactions on 21 (1): 67–79.

16.

Richa

Bó

A.P.L.

and Poignet

(2010). Beating heart motion prediction for robust visual tracking. In Robotics and Automation (ICRA), 2010 IEEE International Conference on, 4579–4584

17.

Yuen

S.G.

Kettler

D.T.

Novotny

P.M.

Plowes

R.D.

and Howe

R.D.

(2009). Robotic motion compensation for beating heart intracardiac surgery. The International journal of robotics research 28 (10): 1355.

18.

Yuen

S.G.

Perrin

D.P.

Vasilyev

N.V.

del Nido

P.J.

and Howe

R.D.

(2010). Force tracking with feed-forward motion estimation for beating heart surgery. Robotics, IEEE Transactions on 26 (5): 888–896.

19.

Bachta

Renaud

Cuvillon

Laroche

Forgione

and Gangloff

(2009). Motion prediction for computer-assisted beating heart surgery. Biomedical Engineering, IEEE Transactions on 56 (11): 2551–2563.

20.

Shiogai

Stefanovska

and McClintock

P.V.E.

(2010). Nonlinear dynamics of cardiovascular ageing. Physics reports 488 (2–3): 51–110.

21.

Liang

Meng

Dong

(2012) Model following control of assisted robotics of CABG surgery, Optics and Precision Engineering, 2012 Vol. 20 (1): 131–139 (in Chinese)

22.

Servant

Logier

Mouster

and Goudemand

(2009). [Heart rate variability. Applications in psychiatry]. L'Encéphale 35 (5) (October): 423–428. doi:10.1016/j.encep.2008.06.016.

23.

Tuna

E. Erdem

. (2011). Heart motion prediction based on adaptive estimation algorithms for robotic-assisted beating heart surgery. M.S. Thesis, Case Western Reserve University, Cleveland, OH, USA, [Online]. Available: http://engr.case.edu/cavusoglu_cenk/Students.html

24.

Priestley

M.B.

(1988). Non-linear and non-stationary time series analysis. Academic Press, London, 1988.

25.

Schäfer

Rosenblum

M.G.

Kurths

, H.H. Abel and others. (1998). Heartbeat synchronized with ventilation. Nature 392 (6673): 239.

26.

Lotri

M.B.

and Stefanovska

2000. Synchronization and modulation in the human cardiorespiratory system. Physica A: Statistical Mechanics and its Applications 283 (3–4): 451–461.

27.

Janson

N.B.

Balanov

A.G.

Anishchenko

V.S.

and McClintock

P.V.E.

(2001). Phase synchronization between several interacting processes from univariate data. Physical Review Letters 86 (9): 1749–1752.

28.

Jamšek

Stefanovska

McClintock

P.V.E.

and Khovanov

I.A.

(2003). Time-phase bispectral analysis. Physical Review E 68 (1): 016201.

29.

Fackrell

J.W.A.

(1996). Bispectral analysis of speech signals. PhD Thesis University of Edinburgh

30.

Jamšek

Stefanovska

and McClintock

P.V.E.

(2004). Nonlinear cardio-respiratory interactions revealed by time-phase bispectral analysis. Physics in medicine and biology 49: 4407

31.

Anderson

B. D.

and Moore

J. B.

, Optimal Control Linear Quadratic Methods[M]. Englewood Cliffs, NJ: Prentice-Hall,1990.

32.

Ali

(2008) Sayed Adaptive Filters, A John Wiley & Sons, Inc., publication

A Quadratic Nonlinear Prediction-Based Heart Motion Model following Control Algorithm in Robotic-Assisted Beating Heart Surgery

Abstract

Keywords

1. Introduction

1.1 Beating heart-assisted robotics system - introduction and requirements

1.2 Related Work

2. Heart Motion Nonlinear Prediction

2.1 Nonlinear dynamics and bi-spectral analysis

2.2 Nonlinear heart motion future reference prediction algorithm

2.2.1 Quadratic sinusoid nonlinear model

2.2.2 EKF parameters estimation

3. Model Following Control Algorithm

4. Experiment and Result Discussion

4.1 Experiment setup

4.2 Experiments and results

5. Conclusions

6. Appendix A

Footnotes

7. Acknowledgments

References