Real-time prediction and gating of respiratory motion in 3D space using extended Kalman filters and Gaussian process regression network

The prediction as well as the gating of respiratory motion have received much attention over the last two decades for reducing the targeting error of the radiation treatment beam due to respiratory motion. In this article, we present a real-time algorithm for predicting respiratory motion in 3D space and realizing a gating function without pre-specifying a particular phase of the patient’s breathing cycle. The algorithm, named EKF-GPRN+ , first employs an extended Kalman filter (EKF) independently along each coordinate to predict the respiratory motion and then uses a Gaussian process regression network (GPRN) to correct the prediction error of the EKF in 3D space. The GPRN is a nonparametric Bayesian algorithm for modeling input-dependent correlations between the output variables in multi-output regression. Inference in GPRN is intractable and we employ variational inference with mean field approximation to compute an approximate predictive mean and predictive covariance matrix. The approximate predictive mean is used to correct the prediction error of the EKF. The trace of the approximate predictive covariance matrix is utilized to capture the uncertainty in EKF-GPRN+ prediction error and systematically identify breathing points with a higher probability of large prediction error in advance. This identification enables us to pause the treatment beam over such instances. EKF-GPRN+ implements a gating function by using simple calculations based on the trace of the predictive covariance matrix. Extensive numerical experiments are performed based on a large database of 304 respiratory motion traces to evaluate EKF-GPRN+ . The experimental results show that the EKF-GPRN+ algorithm reduces the patient-wise prediction error to 38%, 40% and 40% in root-mean-square, compared to no prediction, at lookahead lengths of 192 ms, 384 ms and 576 ms, respectively. The EKF-GPRN+ algorithm can further reduce the prediction error by employing the gating function, albeit at the cost of reduced duty cycle. The error reduction allows the clinical target volume to planning target volume (CTV-PTV) margin to be reduced, leading to decreased normal-tissue toxicity and possible dose escalation. The CTV-PTV margin is also evaluated to quantify clinical benefits of EKF-GPRN+ prediction.


Introduction
The motion of thoracic and abdominal tumours induced by respiratory motion can significantly compromise the conformal dose delivery in radiotherapy (Ozhasoglu and Murphy 2002, Keall et al 2006, Poulsen et al 2010. Motion-adaptive radiotherapy aims to deliver a conformal dose to the target tumour by compensating for its motion in real-time. Prediction of respiratory motion is necessary in real-time for the motion compensation and it has received much attention for over a decade (Sharp et al 2004, Ruan and Keall 2010, Krauss et al 2011, Lee et al 2012, Purdy et al 2012, Li and Sharp 2013, Dürichen et al 2014. Motion prediction approaches can be classified into two groups: model-based and model-free (Ernst et al 2013). The former predicts respiratory motion based on a mathematical motion model. The Kalman filter based on a state-space model is an example and it is computationally simple and responsive to irregular changes in respiratory motion. Model-free approaches do not assume an explicit mathematical model and predict the future position by learning a mapping from previous observations with, for example, support vector regression (SVR) and artificial neural networks (ANNs) (Riaz et al 2009, Murphy andPokhrel 2009). While the model-free algorithms are effective for learning and adapting to slowly time-varying dynamics of respiratory motion, their prediction can largely deviate from the true position over a period of irregular respiratory motion. Recently, a real-time prediction algorithm has been proposed that utilizes both model-based and model-free approaches by combining them in a cascade form (Bukhari and Hong 2015). The algorithm exploits the advantages of both of the approaches in predicting respiratory motion.
In this article, we present a real-time algorithm for predicting respiratory motion in 3D space and realizing a gating function without pre-specifying a particular phase of the patient's breathing cycle. The algorithm, named EKF-GPRN + , first employs an extended Kalman filter (named LCM-EKF) independently along each coordinate to predict the respiratory motion and then uses Gaussian process regression network (GPRN) to correct the error of LCM-EKF prediction in 3D space. The GPRN is a nonparametric Bayesian algorithm for modeling input-dependent correlations between output variables in multi-output regression (Wilson et al 2012, Wilson 2014. We employ variational inference with mean field approximation for GPRN to compute an approximate predictive mean and covariance matrix of the prediction errors of LCM-EKF (Nguyen and Bonilla 2013). The approximate predictive mean is used to correct the prediction error of LCM-EKF. The trace of the approximate predictive covariance matrix is utilized to capture the uncertainty in EKF-GPRN + prediction error and systematically identify breathing points with a higher probability of large prediction error in advance. This identification enables us to pause the treatment beam over such instances. EKF-GPRN + implements a gating function by using simple calculations based on the trace of the predictive covariance matrix. It is demonstrated through extensive numerical experiments that EKF-GPRN + reduces the prediction error to 38%, 40% and 40% in rootmean-square (RMS), compared to no prediction, at lookahead lengths of 192 ms, 384 ms and 576 ms, respectively. It is observed from the experiments that the use of the input-dependent correlations between the prediction error of EKF in 3D space plays a key role in improving the prediction accuracy. The error reduction allows the clinical target volume to planning target volume (CTV-PTV) margin to be reduced, leading to decreased normal-tissue toxicity and possible dose escalation (Vedam et al 2003, Keall et al 2006. We evaluate the CTV-PTV margin as well to demonstrate clinical benefits of EKF-GPRN + prediction. It is also demonstrated through numerical experiments that EKF-GPRN + can improve prediction accuracy and reduce the CTV-PTV margin further by employing the gating function, albeit at the cost of a reduced duty cycle.

Methods
In this section, we present EKF-GPRN and EKF-GPRN + prediction algorithms, starting with a brief introduction to LCM-EKF, EKF-GPR and EKF-GPR + prediction algorithms and GPRN model.

LCM model and LCM-EKF prediction algorithm
The LCM model characterizes one-dimensional respiratory motion locally as a projection of a circular motion to capture the temporal evolution of the respiratory motion (Hong et al 2011).
x k x k y k ,˙,˙, and ( ) Ω k denote the position, the velocity along the x-axis, the velocity along the y-axis, and the angular velocity, respectively. Defining the discrete-time state vector as ( ) [ ( ) ( ) ( ) ( )] = Ω k x k x k y k k x˙˙T, we can represent the local circular motion as a discrete-time state equation with a nonlinear mapping f that characterizes one-step evolution of the state vector in the discrete time. The evolution of position x(k) is the projection of the planar circular motion onto the x-axis. The y-axis is an auxiliary axis augmented to define the circular motion. The process noise vector ( ) k v is a zero-mean white noise sequence. The measurement equation of position x(k) is where z(k) is the position measurement at time k, and w(k) denotes the additive measurement noise. Based on the LCM model, the LCM-EKF algorithm employs the first-order EKF to recursively update the state estimate ˆ( ) | k k x with measurement z(k), and evaluates the l-step lookahead prediction ˆ( ) We use the LCM-EKF algorithm to predict the respiratory motion in 3D space by applying it independently along each coordinate. More details on the LCM-EKF prediction algorithm can be found in Hong et al (2011) and Hong and Bukhari (2014).

EKF-GPR and EKF-GPR + prediction algorithms
The EKF-GPR algorithm employs LCM-EKF to predict the respiratory motion and then uses the Gaussian process regression (GPR) algorithm to estimate and correct the error of the LCM-EKF prediction (Bukhari and Hong 2015). The structure of the EKF-GPR algorithm is illustrated in figure 1. First, let us denote the LCM-EKF prediction error by and let ˆ( ) | k k x denote the state estimate of LCM-EKF at time k. The GPR algorithm predicts the conditional distribution of under Gaussian assumptions. It evaluates the predictive mean (¯(ˆ( )) ξ | k k x LCM ) and the predictive variance as the estimate of ( ) ξ + | k l k LCM and utilize it to correct the LCM-EKF prediction ˆ( ) + | z k l k LCM to obtain the prediction of EKF-GPR, denoted by ˆ( ) The GPR algorithm is a supervised learning algorithm. It evaluates the predictive mean ¯(ˆ( )) ξ | k k x LCM and the predictive variance ( (ˆ( ))) ξ | k k x var LCM based on learning over a set of training samples {ˆ( ) i I LCM as illustrated in figure 1, where I denotes the index set of training samples. The predictive mean is the best point-prediction of the prediction error is also a useful measure for estimating the strength of EKF-GPR prediction error after the error correction. The algorithm named EKF-GPR + utilizes the predictive variance to capture the strength in the EKF-GPR prediction error and systematically identify breathing points with a higher probability of large prediction error in advance. This identification enables us to deactivate the treatment beam over such instances. EKF-GPR + implements a gating function by using simple calculations based on the predictive variance. The EKF-GPR + algorithm effectively reduces the prediction error in RMS by employing the gating function, albeit at the cost of a reduced duty cycle. We use EKF-GPR + (D%) to denote the EKF-GPR + algorithm with a nominal duty cycle parameter set to D. Note that EKF-GPR + (100%) is equivalent to EKF-GPR. More details on the EKF-GPR + algorithm can be found in Bukhari and Hong (2015).  We can use the EKF-GPR + algorithm for predicting and gating respiratory motion in 3D space. First, we apply the EKF-GPR algorithm independently along each coordinate for prediction and error correction. Then we use the sum of all the predictive variances as an approximate measure of the strength of EKF-GPR prediction error in 3D space after the error correction. EKF-GPR + utilize the sum to implement a gating function. The basic structure of EKF-GPR for 3D prediction is presented in figure 2.

Gaussian process regression networks (GPRN)
Modeling correlations between the output variables in the multi-output regression framework is known to improve regression performance (Breiman and Friedman 1997). Gaussian process (GP) models have been developed for the multi-output regression assuming fixed correlations between output variables (Bonilla et al 2008, Alvarez andLawrence 2011). Dürichen et al (2015) utilized a multi-task GP framework for motion compensation in radiation therapy. The fixed correlation assumption is, however, restrictive since correlations between output variables can be input-dependent.
GPRN is a new multi-output regression framework, which combines the structural properties of Bayesian neural networks with the nonparametric flexibility of Gaussian processes (Wilson et al 2012, Wilson 2014). This GPRN model represents each output variable as an adaptive mixture of latent functions shared across all the output variables. This representation enables GPRN to model input-dependent correlations between multiple output variables and improve regression performance. GPRN also underlies a highly expressive nonstationary kernel that learns nonstationary covariance structures between input and output variables. Such a nonstationary kernel can express more complicated structures in the data than stationary kernels that GPR can utilize (Rasmussen and Williams 2006). The following is a brief summary of GPRN. More details can be found in Wilson et al (2012) and Wilson (2014).
Suppose that we are given N samples of input-output pairs, is the corresponding p-dimensional output. Our goal is to predict the mean and the covariance matrix of the predictive distribution ( ( ) ) |D p y x at a new input x. The GPRN represents the output vector ( ) y x as where ( ) W x is a × p q matrix of mixing weights, each being an independent Gaussian process with ( ) ij w , and ( ) f x is a q-dimensional vector of latent functions, each being also an independent Gaussian process with ( ) Here, I q and I p are × q q and × p p identity matrices. This structure of GPRN is depicted in figure 3 where δ aw denotes the Kronecker delta. It can be seen from figure 3 that each output variable is an input-dependent mixture of latent Gaussian process basis functions with the mixing weights being also Gaussian processes. Each of the outputs ( ) y x i , conditioned on the mixing weights W( x), is a Gaussian process with kernel (Wilson 2014) , .
y a w j q ij a f a w ij w a w y 1 2 i j Note that this underlying GPRN kernel is an input-dependent mixture of kernels used in latent functions. This underlying kernel is highly expressive that can extract a nonstationary covariance structure in the data. Inference in GPRN is intractable. Variational inference methods have been usually utilized to approximate posterior distributions for GPRN inference (Wilson et al 2012, Nguyen andBonilla 2013). The approximate distributions are employed in learning hyperparameters , , , f w f y , where θ f and θ w denote the hyperparameters for kernels k f i and k w , respectively. The approximations are also used in deriving the mean and covariance of the predictive distribution. The approximate predictive mean of an output variable ( ) y x i and the approximate predictive covariance between two output variables ( ) y x i and ( ) y x j over a new input x are given, respectively, as¯( for all i, and for all i and j. Note that the approximate predictive mean and covariance depend on the first two moments of latent functions and weight functions.

EKF-GPRN and EKF-GPRN + prediction algorithms
The EKF-GPRN algorithm first employs LCM-EKF independently along each coordinate to predict the respiratory motion and then utilizes GPRN as a multi-output regression in 3D space to estimate and correct the errors of LCM-EKF prediction. The basic structure of EKF-GPRN for 3D prediction is presented in figure 4. In contrast to EKF-GPR in 3D (shown in figure 2), EKF-GPRN utilizes GPRN to jointly estimate and correct the error vector of LCM-EKF pre- to denote the l-step lookahead prediction of LCM-EKF and EKF-GPRN in 3D space, respectively. In addition, we use Then, we can use figure 1 to describe the structure of EKF-GPRN, replacing and GPRN, respectively. Note that the LCM-EKF block of EKF-GPRN in the figure consists of three LCM-EKFs implemented independently along each coordinate.
It has been reported that strong correlations exist between the state estimate ˆ( ) | k k x and the prediction error (Hong and Bukhari 2014). The EKF-GPR algorithm exploits this correlation in estimating , independently for each coordinate, and uses this estimate to correct the LCM-EKF prediction. We observe that strong correlations also exist between the components of the 3D error vector , as we will discuss later in this work. We adopt GPRN to learn and exploit these two different types of correlations in predicting the error Note that the GPRN framework can learn input-dependent correlations between multiple output variables as well as nonstationary covariance structures between input and output variables. While Bayesian models have been developed for multi-output regression to model fixed correlations between the output variables, GPRN is the only known Bayesian model that can represent the input-dependent correlations between output variables.
The GPRN algorithm in EKF-GPRN evaluates the predictive mean of using equations (8) and (9), respectively. We Note that x and ( ) y x in equation (5) The GPRN algorithm evaluates the predictive mean ¯(ˆ( )) ξ | k k x LCM and the predictive covariance to denote the trace of the predictive covariance matrix and employ it to quantify the strength of EKF-GPRN prediction error (after the error correction). The trace approximately corresponds to the mean square error (MSE) of the EKF-GPRN prediction in 3D space. This relation can be made more precise. First, we can express the MSE of LCM-EKF prediction as Suppose that GPRN learning is perfect and the predictive mean , respectively. Then, the EKF-GPRN prediction has zero mean and covariance equal to the predictive covariance as long as the predictive mean and covariance are reasonably accurate. It supports that we can use the trace of the predictive covariance as an approximate measure for estimating the strength of a possible prediction error of EKF-GPRN.
We utilize the trace of the predictive covariance matrix from the GPRN to identify breathing points with a higher probability of large prediction error in advance. This identification enables us to deactivate the treatment beam over such instances. We maintain a set of traces of the matrices sampled over a training window to determine a threshold on the trace. The threshold for a given duty cycle D% is determined as the sample point at which (100-D)% of the traces in the set lie above it. The trace associated with each of the subsequent predictions is then compared to this threshold. If the trace is greater than this threshold, we deactivate the treatment beam for that prediction. We name the EKF-GPRN algorithm with this detection mechanism as the EKF-GPRN + algorithm. The EKF-GPRN + algorithm effectively reduces the prediction error in RMS by employing the gating function, albeit at the cost of a reduced duty cycle. We use EKF-GPRN + (D%) to denote the EKF-GPRN + algorithm with a nominal duty cycle parameter set to D. Note that EKF-GPRN + (100%) is equivalent to EKF-GPRN.

Implementation
A simple grid search is performed over the parameter space of each algorithm to find the optimum parameters that minimize the prediction error of the corresponding algorithm. For the grid search, we utilized the first 10 min of 31 respiratory motion traces, each from a different patient, in a large database of 304 respiratory motion traces that will be described in section 3.1. We implemented the LCM-EKF algorithm as described in Hong and Bukhari (2014). The dynamics of motion along different coordinates, however, are different. We optimized the process noise parameters of LCM-EKF separately for each coordinate. The magnitude of respiratory motion along coordinate 2 is the largest, and the optimum values of process noise parameters for this coordinate are identical with those obtained over PCAprocessed one-dimensional traces in Hong and Bukhari (2014) . The variance of the measurement noise w(k) was set to 10 −4 for all the coordinates.
Inference in GPRN is intractable and we use the mean-field approximation for the inference which is a variational Bayesian method. The method approximates the posterior distribution as a factorized distribution over the parameters of GPRN with each factor distribution being a Gaussian (Nguyen and Bonilla 2013). As input to GPR and GPRN, we use only the velocity components in the state estimate vector of LCM-EKF. To be more specific, we use the estimate of [ ( ) ( )] x k y ki for GPRN. The number of latent functions in GPRN is an important factor which affects the prediction performance. We obtained the prediction performance of EKF-GPRN with the number of latent functions from 1 to 6 over the first 10 min of 31 traces for each lookahead length. The best prediction performance was achieved when the number of latent functions was equal to or greater than 4, and we set it to 4. EKF-GPRN + has been implemented in Matlab and run on a PC (Intel(R) Core(TM) i7-3770 CPU @ 3.40 GHz). In our implementation, we use the first 1 min data from 31 traces, each from a different patient, to train and obtain an initial GPRN model off-line. The GPRN algorithm is initialized with this initial GPRN model for all traces and it is updated every minute. Note that the GPRN algorithm learns the mapping from LCM-EKF state estimates to its prediction error. The mapping does not change fast in time and the GPRN update rate was high enough. The maximum number of iterations in learning hyperparameters for the GPRN update was set to 5 and the computation time for the update was a maximum of 1.5 sec. The update is not time-sensitive and it can be performed in the background. We use the updated GPRN model over a following 1 min interval. The computation time for LCM-EKF prediction in 3D is approximately 150 μsec and the total computation time is 370 μsec for completing one cycle of EKF-GPRN + with a sample data.
The first principal component of a 3D trace has been widely used in studies on the prediction of respiratory motion in 1D space. For a practical assessment of prediction over the first principal component, we implement an online PCA processing and evaluate real-time prediction of 3D respiratory motion in the principal component direction. We obtain the first principal component of a given trace and its direction through the online PCA processing over a 1 min sliding window, and use LCM-EKF and EKF-GPR algorithms to predict the first principal component and project the prediction back to the original 3D space. The prediction error is the difference between the prediction and the corresponding measurement in 3D space. The algorithms are referred to as LCM-EKF(PCA) and EKF-GPR(PCA). The structure of EKF-GPR(PCA) is shown in figure 5, whereas LCM-EKF(PCA) does not include the GPR components.

Materials
A large database of 304 respiratory motion traces is employed to evaluate the proposed algorithms. The 3D traces in the database were obtained during CyberKnife treatment at Georgetown University Hospital from a group of 31 patients (Ernst 2012). The database also includes PCA-processed 1D traces extracted over the longest interval in between the large couch movements. The PCA traces have been employed to evaluate 1D prediction algorithms (Ernst et al 2013, Hong andBukhari 2014). Likewise, we extracted 3D traces from the original traces over the intervals each corresponding to its PCA counterpart and we use them to evaluate our 3D prediction algorithms. The database lists the traces sorted in ascending alphanumeric order. For convenience in our presentation, we assigned the numeric IDs from 1 to 304 to these traces in the same order.
We implemented the proposed algorithms and performed experiments to evaluate their prediction performance over the 3D traces at lookahead lengths of 192 ms, 384 ms and 576 ms. In our experiments, we downsample the traces to 5.2 Hz for conformity to image-guided treatment practices including MRI-guided radiotherapy (Cervino et al 2011, Cho et al 2011, Tehrani et al 2013, Lagendijk et al 2014, Yip et al 2014. The experimental results in this section are based on the data that we obtained over each trace after the first 1 min.

Performance measures
We use the RMS prediction error and the CTV-PTV margin to evaluate the performance of prediction algorithms. First, let ( ) ( ) + k l z i denote the measurement at k + l for a trace i, and let ˆ( denote the magnitude of the 3D prediction error. We evaluate the RMS prediction error ( where K (i) is an index set of samples of trace i and ( ) | | K i denotes the number of indices in the set. Let ( ) j RMSE denote the patient-wise RMS prediction error of patient j, averaged across all traces of the patient. We obtain the patient-wise RMS error as where I( j ) is the index set of traces of patient j and ( ) | | I j is the number of traces of patient j. We evaluate the mean, standard deviation, 90th percentile point and maximum value of the set { ( )} ∈ j RMSE j J , where J is the index set of all patients. These statistics indicate the interpatient variability of RMS prediction error of an algorithm. The RMS prediction error averaged across all of the patients is obtained as where | | J is the total number of patients. The CTV-PTV margin is also evaluated to quantify clinical benefits of EKF-GPRN + . The margin is added around the CTV to ensure that the prescribed dose is actually delivered to the CTV, accounting for the effect of targeting error of the treatment beam (Vedam et al 2003, Keall et al 2006. The targeting error can be reduced by predicting the position of tumour and directing the treatment beam towards the predicted position. The error reduction allows the CTV-PTV margin to be reduced, leading to decreased normal-tissue toxicity and possible dose escalation.
In our experiment, we calculate the CTV-PTV margin that is required to maintain the CTV to lie within the PTV for 95% of each respiratory motion trace at a given lookahead length. Note that the CTV-PTV margin corresponds to the 95th percentile point of the set denotes the magnitude of the l-step lookahead prediction error of trace i as described above. Let m (i) denote the CTV-PTV margin of trace i for a given lookahead length. We obtain the patient-wise CTV-PTV margin averaged across all the traces of each patient. Let m( j ) denote the patient-wise CTV-PTV margin of patient j. For the set { ( )} ∈ m j j J , we evaluate its mean, standard deviation, 90th percentile point and maximum value.

Prediction performance of EKF-GPR + and EKF-GPRN + algorithms
Firstly, we evaluated the RMS prediction error of LCM-EKF, EKF-GPR + and EKF-GPRN + for each patient. Table 1 presents a statistical summary of the RMS error. Note that 'no prediction' in the table corresponds to the case that we employ the current position measurement itself as the lookahead prediction. The values in the parentheses in the table are the percent ratios relative to no prediction. The table shows that LCM-EKF, EKF-GPR and EKF-GPRN algorithms reduce the RMS error of the baseline no prediction approximately by 28%, 45% and 60%, respectively, across all of the lookahead lengths. Similar observations hold for the other statistics. These results confirm that EKF-GPRN reduces the prediction error as well as its inter-patient variability significantly and more effectively than EKF-GPR.
The left column of figure 6 presents scatter plots of the patient-wise RMS error of EKF-GPR and EKF-GPRN against the RMS error of no prediction for lookahead lengths of 192 ms, 384 ms and 576 ms. The two straight lines in each plot, corresponding to the 40% RMS error and 50% RMS error of no prediction, serve as reference lines for comparison against no prediction. The vertical difference between a pair of markers corresponds to the reduction in prediction error achieved by EKF-GPRN over EKF-GPR, and we can see that it increases as the lookahead length increases. In all of the plots, EKF-GPRN markers are always located below the corresponding EKF-GPR markers, indicating that EKF-GPRN reduces the RMS error of EKF-GPR over all of the 31 patients across all of the lookahead lengths. The EKF-GPR markers are populated slightly above the 50% RMS error line, whereas most of the EKF-GPRN markers are located below that line around 40%. This advantage of EKF-GPRN over EKF-GPR is statistically significant. We performed the paired Student's t-test to confirm the statistical significance of EKF-GPRN in prediction accuracy over EKF-GPR. The test rejected the null hypothesis with strong evidence for all the lookahead lengths with p-values less than 10 −6 . Table 1 also presents a statistical summary of the patient-wise RMS prediction error of EKF-GPR + and EKF-GPRN + for the duty cycles of 70%, 80% and 90%. The table shows that EKF-GPRN + (70%), EKF-GPRN + (80%) and EKF-GPRN + (90%) reduce the patient-wise RMS error approximately to 32%, 33% and 35%, respectively, compared to no prediction  across all of the lookahead lengths. These results correspond to an additional error reduction by 10% to 15% over EKF-GPR + in terms of the RMS prediction error of no prediction.
The threshold on the trace of the predictive covariance matrix of EKF-GPRN + was determined by setting the duty cycle parameter D to 70%, 80% and 90% in the above experiments. Recall that we use a set of traces of the approximate predictive covariance matrices from GPRN sampled over a 1 min training window to determine the threshold and employ this threshold over the following minute. As a result, the true duty cycle (D * ) is not exactly the same as set by the nominal parameter D. In our experiments, the true duty cycles have their values in an interval ((D-3)%, D%) for a given duty cycle D% for all the patients across all the lookahead lengths. It indicates that EKF-GPRN + controls the duty cycle with a high accuracy. The average of the true duty cycles were 68%, 79% and 89% for the nominal duty cycles of 70%, 80% and 90%, respectively, across all of the lookahead lengths. Table 2 presents a statistical summary of the patient-wise CTV-PTV margin. It shows that the average CTV-PTV margins required with no prediction are 1.30 mm, 2.35 mm and 3.23 mm, respectively, for lookahead lengths of 192 ms, 384 ms and 576 ms. EKF-GPRN reduces the average margins to 0.49 mm, 0.98 mm and 1.38 mm, respectively, which correspond to 38%, 42%, and 43% of the baseline no prediction. The table also shows that LCM-EKF and EKF-GPR algorithms reduce the margin to approximately 73% and 57% of no prediction, respectively, across all of the lookahead lengths. Note that the average margins in table 2 are between 1.65 and 1.75 times greater than the corresponding RMS prediction errors (RMSE P ) presented in table 1. Similar observations hold for the other statistics in the table.

CTV-PTV margin of EKF-GPR + and EKF-GPRN + algorithms
The right column of figure 6 presents scatter plots of the patient-wise CTV-PTV margins required with EKF-GPR and EKF-GPRN against that of no prediction for lookahead lengths of 192 ms, 384 ms and 576 ms. The two straight lines in each plot correspond to the 40% margin and 50% margin of no prediction. The distributions of the margin are not quite different from those of the RMS prediction error presented in the left column, except that their axis scales are different. As a result, the preceding discussion on the RMS prediction error similarly holds for the margin except that the vertical difference between a pair of markers (corresponding to the reduction in the CTV-PTV margin achieved by EKF-GPRN over EKF-GPR) increased approximately by a factor of 1.7. The advantage of EKF-GPRN over EKF-GPR in the margin is also statistically significant for all the lookahead lengths with p-values less than 10 −6 . Table 2 also presents a statistical summary of the patient-wise CTV-PTV margin required with EKF-GPR + and EKF-GPRN + . The table shows that EKF-GPRN + (70%), EKF-GPRN + (80%) and EKF-GPRN + (90%) reduce the patient-wise margin to 29%, 31% and 33%, respectively, compared to that of no prediction at the lookahead length of 192 ms, 32%, 34% and 37%, respectively, at 384 ms, and 34%, 36% and 38%, respectively, at 576 ms. These results correspond to an additional reduction of the CTV-PTV margin by 12% to 16% over EKF-GPR + relative to no prediction. Note that EKF-GPRN reduces the prediction error and the margin even more than EKF-GPR + (70%). Figure 7 presents scatter plots of the patient-wise CTV-PTV margin of EKF-GPRN, EKF-GPRN + (70%) and EKF-GPRN + (90%) against that of no prediction at the lookahead lengths of 192 ms, 384 ms and 576 ms. The two straight lines in each plot correspond to the 30% margin and 40% margin of no prediction. The scatter plots show that the CTV-PTV margin of EKF-GPRN + decreases as the duty cycle decreases for all of the 31 patients across all of the lookahead lengths, confirming that the gating function satisfies its purpose. Statistical summary of patient-wise CTV-PTV margin (mm) (figures in parentheses are percent ratios relative to no prediction).

Discussion
EKF-GPR and EKF-GPRN algorithms estimate the error of LCM-EKF prediction and correct it by exploiting the correlations between the state estimate and the prediction error of LCM-EKF. Both of the algorithms significantly improve prediction accuracy over LCM-EKF. Prediction of respiratory motion in 3D space provides an additional opportunity for EKF-GPRN to exploit the correlations existing between LCM-EKF prediction errors in different coordinates, resulting in further improvement in prediction accuracy over EKF-GPR. The correlations between the prediction errors are dependent on the state vector estimate ˆ( ) | k k x of LCM-EKF. Note that the latent functions in GPRN are input-dependent and sharing these latent functions across the output variables enables GPRN to represent the input-dependent correlations between the output variables. We utilized GPRN to model the input-dependent correlations between the prediction error of EKF, which is pivotal to improve the prediction accuracy in 3D space.
A showcase example is presented in figure 8 to highlight how effectively EKF-GPRN exploits the input-dependent correlations. Firstly, figure 8(a) illustrates a square grid of the correlation coefficients between the prediction errors of LCM-EKF in coordinates 1 and 2, where x and y-axis correspond to the velocity estimate ˆ( ) | x k k 1 along coordinate 1 and the velocity estimate ˆ( ) | x k k 2 along coordinate 2, respectively. The data was obtained for all sample points over trace 13 at a lookahead length of 384 ms (this trace corresponds to the trace with ID p7_f22_m3 in the CyberKnife treatment database). The grid consists of 129 squares, each with × − 0.16 mm sec 1.6 1 mm sec −1 in the input space and each square includes at least 10 samples. The figure shows that the LCM-EKF prediction errors have relatively strong correlations. The RMS prediction error of LCM-EKF for this trace is 1.93 mm at the lookahead length of 384 ms. EKF-GPR and EKF-GPRN reduce this error to 1.40 mm and 1.20 mm, respectively. Although EKF-GPR significantly reduces the LCM-EKF prediction errors, it cannot model the correlations between the prediction errors in 3D space. Figure 8(b) presents the distribution of input-dependent correlations between the EKF-GPR prediction errors after GPR error correction. We can see from the figure that the correlations become weaker but still maintain a certain level of strength even after the GPR error correction. The prediction errors of EKF-GPRN, in contrast, are almost decorrelated after the GPRN error correction; see figure 8(c). This example illustrates how effectively EKF-GPRN can exploit the input-dependent correlations between LCM-EKF prediction errors in improving its prediction accuracy.
Although not reported in this article, we implemented EKF-GPR and EKF-GPRN with different configurations and performed experiments to evaluate their performances. Through the experiments, we learned that the underlying covariance between LCM-EKF prediction error and its state estimates along each coordinate was not stationary and more complicated than can be reasonably expressed by simple GPR covariance functions. However, the advantage of the expressive underlying kernel of GPRN model contributed only in a small part to overall performance gain of EKF-GPRN over EKF-GPR. We learned in the experiments that the use of the input-dependent correlations between the output variables plays the key role in improving the prediction accuracy in 3D space.
The unscented Kalman filter (UKF) was developed as a nonlinear estimation method to overcome possible limitations arising from linearization that EKF relies on (Julier and Uhlmann 2004). We implemented an UKF based on the LCM model (named LCM-UKF) to compare with LCM-EKF. The comparison shows that LCM-UKF and LCM-EKF are comparable in prediction accuracy, but LCM-EKF is much faster in computation time than LCM-UKF. We also implemented LCM-EKF(PCA) and EKF-GPR(PCA) and evaluated their performances. The results show that LCM-EKF(PCA) and EKF-GPR(PCA) yield significantly worse prediction performance compared to LCM-EKF and EKF-GPR, respectively. It suggests that lateral motion is significant in the plane normal to the first principal component direction and realtime PCA prediction is not a good choice for predicting 3D respiratory motion.
Irregularities in respiratory motion are difficult to predict and prediction accuracy deteriorates over such events, irrespective of prediction algorithms. Traces 223-231 are examples with irregular patterns persisting over the whole period of treatment. The traces are belong to the patient that corresponds to the EKF-GPRN marker singly lying above the 50% lines in figure 6. It can be seen from figure 7 that, despite of such irregularities, EKF-GPRN + (70%) reduces the CTV-PTV margin from 60% to 40% that of no prediction. The marker at the far right (with a largest no prediction error) in figures 6 and 7 corresponds to traces 10-18. The traces represent respiratory motion with a large amplitude and high frequency, resulting in a largest prediction error and CTV-PTV margin.
Finally, we summarize the prediction and clinical advantages of EKF-GPRN + as follows. The patient-wise RMS prediction errors of no prediction are 0.79 mm, 1.43 mm and 1.96 mm, respectively, for lookahead lengths of 192 ms, 384 ms and 576 ms. EKF-GPRN + (70%) reduces the errors to 0.24 mm, 0.44 mm and 0.63 mm. The average CTV-PTV margins of no prediction are 1.30 mm, 2.35 mm and 3.23 mm, respectively, for lookahead lengths of 192 ms, 384 ms and 576 ms. EKF-GPRN + (70%) reduces the margins to 0.38 mm, 0.75 mm and 1.09 mm.

Conclusion
A real-time algorithm for predicting and gating respiratory motion in 3D space has been presented. The algorithm, named EKF-GPRN + , first employs an extended Kalman filter (LCM-EKF) along each coordinate to predict the respiratory motion and then uses a Bayesian GPRN model for multi-output regression to jointly estimate and correct the errors of the LCM-EKF prediction in all the coordinates. The EKF-GPRN + algorithm also implements an automated gating function without pre-specifying a particular region of the patient's breathing cycle. Specifically, EKF-GPRN + systematically identifies breathing points with a higher probability of large prediction error in advance so that we can deactivate the treatment beam over such instances. Extensive numerical experiments were performed based on a large database of 304 traces from 31 patients. The experimental results have shown that EKF-GPRN + improve prediction accuracy significantly and more effectively than EKF-GPR + , allowing the CTV-PTV margin to be reduced further. The experiments also confirmed that EKF-GPRN + controls the duty cycle with a high accuracy.
The extended Kalman filter seems more flexible and easier to implement information fusion from multiple imaging sources to guide radiation treatment compared to methods relying on learning. The EKF-GPRN + algorithm has an attractive structure for this purpose with an extended Kalman filter at its first stage. We plan to investigate this direction in the future.