A Respiratory Motion Prediction Based on Time-Variant Seasonal Autoregressive Model for Real-Time Image-Guided Radiotherapy

In radiation therapy, to deliver continuously a sufficient radiation dose to target volume yields a better therapeutic effect. While, avoiding an exposure to healthy tissues surrounding the target volume is also an important requirement for suppressing the adverse effect. Image-guided radiation therapy (IGRT) has potential to achieve the two requirements and as it’s application, stereotactic body radiation therapy (SBRT) has been used in clinic. In SBRT, the irradiated field is positioned with millimeter accuracy by proper daily setup. The accurate irradiation can allow the increase of radiation dose by ignoring the irradiation to the healthy tissues. Indeed, it has been reported that the treatment result of SBRT is comparable to the outcome from surgery [1].


Introduction
In radiation therapy, to deliver continuously a sufficient radiation dose to target volume yields a better therapeutic effect. While, avoiding an exposure to healthy tissues surrounding the target volume is also an important requirement for suppressing the adverse effect. Image-guided radiation therapy (IGRT) has potential to achieve the two requirements and as it's application, stereotactic body radiation therapy (SBRT) has been used in clinic. In SBRT, the irradiated field is positioned with millimeter accuracy by proper daily setup. The accurate irradiation can allow the increase of radiation dose by ignoring the irradiation to the healthy tissues. Indeed, it has been reported that the treatment result of SBRT is comparable to the outcome from surgery [1].
As mentioned above, the higher accuracy as irradiation allows the higher planning dose to the target. However, intra-fractional organ motion, such as respiratory motion of lung, often makes misalignment of the isocenter and the target volume during treatment fraction. For example, respiratory motion moves lung tumor over 10 mm per second [2,3]. In this case, the irradiation error caused by intra-fractional motion is not negligible in terms of adverse effect. Therefore, the respiratory motion management is a very important task in the field of radiation treatment [4].
To take into account the intra-fraction motion of lung tumor, some irradiation techniques have been investigated. The technique most widely used is the use of internal target volume (ITV) [5]. ITV includes internal margin determined by the intra-fractional organ motion. Consequently ITV covers the target volume CTV without misalignment as shown in Figure  1). However, radiation dose must be lower than the case of the irradiation without internal margin, because the irradiated area also covers normal tissues surrounding the target volume. On the other hand, a gating technique can give a high radiation dose to the lung tumor [6,7].    This is on-and-off irradiation to static region as shown in Figure 2. The radiation dose is delivered when the target volume is within the area planned preliminarily, and the irradiation interrupts if the target volume is out of the planned area. The gating technique can suppress the exposure to healthy tissues and allows high-dose-rate irradiation which is sufficient to yield a better therapeutic effect. However, instead of high-dose exposure, the gating takes longer treatment time to yields same therapeutic effect to SBRT due to the interruption.
An ideal irradiation to the moving target is to continuously irradiate a sufficient dose to the tumor which can be achieved by controlling the radiation beam to chase the moving target [8]. Such real-time tumor following (or chasing) irradiation yields an ideal therapeutic effect and can shorten the treatment duration. A schematic diagram of the tumor following irradiation is shown in Figure 3.
The tumor following irradiation requires the following two key techniques at least.    Figure 4. A schematic diagram of the tumor following irradiation with system latency. In this example, the total delay is 1 s and the delay distances the target tumor from the irradiated area over 10 mm.
The first technique for tumor localization has been achieved by using an X-ray fluoroscopic imaging system such as real-time tumor tracking system [7]. The second technique of beam-positioning systems can be realized by using dynamic multi-leaf collimator, robotic couch, robotic manipulator, and so on. However, the tumor following irradiation has not been developed clinically yet. This is because even we can obtain the precise lung tumor position in real-time, current radiotherapy machine has mechanical and computational time delays of up to about 1 s for controlling the irradiation field and image processing [9]. The latency can definitely affect badly on the irradiation accuracy [10], and must be compensated. Figure 4 shows a schematic diagram of the tumor following irradiation without delay compensation.
Prediction of future lung tumor motion is a typical solution for compensating the system latency and thus is a hot research topic for the tumor following irradiation. Then, the third technique is needed to develop the tumor following radiation therapy.
(iii) Real-time motion prediction Several prediction methods have thus been proposed for the respiratory motion. These include linear regression [11,12], extrapolation [11], artificial neural network [11,13], Kalman filter [11], nonlinear regression based on the Takens theorem [12,14], probabilistic modeling [15,16] and so on. More details can be found in a literature such as a survey study of the prediction of lung tumor motion [17]. A variety of the prediction approaches indicates that there is currently no best prediction method in clinical use because of the insufficient accuracy.
In general, the prediction is realized as an application of time series analysis. General prediction methods, such as autoregressive moving average (ARMA) model [18], require a stationarity of the target time series. However, the respiratory motion essentially has complex and time-varying characteristics. For example, the repetition of inhalation and exhalation naturally involves periodicity, but the periods are time-varying. Nevertheless the periodic component can still help to predict the motion because the past observed motion will arise again at periodic intervals. It means that we can predict the motion accurately if the period of target motion is obtained. Thus, the use of the periodicity can be a good approach to predict the respiratory motion.
The methods that focus on periodic components in respiratory motion have been developed as seasonal autoregressive (SAR) model-based method [19] and periodic ARMA model-method [20]. However, the periodicity in the respiratory motion fluctuates with time, and the adaptation of these methods to time-variant periodic variation seems insufficient yet. For example, SAR model-based method converts the target motion with time-variant period to a new motion with constant period to use the SAR model properly. But if the conversion from time-variant to constant period is incomplete, the prediction accuracy can decrease. On the other hand, the periodic ARMA estimates the period of the target motion by using long historical samples of 60 s. However, the use of long historical samples causes a hysteresis and cannot trace the change of the periodicity. Therefore, to adapt to the time-varying periodicity still remains as a challenge to accurately predict the respiratory motion.
In this chapter, to predict the respiratory motion by use of its periodicity, a time-variant seasonal autoregressive (TVSAR) model-based method is proposed [21]. The proposed method can model the time-variant periodicity more precisely, and we adopt only a small number of samples to suppress the hysteresis. In the next section, we briefly describe the characteristics of the lung tumor motion. Then, in section 3, TVSAR model is explained. The model is based on SAR, but is newly developed to adapt to the fluctuated periodicity by using unequally-spaced intervals instead of multiples of the constant period. To show the prediction performance of TVSAR model-based method, experimental result by using some clinical data sets of lung tumor motion is described in section 4. Section 5 provides concluding remarks. Figure 5 shows an example of time series of respiratory motion of lung tumor. The time series of tumor motion was observed as a location of the golden fiducial marker implanted into around the tumor, by using a kV X-ray fluoroscopic imaging system known as RTRT system, with sampling frequency 30 Hz, at Hokkaido University Hospital, Sapporo, Japan [7].

Respiratory motion of lung tumor
In general, a lung tumor motion involves a periodical component because breathing is composed of repetition of inspiration and expiration. corresponding to the respiratory cycle can be found at 0.33 Hz approximately in this example. This means that the tumor motion has a periodical component induced by the 3 s respiratory cycle in average.
As is clear from the analysis, one of dominant components in respiratory motion is periodic variation due to repetition of inhalation and exhalation. However, the period of the motion is not constant, but a time varying function. Such quasi-periodical nature of the motion is simply found as fluctuation of peak-to-peak intervals. For example, time interval between a peak and next peak differs each other even if the amplitude variation is sufficiently small as shown in Figure 7. This suggests that the dominant frequency is time-varying.
Recalling that the breathing period is time-variant, short-time Fourier transform (STFT), instead of the normal long-time Fourier transform, was performed on the motion example. The time-variation of the frequency spectrum is shown in Figure 8

Prediction methods
The periodicity found in the respiratory motion is a typical one of nonlinear and complex nature but also it can be useful to predict itself. That is, the past observed patterns in the motion will repeatedly arise at constant period. In this section, seasonal autoregressive (SAR) model is explained as a method with use of periodicity. Also, a limitation of general SAR model for lung tumor motion prediction is shown. Then, time-variant SAR model is introduced as a method designed for tumor motion prediction.

Seasonal autoregressive (SAR) model
Seasonal autoregressive integrated moving-average (SARIMA) model [22] is a general expression of the time series which changes periodically (i.e., a periodical function of time such that The SARIMA model of the time series {x(0), x(1), . . . , x(T)} with period s sample can be expressed as follows.
where d and D are respectively the order of local and seasonal integrated components, (t) is the Gaussian noise of which mean and variance are 0 and σ 2 , respectively, and B is a delay operator defined by Then, each components of the SARIMA model are given as follows. Moving-average: where p, q, P and Q are the orders of four components in (3)-(6) respectively.
The SARIMA model can express various periodical time series by designing the model parameters. In the followings, let us consider only the seasonal autoregressive (SAR) component to avoid the over fitting problem and to simplify the explanation of the prediction method. That is, let (1), then for this special case, we can obtain SAR model for the time series y(t) as follows.
where P is the order of the SAR components, Φ ρ , ρ = 1, 2, . . . , P are the SAR coefficients, and s are the constant period of the target time series.
Then, to substitute t + h for t, the prediction equation by using (7) can be given bŷ whereŷ(t + h|t) is the predicted value of future time t + h samples with h samples ahead of current time t, and the term −ρ × s + h must be not longer than 0 for composing the prediction by using only past observations.
The equations (7) and (8) show that an essential core of the general SAR depends on an assumption that each values of the same phase correlate each other. These equations can provide a simple description with a single model for the periodical variation with period s. For example, if target time series are given by the following deterministic and fully periodical function: where A n , n = 0, 1, 2, . . . and ϕ n , n = 1, 2, 3, . . . are amplitudes and initial phases for each n-th harmonic component, respectively. Then, let Φ ρ = 1/P in Eq. (8)  predict the future value of the periodical function x(t + h) at h-sample ahead future as follows.
x(t + h|t) Figure 9 shows a schematic diagram of SAR model-based prediction for simple periodic time series with period s. In this example, SAR model can predict the target value (drawn as star) by referring the past observed values (drawn as circles) which are similar to the target value by using constant period s.
As shown by the prediction example, the SAR model-based equation (8)  A simple solution for adapting the SAR model to the fluctuated periodicity is to substitute a time-variant period s(t) f orthe constant period s given bŷ Then the new prediction equation is written aŝ The equation above is about the same to the SAR model-based method [19] and it seems that it can adapt to the fluctuated periodicity, but this simple extension of SAR model has a limitation.
For explanation of the limitation, let us consider time series with a time-varying period as a function of time t. For this situation, the instantaneous phase ϑ(t) of the time series is given as a time-integration of the angular velocity ω(t), such as ϑ(t) = t 0 ω(t)dt. If the angular velocity is given as a linear function of time t, where ω a and ω 0 are constants. And the instantaneous phase is given a nonlinear (the second-order in this case) function as Then let us assume that s(t) is known or estimated correctly. That is, the term s(t) can refer the instantaneous phase same to the target time t, i.e., ϑ(t) = ϑ(t − s(t)) However, the terms ρ · s(t), ρ = 2, 3, . . . , P in (11) cannot refer the same instantaneous phase, i.e., ϑ(t) = ϑ(t − ρ · s(t)) + 2πρ, even if the s(t) is given as a true value.
Thus, the simple SAR model-based solution by using the time-varying period s(t) cannot express the fluctuated periodic nature suitably, even if the change of the periodicity is sufficiently simple such as a linear function of time. This is one of the limitations of SAR model in the prediction of fluctuated periodical time series such as respiratory tumor motion.

Time-variant SAR (TVSAR) model
To overcome the limitation of the general SAR model, a time-variant SAR (TVSAR) model for prediction of the lung tumor motion has been proposed [21]. The TVSAR equation corresponding to (7) can be expressed as where r ρ (t) are reference intervals of the ρth order at time t.
The reference intervals r ρ (t) can be defined as a time interval between the current time and the corresponding past time which has the same phase to the current one. In other words, if an instantaneous phase ϑ(t) of the time series is given, the reference intervals can be defined as follows.  Also, if the period is a constant, the reference interval can be expressed by Note that, the constant reference interval shown by (17) corresponds to that used in the general SARIMA equations of (7) and (8).
Then, an ideal prediction equation of the time-variant SAR is expressed by substituting t + h for t as follows.ŷ where the term h − r ρ (t + h) ≤ 0. Figure 11 shows a schematic diagram of TVSAR model-based prediction. The reference intervals can refer the past values observed which are corresponding to the target value, on the time series with fluctuated period. Also, Figure 12 shows TVSAR model-based prediction system.
Note that, in (18), we have to know the reference interval at h sample ahead future, but it is unknown in practice. In this case, we need to estimate it. The estimation method will be explained in the next section.

Online estimation of reference intervals: Short-time correlation analysis
In TVSAR model, the reference interval r ρ (t) is an important factor to predict the lung tumor motion accurately, and must be estimated on-line. In this study, a correlation analysis is adopted to estimate the reference interval.
The correlation analysis based estimation procedure of the reference interval is as follows.

Calculate a correlation function between the latest subset
{y(t − w), y(t − w + 1), . . . , y(t − 1), y(t)} (19) and k-sample lagged subset where w is a window length for subset time series. The correlation function is given by Here μ t and σ t are the sample mean and standard deviation of the subset time series. 2. The estimated reference intervalsr ρ (t|t) can be obtained by the intervals between k = 0 and the peak points of the correlation function γ(t, k) corresponded to each seasonal order ρ:r where l defines the search area for the peak of γ(t, k). 3. The window length is updated at each time by using the latest estimation of the first order reference interval as w = round(0.5r 1 (t|t)).
(23) Note that the windows length can affect the accuracy and response time to estimate the fluctuated periodicity. In this study, the window length was empirically set as almost half of single wave.
The initial condition is defined as follows.
wheres is the average of pre-observed time variant periods of the time series.s = 90 sample was used for this study.
An example of the correlation function and estimated reference intervals are shown in Figure  13. The reference intervals are thus estimated as intervals from lag 0 to ρ-th local maximum points of the correlation function.
Time-variation of the correlation function and estimated reference intervals are shown in Figure 14. As is clear from this figure, reference intervals of lung tumor motion intricately change with time evolution.
We can now estimate the reference intervals, but we need future value of them for (18). According to the figure 14, these reference intervals also fluctuate intricately and those prediction is difficult. So r ρ (t + h) cannot be directly used for the prediction in (18). As a realistic way, we simply extrapolated r ρ (t + h) from the current estimationr ρ (t|t) with zero-order hold:r ρ (t + h|t) =r ρ (t|t). Then, (18) can be rewritten aŝ

Experimental results
To evaluate the prediction performance of the methods described in Section 3, a prediction experiment by using clinical data sets was performed.

Clinical data sets of lung tumor motion
Three data sets of respiratory tumor motion time series were used for the experiment. All the data sets were observed with sampling frequency F s = 30 Hz and provided by Hokkaido University Hospital. Note that, a part of the data sets and measuring condition have been already shown in Section 2. Each data sets include three time series that correspond to spatial axes; left-right (LR), superior-inferior (SI), and anterior-posterior (AP) axes, respectively. The observational noises included in original data sets were preliminarily eliminated by using the statistical and low-pass filters. Figure 15 shows the three data sets and Table 1 summarizes data sets characteristics.

Tested methods
For comparison of the prediction performance, the following four methods were tested.
1. Zero-order hold (ZOH) method assumes that the latest position of tumor will not change in future. The prediction equation is given bŷ Note that the use of ZOH corresponds to the case that the system latency in radiotherapy machine is not compensated.
2. First-order hold (FOH) method assumes that the latest position and velocity will not change in future. The prediction equation is given bŷ 3. Seasonal autoregressive (SAR) model-based prediction given in (12). Note that the time-variant period is given as s(t + h) =r 1 (t|t).

Time-variant SAR (TVSAR) model-based prediction given in (26).
Parameters for SAR and TVSAR models are summarized in Table 2.

Performance indexes
To evaluate the prediction performance, mean absolute error (MAE) was calculated. MAE is given as a function of prediction horizon h as follows.
Here t s and t e are lower and upper bounds for the evaluation shown in Table 1, and e(t + h|t) is a prediction error defined by the Euclidean distance between the actual position and the predicted position.
The prediction error is calculated as follows.
In addition to MAE, prediction success rate was adopted to evaluate the prediction efficiency for short treatment time. While treatment, the irradiation must be interrupted when unacceptable distance between the irradiated field and the target volume is detected for patient safety. Then, frequent prediction failure may prolong the treatment fraction. Thus, the rate of prediction success can discover better method. The prediction success rate (PSR) is defined as follows.
PSR(h) = Number of successfully predicted samples Number of all evaluated samples (31) where the numerator is calculated as the number of the prediction errors within a threshold for tolerant accuracy. In this study, the threshold was set as 1 mm. Order P 2 Coefficients Φ ρ , ρ = 1, 2, . . . , P 1/P Maximum lag for correlation function 400  difference between SAR and TVSAR predictions at 87 to 89 s. Thus, TVSAR has provided better prediction. On the other hand, ZOH and FOH predictions are noisy and not fitted into the actual position. Their errors are frequently larger than 5 mm. Figures 17 and 18 show MAE and PSR as a function of prediction horizon h/F s (s). Each curves of prediction performances drawn in the figures are averaged over the data sets. Also, error bars indicates standard deviation of each performances for data sets variation.

Results
As is clear from the figures, TVSAR is superior to SAR constantly. This suggests that the basic concept of the TVSAR is more proper than SAR for the lung tumor motion prediction, i.e., reference intervals can work to improve the prediction accuracy for quasi-periodical nature. Then, MAE and PSR curves of TVSAR are the least and the highest for prediction horizon h/F s > 0.2 s, respectively. It means that TVSAR is the first best method in the tested methods for the radiotherapy machine with system latency of 0.2 s or longer. The average accuracy and 70 % of the predicted samples respectively are sub-millimeter. This may suggest that the TVSAR can perform tumor following irradiation during over 70 % of a single treatment fraction.
SAR is the second-best, but the average error for prediction horizon longer than 0.5 s is over 1 mm. Note that the difference from TVSAR depends on the use of the second reference interval in this experiment. Thus, unsuitable referring the past values by using the latest period, ρ · s(t), increases the considerable prediction error.
ZOH and FOH are superior to SAR and TVSAR predictions for shorter prediction horizons < 0.2 s. It is not remarkable because the position and velocity don't change drastically in short term in general. The periodicity in the respiratory motion quickly changes the position  and velocity at each peaks. Then ZOH and FOH cannot take into account the changes arisen. Therefore, their prediction errors drastically increase to unacceptable level with increase in prediction horizon. Those frequent and large prediction errors cause prolong the treatment fraction.

Conclusion
In this chapter, a prediction method of respiratory induced tumor motion for tumor following radiotherapy were introduced. The respiratory motion involves quasi-periodic nature as a dominant component and the motion causes unacceptable error of irradiation due to the system latency between tumor localization and beam-repositioning. To compensate the latency, prediction method of respiratory motion is very important and its prediction accuracy directly affects irradiation accuracy. Then, as a potent prediction method, time-variant seasonal autoregressive (TVSAR) model was introduced. TVSAR is an extended model of seasonal autoregressive model to be adaptable to the quasi-periodical nature. An experimental result by using clinical data sets of lung tumor motion showed the average accuracy of TVSAR is less than 1 mm for up to 1 s forward prediction and TVSAR is superior to other methods tested for <0.2 s ahead future prediction. Thus, TVSAR model-based prediction method achieved accurate prediction of the respiratory tumor motion and can help the tumor following irradiation.