Heart Motion Prediction in Robotic-Assisted Beating Heart Surgery: A Nonlinear Fast Adaptive Approach

Abstract

Off-pump Coronary Artery Bypass Graft (CABG) surgery outperforms traditional on-pump surgery because the assisted robotic tools can alleviate the relative motion between the beating heart and robotic tools. Therefore, it is possible for the surgeon to operate on the beating heart and thus lessens post surgery complications for the patients. Due to the highly irregular and non-stationary nature of heart motion, it is critical that the beating heart motion is predicted in the model-based track control procedures. It is technically preferable to model heart motion in a nonlinear way because the characteristic analysis of 3D heart motion data through Bi-spectral analysis and Fourier methods demonstrates the involved nonlinearity of heart motion. We propose an adaptive nonlinear heart motion model based on the Volterra Series in this paper. We also design a fast lattice structure to achieve computational-efficiency for real-time online predictions. We argue that the quadratic term of the Volterra Series can improve the prediction accuracy by covering sharp change points and including the motion with sufficient detail. The experiment results indicate that the adaptive nonlinear heart motion prediction algorithm outperforms the autoregressive (AR) and the time-varying Fourier-series models in terms of the root mean square of the prediction error and the prediction error in extreme cases.

Keywords

Heart Motion Prediction Volterra Model Quadratic Phase Coupling Beating Heart Surgery

1. Introduction

Beating heart surgery avoids cognitive loss [1] caused by the application of a cardio-pulmonary bypass machine and reduces time spent hospitalized and cost [2]. However, it is more demanding than traditional methods because beating heart surgery is unable to be carried out manually due to the quick motions of the heart [3]. The advancement of robotic technology provides an effective alternative to resolve this challenge by eliminating the relative motion between the beating heart and robotic tools to achieve a relatively stationary surgical environment.

1.1 Beating Heart Assisted Robotics System and Heart Motion Prediction

As an assisted surgical robotic tool, a beating heart assisted robotics system was designed to allow surgeons to perform operations as if the heart was stationary. The relatively stationary motion is realized through synchronizing the beating heart with surgical equipment, thus cancelling the relative motion between the Point Of Interest (POI) on the beating heart surface and the end-effecter of the robot. A conceptual model of the robotic tele-operational beating heart surgery system is shown in Figure 1. 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 faces the still heart view on the screen and manipulates the joysticks to conduct surgical operations. With the assistance of the robotic tool, surgeons can operate on a heart surface with high-frequency motions.

Figure 1.

Robotic tele-operational system concept model[4]

The challenges in identifying the solution to the beating heart motion-tracking problem are two-fold: sufficient accuracy and fast tracking. On the one hand, in off-pump CABG surgery, surgeons can operate on blood vessels ranging in diameter from 0.5–2.0mm. The dynamic position error has to be controlled within 1% of the blood vessel diameter or in the order of 100–250μm root mean square in order to operate safely [5]. On the other hand, because the motion of the POI on the heart surface has a relatively high bandwidth and associates with non-stationary characteristics, the robot tools need to predict and follow the POI in a very short period. A traditional causal feedback controller alone is not sufficient for beating heart surgery, since it cannot follow the moderate to high frequency component of the heart motion. A predictive controller with a receding horizon prediction of the heart motion in the feed forward path is necessary [5]. The heart motion-tracking problem could be reduced into a heart motion prediction problem if a satisfactory prediction could be provided [6]. In this paper, we focus on alleviating heart motion prediction error.

1.2 Related Work

Literature discussing heart motion modelling and prediction can be roughly categorized into three groups: model free methods, time domain methods and frequency domain methods.

Ortmaier in [7] proposed a global robust prediction scheme in which the algorithm reconstructs the underlying dynamics of the heart motion through past measurements. This is a model-free method based on Taken Theory, which can predict the heart motion behaviour even though the visual sensor is obscured for a period of time.

Bebek and Cavusoglu in [4] presented a time domain method to predict the current circle of heart motion based on the previous circle of heart motion, with the help of an Electrocardiogram (ECG) signal to correct the heart motion period. Frank et al. proposed a new Autoregressive (AR) model based adaptive prediction technique [8, 9] that is less susceptible to the variation of the heart motion period and robust to noise. However, this AR model is still a linear model that cannot describe some of the nonlinear dynamics of heart motion. [10] introduces a method using recursive Bayesian estimation based on multiple motion models that represent the motion of the target organ during cycles of different amplitudes. Nonlinearity is overcome by interpolating between the motion estimates of these models based on the amplitude of the current cycle. However, the models are just for respiration only.

For the frequency domain method, researchers separate the heart motion into cardiac and respiration [11, 12] and describe them as the summation of the truncated Fourier series. Richa et al. [13] model the quasi-periodic beating heart motion as a time-varying dual Fourier series and use an Extended Kalman Filter (EKF) to estimate the parameters. Likewise in [14, 15] a truncated time varying Fourier series with an offset plus EKF is used to model the motion and estimate the characteristic parameters recursively.

The limitations of current methods include: 1) A lack of methods to investigate and analyse the nonlinearity in the dynamics of beating heart motion. 2) A lack of methods to mathematically model the nonlinearity. The deficiency of explicit modelling of the nonlinear part and the coupling between cardiac motion and respiration motion in the beating heart motion would lead to insufficient mathematical descriptions of the dynamics of beating heart motion and thus the increase of prediction errors. Although Bachta et al. [16] address the coupling of cardiac motion and respiration motion and model the coupling as amplitude modulation between two period signals, nonlinear analysis of the beating heart motion is still needed to identify the cause of the nonlinearity [17].

The advantages of the method proposed in this paper include: 1) In comparison to the model free method, the proposed method explicitly models the dynamics of the beating heart motion using the Volterra Series. 2) In comparison to the model based methods, the proposed method is extended by the addition of describing the nonlinear motion mode. It is the newly produced frequency components due to nonlinearity that make the beating heart motion hard to predict and follow using a robot tool. Therefore, it is a good investment to explore nonlinear modelling of the beating heart motion.

Volterra Series expansion is widely used in physiology signal estimation and modelling [18 –23]. For example, Mitsis et al. [19] present the insulin–glucose relationship by using the Volterra Model. Chen at al. [23] use the Volterra Series to model the instantaneous mean of the heartbeat interval (R-R interval) extracted using ECG. It shall be pointed out that the second order nonlinear term in the Volterra Series in our paper is employed to describe the important physical phenomenon of the coupling between cardiac motion and respiration motion in the beating heart motion, which is the mechanic motion on the surface of the heart and is different from the electrical related signal mentioned in [23]. To the best of our knowledge, the Volterra Series model was first introduced to characterize the beating heart motion for robot-assisted surgery. Attempting to improve tracking performance based on the second order Volterra Series model and extended time domain method AR model in [8, 9] by adding the nonlinear quadratic term, we propose a novel nonlinear prediction method in this paper. The main reasons for choosing this approach are multifaceted: 1) this nonlinear prediction method is naturally derived from the nonlinear dynamics of the heart motion. 2) The nonlinear term in the Volterra Series is the most appropriate for a good estimation of the phase-coupling phenomenon in heart motion dynamics. The parallel structure of the proposed algorithm allows the separation of linear and nonlinear elements of the process. To be more specific, the first order in the expansion is the linear part, which works identically to the AR filter to model basic motion in time as in previous studies [8, 9]. The second order is the quadratic part, which can model the coupling between cardiac motion and respiration motion [17]. 3) The computation complexity of this method can be reduced by applying a nonlinear Lattice structure from the fast recursive least square theory.

2. Heart Motion Nonlinear Analysis

We use Fourier analysis to obtain the prior knowledge of heart motion and use Bi-spectral analysis to reveal the nonlinear nature of heart motion.

2.1 Fourier analysis

We obtain the heart motion information through Fourier Analysis (FA) based on collected data [24]. In a 200s data collection period, the average heart rate of the animal model is 110 beats per minute. Figure 2 shows the Power Spectrum Density (PSD) of beating heart motion in 3 Cartesian axes. We observe two dominant modes, the first one is at 0.33Hz, which corresponds to the respiration motion and the second one is at 1.82Hz, which is the main mode of the whole beating heart and corresponds to the heartbeat. In the time domain, the peak displacement of POI is 11.5 mm from its mean location with a Root Mean Square (RMS) value of 5.3mm. The motion below 1Hz has a 2.67mm RMS value whereas the motion above 1Hz has a 4.27 mm RMS value. In the frequency domain, the bandwidth of the heart motion is up to 20Hz, since there are still some components with a PSD magnitude value of about 7E-6 at around 20Hz. Several harmonics that hold significant energy between 1Hz and 5Hz not only contain the integral times of the two dominant modes' frequency but also include the other mode frequency, which could be explained as the modulation effects among respiration, cardiac motion and other cardinal activities [17, 25]. Fourier analysis is not sufficient to interpret and represent this nonlinear nature since it could not clearly show the newly produced coupling frequencies caused by nonlinearity. Therefore another analysis tool, Bi-Spectrum, is introduced.

Figure 2.

PSD of the 3D heart motion data

2.2 Nonlinear dynamics and Bi-spectral Analysis

The dynamics of the beating heart can be described as the coupled oscillators mainly between the cardiac and respiratory oscillations [17] which are formulated as

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

(1)

Where F(P)is the common structure for the N oscillators (N=2 for this study) and 5F(Pi) is the derivative of the F(P_i), V_ij(P_i, P_j) is the interacting term.

It is the general nonlinear model for a beating heart, which represents the N dominant oscillators and their mutual modulation. We can refer to the beating heart motion as a process driven by coupled respiration motion and cardiac motion, each with individual characteristic frequency. Coupled oscillators interact with each other through their amplitudes and phases. Thus, the coupled phase phenomenon can alternatively be called phase coupling (PC). In PC, a new frequency component will be produced and its phase angles will be dependent on the phase angles of the original components [26]. Some researchers [27 –29] emphasise the beating heart synchronization properties of the phase dynamics in revealing the phase coupling nature. Quadratic phase coupling (coupling at sum and difference frequencies) occurs when a signal is passed through a square law device [26], which is the first motivation to model the heart motion by a quadratic term.

By observing Figure 2 again we find there are new frequency components other than dominate mode frequency components. Actually, the newly produced frequency component is in a certain format instead of being distributed randomly. In order to reveal the nonlinear phase coupling nature of the beating heart motion, we use the higher order statistics analysis tool bi-spectrum.

Both the amplitude and phase of the interactions can be detected by analysing the recorded time series using the bi-spectrum tool, which quantifies relationships among the underlying oscillatory components of the observed signals [30, 31]. Specifically, bi-spectral analysis examines the relationships between the oscillations at two basic frequencies, fr (respiration mode frequency), fh (heart beat mode frequency), 2fr, 2fh and a modulation component at the frequency fr ± fh [32].

Figure 3 presents a typical bi-spectrum for the whole frequency domain for the beating heart motion signal. The high peaks are located at the bi-frequency (0.16 Hz, 0.16 Hz) and (1.82Hz, 1.82Hz) belonging to the respiratory self-coupling and cardiac motion self-coupling, respectively. The two sets of peaks located around (2Hz, 0.33Hz) and (0.33Hz, 2Hz) are attributable to cardio-respiratory coupling between frequency components fr and fh. Moreover, there is also coupling between the respiratory component fr and the difference fr+/− fh. All of the newly produced bi-frequencies are due to a nonlinear coupling mechanism. The distinguished peaks and their bi-frequencies are collected in Table 1. Due to the symmetry of the bi-spectrum, only half of the main nonlinear coupling bi-frequencies components are shown.

Figure 3.

Bi-spectrum of beating heart motion

Table 1.

Bi-frequencies of beating heart motion

Peak	Nonlinear coupling Bifrequency
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)

We derive two results from the bi-spectrum analysis of the heart motion.

It is necessary to import the adaptive scheme in the estimation model. Based on the bi-spectrum in Figure 3, we observe that the cardiac frequency span is 1.75Hz to 1.85Hz. Consequently, the coupling component has a wide frequency range result from the variation of the heart rate frequency. The width of the peaks in Figure 3(b) also indicates the time varying frequency of the beating heart motion. If this variation is slowing changing statistics, the adaptive estimator can converge to the stationary signal. Furthermore, if the changing statistics are slow relative to the algorithm adaptation ability, the predictor provides the ideal time varying solution, which can track the heart motion exactly.

It is necessary to import the quadratic term in the estimation model. As shown in Figure 3, the peaks at the nonlinear coupling bi-frequencies are mainly distributed in the middle to high frequency band, which indicates that the nonlinear quadratic coupling regulates the dynamics of the beating heart. The quadratic term in the second order Volterra Series is appropriated to model this cardiac-respiration interaction such that with the quadratic term the model includes more motion details and approximates the frequency response better

3. Nonlinear adaptive Volterra Lattice Prediction Algorithm

In order to track the beating heart accurately, we design a future heart motion prediction algorithm for the robot control architecture. The prediction is to approximate the beating heart motion by the second order Volterra model with historical heart motion measurements. The Volterra Model is updated recursively whenever a new piece of data comes to update the coefficient vector, which is achieved through an adaptive Volterra filter in a nonlinear way. Since the computation complexity is increased significantly by the nonlinear model, we import a Fast RLS algorithm to reduce the computation burden and meet the real time needs.

3.1 Second order Volterra model

The general formulation of the second-order Volterra Series Model is given by:

\begin{array}{l} y [n] = \sum_{m_{1} = 0}^{N - 1} h_{1} (m_{1}; n) x (n - m_{1}) \\ + \sum_{m_{1} = 0}^{N - 1} \sum_{m_{2} = 0}^{N - 1} h_{2} (m_{1}, m_{2}; n) x (n - m_{1}) x (n - m_{2}) \\ = H^{T} [n] X [n] \end{array}

(2)

where X[n] is the input vector of the second-order Volterra Model

\begin{array}{l} X [n] = [x (n), x (n - 1), \dots, x (n - N + 1), x^{2} (n), \dots \\ , x (n) x (n - N + 1), x^{2} (n - 1), \dots, x^{2} (n - N + 1)]^{T} \end{array}

(3)

And the coefficient vector H[n] is defined as

\begin{array}{l} H [n] = [h_{1} (0; n), h_{1} (1; n), \dots, h_{1} (N - 1; n), h_{2} (0, 1; n), \\ \dots, h_{2} (0, N - 1; n), h_{2} (1, 1; n), \dots, h_{2} (N - 1, N - 1; n)]^{T} \end{array}

(4)

3.2 Fast Second Order Volterra Lattice -Recursive Least Square algorithm

The object of the algorithm is to approximate the desired future step heart motion signal, d[n], adaptively through the second order Volterra series expansion in the most recent N samples, denoted as d̂[n] in (6), by minimizing the accumulated error signal defined as a cost function

J [n] = \sum_{k = 0}^{n} λ^{n - k} {| d [k] - {\hat{H}}^{T} [n] X [k] |}^{2}

(5)

where

\hat{d} [n] = {\hat{H}}^{T} [n] X [n]

(6)

By minimizing the cost function at each time instant n, the estimated coefficient vector Ĥ[n] is optimally solved as

\hat{H} [n] = C^{- 1} [n] M [n]

(7)

Where

$C [n] = \sum_{k = 0}^{n} λ^{n - k} X [k] X^{T} [k]$ is the autocorrelation matrix of the input vector and $M [n] = \sum_{k = 0}^{n} λ^{n - k} X [k] d [k]$ is the cross

correlation vector between the input vector and the desired response.

The development of the conventional recursive least square (RLS) adaptive Volterra estimator can be achieved easily through the same procedure employed to derive the RLS adaptive filters employing linear system models [26]. However, the computational complexity of this algorithm requires O(N⁴) multiplications per iteration. Matthews et al. [33 –35] proposed nonlinear lattice-based fast RLS could simplify the computation to O(N³) multiplications per iteration. This algorithm is not only computationally efficient but also numerically stable. It is the extension of the linear lattice recursive least square algorithm. From [34, 35] we find the complete derivation of this algorithm. In addition, [36, 37] provide a complete discussion of the least square lattice algorithm in general and their properties.

We develop the Volterra lattice filter in the same way the linear lattice filter was developed. First, in order to develop the lattice parameterization of the second order Volterra model, it is convenient to visualize the nonlinear prediction problem as a multichannel linear prediction problem by reconstructing the input vector X[n] into the matrix form X_N[n]. Each column of the matrix can be viewed as a different channel. Data in each channel is the orthogonal basis for the vector space spanned by input data.

X_{N} [n] = [\begin{matrix} x (n) & x (n - 1) & x (n - 2) & \dots & x (n - N + 1) \\ x^{2} (n) & x^{2} (n - 1) & x^{2} (n - 2) & \dots & x^{2} (n - N + 1) \\ ​ & x (n) x (n - 1) & x (n - 1) x (n - 2) & ​ & x (n - N + 2) x (n - N + 1) \\ ​ & ​ & x (n) x (n - 2) & ​ & x (n - N + 3) x (n - N + 1) \\ ​ & ​ & ​ & ​ & ⋮ \\ ​ & ​ & ​ & ​ & x (n) x (n - N + 1) \end{matrix}]

(8)

Second, the main difference between a linear Lattice filter and a Volterra Lattice Filter is the number of elements in the backward and forward prediction vectors during order updates. The Volterra Lattice Filter propagates its order in the same way as a linear Lattice filter (Table 2 and Figure A.1 in Appendix A show us the distinguished facts). The vector dimension in the linear lattice filter is always the same during the order update, while in the Volterra Lattice Filter, the dimension of the backward and forward prediction vector is increased by 1 after each order update until the dimension equals N (N is the order number in the algorithm). The order numbers are 4 and 5 for the two datasets, which indicates that at the time instant n, the dimension of the backward and forward prediction vectors increases up to 4 or 5 in order to collect all the internal modes of the beating heart motion for prediction purposes.

Table 2.

Comparison of linear Lattice filter and Volterra Lattice Filter [33]

Linear Lattice Filter
f₀(n) = b₀(n) = x(n)
$b_{i}^{(n)} = b_{i - 1} (n - 1) - k_{i}^{b^{T}} f_{i - 1} (n)$
$f_{i}^{(n)} (n) = f_{i - 1} (n) - k_{i}^{f^{T}} (n) b_{i - 1} (n - 1)$
Volterra Lattice Filter
$\begin{matrix} {\underline{f}}_{0} (n) = {\underline{b}}_{0} (n) = {[x^{2} (n) x^{2} (n)]}^{T} \\ f_{0}^{i} (n) = x (n) x (n - j + 2), j = 3, \dots, N + 1 \end{matrix}$
${\underline{b}}_{i} (n) = [\begin{matrix} b_{i} (n) = {\underline{b}}_{i - 1} (n - 1) - k_{i}^{b^{T}} {\underline{f}}_{i - 1} (n) \\ b_{i}^{i + 2} (n) = f_{i - 1}^{(i + 2)} (n) - k_{i - 1}^{b {(i + 2)}^{T}} (n) {\underline{f}}_{i - 1} (n - 1) \end{matrix}]$
${\underline{f}}_{i} (n) = [\begin{matrix} f_{i} (n) = {\underline{f}}_{i - 1} (n - 1) - k_{i}^{b^{T}} {\underline{b}}_{i - 1} (n - 1) \\ f_{i}^{i + 2} (n) = f_{i - 1}^{(i + 2)} (n) - k_{i}^{f {(i + 2)}^{T}} (n) {\underline{b}}_{i - 1} (n - 1) \end{matrix}]$

The derivations of the recursive equations to update the forward and backward reflection coefficient matrices and the auxiliary coefficient vectors are the same as the method used to derive the linear lattice filter. We formulize the complete fast Volterra Lattice algorithm in Appendix A

4. Experiment and results discussion

In this section, we evaluate three types of predictive filters: adaptive AR Filters (AR-F), time varying Fourier Series with Extended Kalman Filters (FS-EKF) and Fast second order Volttera Lattice recursive least square Filters (FVL-F). The predictive filters are tested using two different sets of pre-recorded data. We attempt to investigate the influence of different heart motion models on prediction accuracy. We derive the prediction results through fine-tuning the parameters to minimize the prediction errors for the three algorithms.

We summarize the RMS prediction errors in Table 3. Figure 4 shows the prediction results from the first dataset and Figure 5 shows the prediction results from the second dataset.

Figure 4.

(a) Heart Motion Prediction Algorithm Results in X-axis using First Dataset; (b) Heart Motion Prediction Algorithm Results in Y-axis using First Dataset; (c) Heart Motion Prediction Algorithm Results in Z-axis using First Dataset

Figure 5.

(a) Heart Motion Prediction Algorithm Results in X-axis using Second Dataset; (b) Heart Motion Prediction Algorithm Results in Y-axis using Second Dataset; (c) Heart Motion Prediction Algorithm Results in Z-axis using Second Dataset

Table 3.

Prediction results

RMS prediction errors [μm]	AR-F	FS-EKF	FVL-F
The 1st dataset
X axis	135	99.4	27.5
Y axis	84	56.1	49.4
Z axis	165	64.4	19.6
The 2nd dataset
X axis	123	120	87
Y axis	157	89.8	13.1
Z axis	115	63.7	27.4

We find that the FVL-F algorithm has the best prediction accuracy when the nonlinearity of the phase coupling in the heart motion is relatively high. The performance of the FVL-F algorithm is evaluated through two distinguished characteristics of the beating heart motion: time-varying heart rate and phase coupling representation

Time varying heart rate; All three methods have the ability to compensate for to some extent the heart motion changes through different ways. The FS-EKF adapts the cardiac motion frequency variation through changing the frequency state in the model to update from one time sample to another. The AR-F model varies its filter parameters to fit a new block of data if the heart motion frequency varies. The FVL-F changes both its linear part kernel coefficients and its quadratic part kernel coefficients to take into account a change in heart motion frequency.

Phase coupling representation; The FVL-F overcomes the shortcomings of the AR-F and FS-EKF in the coupling between cardiac motion and respiratory motion, which are the two components of the heart motion. FS-EKF, a method in the frequency domain, is the linear combination of the sinusoid, excluding the quadratic part. The AR-F model is also a linear combination of past observations and is not able to adapt itself to newly produced frequency components. Judging from the prediction results, we can conclude that FVL-F is the best predictor among the three algorithms, because it provides the mechanism to conquer the phase coupling between the two dominant parts in the beating heart motion such that it can make an accurate prediction when the heart motion changes quickly. This is the property needed in robot tracking of beating heart surgery in the future.

Figure 6 shows the prediction of the most abrupt peak of the heart motion. The accurate prediction of this peak indicates the proposed algorithm has the ability to precisely predict the heart motion even in the extreme irregular case. A good property of the proposed prediction method is that it lays a good foundation for tracking control.

Figure 6.

Heart Motion Peak prediction comparison results

Figure A1.

Signal flow Diagram of Volterra Lattice Filter

5. Conclusions

In this paper, we aim to improve tracking performance by increasing the heart motion prediction accuracy using a nonlinear model with fast adaptive techniques to guarantee online real time computation. In this research, a bi-spectral method is used to analyse the phase-coupling nature between the respiration component and the heart beat component. The nonlinearity phenomenon dominates the beating heart motion. A novel nonlinear prediction method is proposed based on the Volterra Series Model. The fastest algorithm for the predictor is realized through the nonlinear lattice structure to satisfy the real time needs. A comparison of the experiment results is provided to evaluate the proposed algorithm using two different vivo datasets. The FVL-F algorithm shows convincing results and properties. When heart rate varies, FVL-F will change its kernel coefficients to adaptively follow the time-varying statistics of the beating heart motion. When considering phase coupling, FVL-F has a quadratic part in its structure to model the phase coupling effect and make better predictions through it.

In the future, the control strategy with the nonlinear Second Order Volterra Lattice Filter prediction will be implemented in a 3D test bed. More precise tracking performance by assisted robotics in beating heart surgery is expected by using the proposed prediction algorithm. Secondly, the effect of different positions of POI on the surface of the heart on the prediction algorithm will be investigated.

6. Appendix A

Table A1.

Pseudo code of FSOVL-RLS algorithm

1. Time initialization:
$α_{m} (0) = 1, f_{m}^{f} (0) = r_{m}^{b} (b) = {δ I}_{m + 2}, δ > 0 m = 0, 1, ... N - 1$
2. Order initialization: n = 1
$α_{0} (0) = 1, e_{0} (n) = d (n) = {\underline{f}}_{0} (n) {\underline{b}}_{0} (n) = {[x (n), x^{2} (n)]}^{T}$
$f_{0}^{(j)} (n) = x (n) x (n - j + 2), j = 3, ..., N + 1, r_{0}^{f} (n) = r_{0}^{b} (n) = {λ r}_{0}^{f} (n - 1) + {\underline{f}}_{0} (n) {\underline{f}}_{0}^{T} (n)$
3. For n ≥ 0 Set i = 0 to i = M −1
$α_{i} (n) = α_{i - 1} (n) - {\underline{b}}_{i - 1}^{T} (n) r_{i - 1}^{- b} (n) {\underline{b}}_{i - 1} (n)$ ,
$k_{i}^{f} (n) = k_{i}^{f} (n - 1) + \frac{r_{i - 1}^{- b} (n - 1) {\underline{b}}_{i - 1} (n - 1)}{α_{i - 1} (n - 1)} [{\underline{f}}_{i - 1}^{T} (n) - {\underline{b}}_{i - 1}^{T} (n - 1) k_{i}^{f} (n - 1)]$
$k_{i}^{b} (n) = k_{i}^{b} (n - 1) + \frac{r_{i - 1}^{- f} (n - 1) {\underline{f}}_{i - 1} (n)}{α_{i - 1} (n - 1)} [{\underline{b}}_{i - 1}^{T} (n - 1) - {\underline{f}}_{i - 1}^{T} (n) k_{i}^{b} (n - 1)]$
$f_{i} (n) = {\underline{f}}_{i - 1} (n - 1) - k_{i}^{b^{T}} {\underline{b}}_{i - 1} (n - 1), b_{i} (n) = {\underline{b}}_{i - 1} (n - 1) - k_{i}^{b^{T}} {\underline{f}}_{i - 1} (n)$
$k_{i}^{y} (n) = k_{i}^{y} (n - 1) + \frac{r_{i - 1}^{- b} (n - 1) {\underline{b}}_{i - 1} (n)}{α_{i - 1} (n - 1)} [e_{i - 1} (n) - {\underline{b}}_{i - 1}^{T} (n) k_{i}^{y} (n - 1)]$
$e_{i} (n) = e_{i - 1} (n) - k_{i}^{y^{T}} (n) {\underline{b}}_{i} (n)$
For j = 3,…,N-1
$k_{i}^{f (j)} (n) = k_{i}^{f (j)} (n - 1) + \frac{r_{i - 1}^{- b} (n - 1) {\underline{b}}_{i - 1} (n - 1)}{α_{i - 1} (n - 1)} [f_{i - 1}^{(j)} (n) - {\underline{b}}_{i - 1}^{T} (n - 1) k_{i}^{f (j)} (n - 1)]$
$f_{i}^{(j)} = f_{i - 1}^{(j)} (n) - k_{i}^{f {(j)}^{T}} {\underline{b}}_{i - 1} (n - 1)$
$k_{i}^{b (i + 2)} (n) = k_{i}^{b (i + 2)} (n - 1) + \frac{r_{i - 1}^{f} (n) {\underline{f}}_{i - 1} (n)}{α_{i - 1} (n - 1)} [f_{i - 1}^{j + 2} (n) - {\underline{f}}_{i - 1}^{T} (n) k_{i}^{b (i + 2)} (n - 1)]$
$b_{i}^{(i + 2)} (n) = f_{i - 1}^{(i + 2)} (n) - k_{i}^{b {(i + 2)}^{T}} (n) {\underline{f}}_{i - 1} (n)$
${\underline{f}}_{i} (n) = [\begin{matrix} f_{i} \\ f_{i}^{(i + 2)} (n) \end{matrix}], {\underline{b}}_{i} (n) [\begin{matrix} b_{i} \\ b_{i}^{(i + 2)} (n) \end{matrix}]$
$r_{i}^{- f (n)} = λ^{- 1} r_{i}^{- f} (n - 1) - \frac{λ^{- 2} r_{i}^{- f} (n - 1) {\underline{f}}_{i}^{T} (n) r_{i}^{- f} (n - 1)}{α_{i} (n - 1) + λ^{- 1} {\underline{f}}_{i}^{T} (n) r_{i}^{- f} (n - 1) {\underline{f}}_{i}^{T} (n)}$
$r_{i}^{- b (n)} = λ^{- 1} r_{i}^{- b} (n - 1) - \frac{λ^{- 2} r_{i}^{- b} (n - 1) {\underline{b}}_{i} (n) {\underline{b}}_{i}^{T} (n) r_{i}^{- b} (n - 1)}{α_{i} (n - 1) + λ^{- 1} b_{i}^{T} (n) r_{i}^{- b} (n - 1) b_{i}^{T} (n)}$
For i = M
$k_{M}^{y} (n) = k_{M}^{y} (n - 1) + \frac{r_{M - 1}^{- b} (n - 1) {\underline{b}}_{M - 1} (n)}{α_{M - 1} (n)} [e_{M - 1} (n) - {\underline{b}}_{M - 1} (n) k_{M}^{y} (n - 1)]$
$e_{M} (n) = e_{M - 1} (n) - k_{M}^{y^{T}} (n) {\underline{b}}_{M - 1} (n)$

f_m(n) forward prediction error vector	e_m(n) joint process estimation error at the mth stage
r_b (n) Backward prediction error correlation matrix
b_m(n) backward prediction error vector	f_m^(j)(n) Auxiliary prediction error
k_f^m(n) forward reflection coefficient matrix	k^y_m(n) Joint process estimation coefficient error
α_m(n) likelihood variable	k^b_m^(j)(n) Auxiliary reflection coefficient vector
k_m^b (n) backward reflection coefficient matrix	r_m^f (n) Forward prediction error correlation matrix
k^f_m⁽ ^j)(n) coefficient vector(j=m+2,…, N+1)	b⁽_m^j ⁾(n) Auxiliary prediction error

Footnotes

7. Acknowledgments

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

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

References

Newman

M. F.

Kirchner

J. L.

Phillips-Bute

Gaver

Grocott

and 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, 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, vol. 23, no. 3, pp. 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, Advanced Robotics, ICAR'05, Proceedings. 12th International Conference on, pp. 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

Ortmaier

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

Franke

T. J.

Bebek

and Cavusoglu

M. C.

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

Franke

T. J.

Bebek

and Cavusoglu

M. C.

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

10.

King

A. P.

Rhode

K. S.

Razavi

R. S.

and Schaeffter

T. R.

(2009), An adaptive and predictive respiratory motion model for image-guided interventions: Theory and first clinical application, Medical Imaging, IEEE Transactions on, vol. 28, no. 12, pp. 2020–2032.

11.

Thakral

Wallace

Tomlin

Seth

and Thakor

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

12.

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 vol. 21, no. 1, pp. 67–79.

13.

Richa

Bó

A. P. L.

and Poignet

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

14.

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, vol. 28, no. 10, pp. 1355.

15.

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, vol. 2, no. 5, pp. 888–896.

16.

Bachta

Renaud

Cuvillon

Laroche

Forgione

and Gangloff

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

17.

Shiogai

Stefanovska

and McClintock

P. V. E

(2010), Nonlinear dynamics of cardiovascular ageing, Physics Reports, vol. 488, no. 2–3, pp. 51–110.

18.

Akay

, editor (2000), Nonlinear Biomedical Signal Processing, Volume II: Dynamic Analysis and Modeling. Wiley-IEEE Press; New York

19.

Mitsis

G. D.

Markakis

M. G.

and Marmarelis

V. Z.

(2009), Nonlinear Modeling of the Dynamic Effects of Infused Insulin on Glucose: Comparison of Compartmental With Volterra Models, IEEE Trans Biomed Eng[J], vol. 56, no. 10, pp. 2347–58

20.

Korenberg

M. J.

(1991), Parallel cascade identification and kernel estimation for nonlinear systems, Ann. Biomed. Eng, vol. 19, no. 4, pp. 429–455

21.

Korenberg

M. J.

and Paarmann

L. D.

(1991), Orthogonal approaches to time series analysis and system identification, IEEE Signal Proc., vol. 8, no. 3, pp. 29–43

22.

Marmarelis

P. Z.

(1993), Identification of nonlinear biological systems using Laguerre expansions of kernels, Ann. Biomed, vol. 21, pp. 573–589

23.

Zhe

Brown

E. N.

and Barbieri

(2010), Characterizing Nonlinear Heartbeat Dynamics Within a Point Process Framework, IEEE Trans. on Biomedical Engineering, vol. 57, no. 6, pp. 1335–1347

24.

Tuna

E. E.

(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://robotics.case.edu/publications.html#thesis

25.

Servant

Logier

Mouster

and Goudemand

(2008), Heart rate variability. Applications in psychiatry, L'Encéphale, vol. 35, no. 5, pp. 423–428. doi:10.1016/j.encep.2008.06.016.

26.

Priestley

M. B.

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

27.

Schäfer

Rosenblum

M. G.

Kurths

and Abel

H. H.

(1998), Cardiac motion synchronized with ventilation, Nature vol. 39, no. 6673, pp. 239.

28.

Lotri

M. B.

and Stefanovska

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

29.

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, vol. 86, no. 9, pp. 1749–1752.

30.

Jamšek

Stefanovska

McClintock

P. V. E.

and Khovanov

I. A.

(2003), Time-phase bispectral analysis, Physical Review E, vol. 68, no. 1, pp. 016201.

31.

Jamšek

Stefanovska

and McClintock

P. V. E.

(2004), Nonlinear cardio-respiratory interactions revealed by time-phase bispectral analysis, Physics in Medicine and Biology, vol. 49, pp. 4407.

32.

Fackrell

J. W. A

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

33.

Mathews

V. J.

(2000), Polynomial Signal Processing, John Wiley Sons Inc.

34.

Lee

and Mathews

(1993), A fast recursive least squares adaptive second order Volterra filter and its performance analysis, IEEE trans. Signal processing, vol. 41, no. 3, pp. 1087–1873

35.

Sayed

and Mathews

(1994), Lattice Algorithms for Recursive Least Squares Adaptive Second Order Volterra Filtering, IEEE transactions on circuits and systems analog and digital signal processing, vol. 41, no.3

36.

Haykin

(2001), Adaptive Filter Theory, 4th ed. Upper Saddle River, New JerseyPrentice Hall

37.

Ali

(2008), Sayed Adaptive Filters, A JOHN WILEY & SONS, INC., PUBLICATION