Internal Motion Estimation by Internal-external Motion Modeling for Lung Cancer Radiotherapy

The aim of this study is to develop an internal-external correlation model for internal motion estimation for lung cancer radiotherapy. Deformation vector fields that characterize the internal-external motion are obtained by respectively registering the internal organ meshes and external surface meshes from the 4DCT images via a recently developed local topology preserved non-rigid point matching algorithm. A composite matrix is constructed by combing the estimated internal phasic DVFs with external phasic and directional DVFs. Principle component analysis is then applied to the composite matrix to extract principal motion characteristics, and generate model parameters to correlate the internal-external motion. The proposed model is evaluated on a 4D NURBS-based cardiac-torso (NCAT) synthetic phantom and 4DCT images from five lung cancer patients. For tumor tracking, the center of mass errors of the tracked tumor are 0.8(±0.5)mm/0.8(±0.4)mm for synthetic data, and 1.3(±1.0)mm/1.2(±1.2)mm for patient data in the intra-fraction/inter-fraction tracking, respectively. For lung tracking, the percent errors of the tracked contours are 0.06(±0.02)/0.07(±0.03) for synthetic data, and 0.06(±0.02)/0.06(±0.02) for patient data in the intra-fraction/inter-fraction tracking, respectively. The extensive validations have demonstrated the effectiveness and reliability of the proposed model in motion tracking for both the tumor and the lung in lung cancer radiotherapy.

Organs and tumors in the thoracic and abdominal region can move and deform significantly due to respiration 1 . The respiration-induced tumor motion can be up to 3 cm in the superior-inferior (SI) direction 2 and 2-4 mm in the anterior-posterior (AP) direction 3 . Respiration motion-induced translation, rotation, and deformation of the tumor and surrounding organs at risk (OARs) can cause significant geometric and dosimetric errors. Gierga et al. 4 reported that the planned target dose-volume histogram (DVH) was significantly degraded where the received CTV dose was reduced by 2-28% with tumor motion of 7.4 mm and 3.8 mm in the SI and AP directions, respectively. In a 4D Monte Carlo study, 3-5% dose differences (9.3 Gy tumor under-dosage) was observed between a 3DCT plan and a 4DCT plan where respiratory motion was considered 5 . Even with additional planning target volume (PTV) margin compensations for breathing motion, dose deviation of PTV D95 can be up to 26% for fractional dose and 14% for total dose with tumor motion 6 .
Respiration motion, on the other hand, can be utilized as an additional degree of freedom besides conventional 3D spatial domain to achieve 4D optimized treatment plan with greater OAR-sparing while maintaining PTV coverage and delivery efficiency 7 . The primary step to take advantage of respiration motion is to accurately track both the tumor and OARs motion. Many motion tracking strategies have been investigated in the past two decades. Some have been already successfully implemented in clinics [8][9][10][11][12] . Intuitively, directly imaging of tumor and OARs to obtain real-time positions is the most accurate tracking method [13][14][15][16] . However, real-time imaging techniques often have their own limitations. For instance, the X-ray fluoroscopy has low image contrast in soft tissue 15,16 . MRI can provide high-quality soft tissue images and is promising for tumor and OAR tracking 13,14 ; however, it is expensive and not widely available in clinics yet. While the clinically available MRI guidance radiotherapy units only have 2D planar tracking capability 14,[17][18][19] . Instead of tracking tumor and OAR directly, tracking

Methods and Materials
Ethics statement. This retrospective patient study was approved by Human Research Protection Program Office (HRPPO)/Institutional Review Board (IRB) of The University of Texas Southwestern Medical Center. All methods in this study were conducted in accordance with the relevant guidelines and regulations. Considering that this is not a therapeutical treatment study, our institutional review board waived the need for obtaining written informed consent from the participants.
Surface meshing. In this study, external and internal motions were represented by deformation of the external surface and internal organ surface, respectively. The surfaces used for modeling and validation were extracted from the 4DCT, on which the contours were first delineated and then the superior and lateral portion of the body contours, and the internal organ contours were converted to meshes using a particle-based surface meshing approach 39 . Given the contour points of the segmented organ masks (or particles in this algorithm), a high quality isotropic triangular surface meshing can be obtained by solving an inter-particle energy function with the quasi-Newton L-BFGS optimizer. The algorithm implementation details have been described by Zhong et al. 39 .
Motion and deformation tracking model. The internal and external DVFs were denoted as I j and S j to characterize the motion of internal organs and external surface on phase j of the planning 4DCT images. The I j and S j can be represented as Equations (1) and (2). I  I  I  I  I  I  I  [  ,  , ;

= …
Here, x, y, z represented three cartesian coordinates, M was the number of vertices of the internal organ surface, and N was the number of small patches on the external surface. Here, N = 154 patches (14 × 11, 14 in the SI direction and 11 in the ML direction) were uniformly extracted from the external surface via the strategy detailed in a previous work 33 . Both I j and S j were estimated by registering the surface points on phase j (j ∈ [1,10]) of the 4D planning CT to those on a middle position (MidP) CT, which is a time-averaged CT image representing mean position of patient's anatomy in the breathing cycle 40 . The DVFs (S n,x,j , S n,y,j , S n,z,j ) on patch n were calculated by averaging the DVFs on vertices inside the corresponding patch. Surface point registration was accomplished using a recent developed local topology preserved non-rigid point matching algorithm(TOP-DIR) 41 . Unlike internal motion, which was the combination effect of respiration, heart beating and organ deformation, the external surface motion is mainly caused by breathing. However, respiration induced motion may diminish when propagating from internal to external, and is weakly reflected on the external surface. To capture subtle external motion difference between different breathing process (e.g., from inhale to exhale vs. from exhale to inhale), we defined an A j to describe respiration-induced directional external motion between consecutive respiratory phases: A j was calculated by subtracting the DVFs S j−1 from DVFs S j . Vector d j was constructed by combining I j , S j and A j to describe the internal and external motion on phase j: The internal-external motion pattern for all phases was described by a composite matrix: The principle component analysis (PCA) was used to extract motion characteristics from D. Instead of calculating the eigenvalues and the corresponding eigenvectors of the covariance matrix DD T directly with a standard PCA procedure, we adopted the approach of Fayad et al. 37 to reduce computational expense. Let Xand λ be the eigenvectors and eigenvalues of matrix D T D, T multiplying D on both side of Equation (6) leads to: T Equation (7) indicates that DX and λ are the eigenvectors and corresponding eigenvalues of the covariance matrix DD T . Note that λ was the eigenvalues of both D T D and DD T , and computation of eigenvalues and corresponding eigenvectors from D T D was more efficient (D T Dwas with size J × J). Finally, the internal-external motion at time t can be approximated as a weight sum of the obtained eigenvectors E = [e 1 , …, e K ] of the largest K(K ≤ J − 1)eigenvalues as the Equation (8) where E I (size 3M × K) was constructed from the first 3M rows of E and E S (size 6N × K) was constructed from the rest 6N rows of E. Essentially,  I and  S were the internal and external DVFs, and E I and E S correspond to the internal and external components of the eigenvectors. By assuming E S was invertible and eliminating the unknown weight matrix W, the internal DVFs  I can be predicted using the external DVFs  S as: where B (size 3M × 6N) was a matrix correlating the internal and external motion.
Quantification of tracking accuracy. The accuracy of the proposed model was assessed by comparing the tracked internal organ contours with the manually delineated ground truths. Four similarity metrics were used including the center of mass (COM) error, the Dice's coefficient (DC) 42 , the percent error (PE) 43 and the Housdourf 's distance (HD) 44 . Given two volumes A (manual contoured volume/ground truth) and B(predicted volume), and their corresponding boundary . COM error was defined as the 3D Euclidean distance between the mass center of ground truth volume A and that of predicted volume B, to measure the tracking accuracy of the motion trajectory. DC, PE and HD were used to measure the agreement between the predicted tumor and organ contours and the manual delineated ground truths. The COM, DC and PE were defined as: In Equation (10), c(A) and c(B) represented the mass center coordinate of volume A and B and ||·|| was the L 2 norm. DC ranges from 0 to 1, corresponding to the worst and the best agreement, respectively. PE ranges from 0 to infinity, with 0 represents the best agreement. The HD was defined as: . The ground truth organ masks were contoured by an experienced physician. Better tracking results were indicated by lower COM, PE, and HD, or higher DC.
The developed model using the phasic DVFs and directional DVFs on the entire external surface (including all discrete surface patches, termed as SurMod) was compared with 1) using phasic DVFs (without directional DVFs) on partial external surface (10 selected patches, termed as RoiMod) 37 , and 2) using phasic DVFs (without directional DVFs) on the entire external surface (termed as SurphaMod).
Independent samples Kruskal-Wallis (K-W) test was adopted for organ tracking performance comparisons among RoiMod, SurphaMod, and SurMod. Inter-group comparisons were performed with Mann-Whitney tests. K-W tests and one-way analysis of variance (ANOVA) were conducted for non-parametric and parametric data, respectively, to analyze tracking performance difference. All statistical analyses were implemented using SPSS 19.0 software (SPSS Inc., Chicago, IL), and the statistical significance level was set at p = 0.05. For multiple comparisons, the p-value was adjusted accordingly using the Bonferroni correction method for individual comparison tests in SPSS.

Synthetic cases.
4DCT images with five breathing cycles (Cycles 1-5) simulated by using the 4D NCAT phantom 45 were employed for evaluation. As detailed in Table 1, the changes of breath period, amplitude and the atrophy of tumor response to the treatment were synthesized in the respiration motion in cycles 1-5. The tumor motions in the ML (Medio-Lateral), AP and SI directions for each respiratory cycle were given by: ML AP SI where H was the maximal amplitude of the surface motion in AP direction, T was the respiratory period, H ML , H AP and H SI were the motion trajectories of the tumor centroids in the ML, AP and SI directions, respectively. The resolution and voxel size of the synthetic cases were 256 × 256 × 120 and 2.0 mm × 2.0 mm × 2.5 mm.
In this NCAT phantom evaluation, Cycle 1 was used for motion modeling. The established correlation model was then applied for motion tracking, where Cycles 2 and 3 were used for intra-fraction validation (with different periods but same tumor diameter), and Cycles 4 and 5 were used for inter-fraction validation (with different periods and different tumor diameter).
Clinical cases. Clinical 4DCT images from five lung cancer patients (4 males and 1 female, ages range from 53 to 78 with median of 63) were collected for validations ( Table 2). All the 4DCT images were acquired on a Brilliance Big Bore-16 (Philips) CT scanner. The 4DCT images were sorted into 10 phases. Patients 1-3 have one set 4DCT image, and patients 4 and 5 have two sets 4DCT images. Both intra-4DCT and inter-4DCT evaluations were conducted. For intra-4DCT evaluation, the leave-one-out method was used, i.e., nine out of ten phases were  Data availability. The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Results
NCAT phantom. The quantitative comparisons, in terms of COM error, DC, PE and HD, between the RoiMod, SurphaMod, and SurMod are illustrated in Fig. 1. Only slight improvement is observed when comparing the SurphaMod with the RoiMod, while the proposed SurMod achieved the highest tracking accuracy for all the four metrics except for the HD in the lung tracking. The limited improvement of SurphaMod over RoiMod can be explained by the almost consistent surface motion pattern observed in the NCAT phantom, where employing the entire surface does not necessarily provide additional information in motion modeling. While further adding directional DVFs, as demonstrated by the SurMod, can offer more useful motion information in modeling the internal-external correlation. The advantage of utilizing complete surface and directional DVFs in motion modeling is visualized as an example case shown in Fig. 2, where the SurMod and the RoiMod are compared at the end of inspiration ( Fig. 2(a), (b)) and at the end of expiration (Fig. 2(c), (d)). It is observed that the SurMod can produce more accurate The RoiMod and SurphaMod both failed in tracking the tumor motion in the ML direction ( Fig. 3(a)), even though the simulated motion is relatively small (~3 mm in the ML direction). In contrast, the proposed SurMod is able to capture those subtle changes in the ML direction tracking than the RoiMod for both the tumor and the lung, especially in the apex of the lung (Fig. 2(a-d)-2). Figure 3 shows the intra-fraction and inter-fraction tracking comparisons in the ML, AP and SI directions as well as the COM trajectories and COM differences. The RoiMod and SurphaMod both fail in tracking the tumor motion in the ML direction ( Fig. 3(a)), even though the simulated motion is relatively small (~3 mm in the ML direction). In contrast, the proposed SurMod is able to capture those subtle changes in the ML direction, and yields superior tracking accuracy in all directions (Fig. 3(a) The quantitative evaluations on the synthetic phantom are illustrated in Table 3 and Fig. 4. For both intra-/ inter-fraction tracking, the SurphaMod achieves significant improvements over RoiMod in lung tracking accuracy, except for the HD of the intra-fraction tracking and the COM of the inter-fraction tracking. However, no significant improvement is observed in tumor tracking from RoiMod to SurphaMod. In contrast, the proposed SurMod achieves significant improvements over RoiMod in all metrics for both intra-/inter-fraction tumor and lung tracking. Furthermore, the SurMod also achieves significant improvements over the SurphaMod in tumor tracking except for the PE of inter-fraction tracking. After all, both the SurphaMod and SurMod achieve improvements over the RoiMod, and the SurMod performs the best. Quantitatively, the tracking accuracies achieved by the proposed SurMod in terms of COM, DC, PE and HD are 0. 8  We can see that the SurphaMod is able to produce more accurate tracking results than the RoiMod, though slight inferior COM errors are observed in the lung. Among the three models, the proposed SurMod yields the best results.
Visual comparisons of the intra-4DCT tracking of the clinical case 1 are shown in Fig. 6. More accurate tumor and lung contours are predicted by the SurMod compared with the RoiMod, indicates by better agreements between the tracked contours with the physician delineated ground truths. However, the COM error of SurMod is larger than RoiMod and SurphMod in phase 50% and 60%. The reason for this is that the COM error is a biased metric when the target involves large deformation in motion, especially in scenarios(such as the lung at the end of exhale). The four metrics combined results show the three model behave similarly in phase 50% and 60%.
The tracking accuracy comparisons for all clinical cases (Cases 1 to 11) are shown in Table 4 and Fig. 8. In the intra-/inter-4DCT tracking, significant improvements are observed from SurphaMod over RoiMod in COM error of intra-/inter-4DCT tracking and HD of the inter-4DCT tracking for tumor, and in all metrics of intra-/inter-4DCT lung tracking except for the COM error of intra-4DCT tracking and the HD of inter-4DCT tracking. The SurMod also achieves significantly better tracking accuracies over RoiMod except for the COM Efficiency. All experiments in this study were conducted on a CPU platform equipped with 8 GB memory, and the proposed model was coded and implemented on the software platform of MatlabR2011a. The average computational time for motion modeling and prediction were ~9 minutes and ~13 seconds. The computational time was closely related to the number of the generated mesh vertices of the internal organs and external surface. Intuitively, more vertices can depict more anatomical details on organ surface, however, it does not necessarily imply higher registration accuracy 41 . In this study, therefore, the quantity of mesh vertices was empirically set as 1500 to balance the accuracy and efficiency, which was proved to be adequate to yield satisfactory registration result.

Discussion
As the evaluation results on synthetic NCAT phantom data and clinical data illustrated, the proposed SurMod achieved superior performance over the RoiMod and SurphaMod in the internal organs tracking for lung cancer radiotherapy. It also can be observed that, compared with the performance in clinical cases, the proposed model achieved better in the NCAT phantom. The reasons are two folds: firstly, the breathing pattern is constant in the NCAT phantom but might be irregular in real patient data. The proposed model is established using one set of the 4DCT image with the assumption that change in breathing pattern is small; however, this is an ideal hypothesis which might not happen in more complicated clinical cases. Secondly, since the target shape may vary considerably between different treatment fractions, thus, delineation consistency of the target will introduce another source of error for tracking accuracy evaluation. Compared with the inter-fraction tracking, higher accuracy was seen in intra-fraction tracking. This may attribute to the large anatomical changes between different treatment fractions, even though the model had been adapted by registering the patient anatomy between fractions, the anatomical differences were unlikely to eliminate. Furthermore, the respiration pattern tends to be stable in intra-fraction treatment. However, it may vary considerably between treatment fractions. With effective model update strategies, this tracking accuracy discrepancy between intra-fraction and inter-fraction should be reduced.
The effectiveness of the proposed SurMod attributes to the utilization of motion information from the entire surface, as well as directional respiration-induced DVFs A j . The partial ROIs based model was proved to be superior to the conventional approaches using phasic, amplitude, and phasic & amplitude as external surrogates 37 .  Table 3. Quantitative comparisons (Mean ± STD(p-value)) between the RoiMod, SurphaMod and the SurMod in intra-fraction tracking and inter-fraction tracking for the NCAT phantom. However, the location of the ten proposed ROIs were identified from ten patients, which might not be able to generalize to a broader patient cohort because of inter-patient respiration variations. The generally superior performance of the SurphaMod over the RoiMod implied the robustness and effectiveness of using entire surface information. Moreover, benefited from the utilization of directional DVFs A j , the proposed SurMod outperformed SurphaMod in scenarios when large inspiration-expiration differences occurred (e.g., synthetic cases, case 1 and case 9, as results shown in Figs 1, 5 and 7, respectively). In general, compared with the RoiMod 37 , the proposed SurMod improved the mean DC from 0.84/0.95 to 0.89/0.97 for inter-fraction tumor/lung tracking  For the proposed SurMod, the modeling procedure costs <10 mins in a CPU platform and the motion prediction costs <15 s, where most of the computation time is spent on surface meshes registration. Though 8~10 frames per second (fps) tracking speed 46 is usually expected in a clinical context for real-time tracking, 10 × or more acceleration, as reported by previous studies [47][48][49][50] , is possible and not effort-demanding. This can be achieved by simply coding the current model in parallel in a GPU environment equipped with high-end computational platform.
Though promising results were achieved on synthetic phantom data and clinical patient data, there are many challenges needed to be addressed before applying the proposed model in a clinical setting. Firstly, the proposed model is based on the stable respiration assumption, which may not apply to all patients, especially for those late stage lung cancer patients with poor lung function. Model updating is necessary to adapt the model to a new on-treatment scenario 51,52 , especially for inter-fraction tracking. Therefore, model validation during treatment and a regular model update scheme will be essential for renewing the model established on the pre-treatment 4DCT. The model update can be achieved by registering the patient anatomy on the treatment day with that from the planning 4DCT, which is also used to accommodate the inter-fractional baseline drift in this study. The efficiency of the model update is closely related to the registration algorithm used. The registration process is expected to be done within half a minute with the aid of GPU acceleration 53 . With current technologies, the kV X-ray verification is the most practical tool for the tumor tracking model validation, which is also easy for model update to accommodate the baseline shift estimation. Since the purpose of the proposed model was to track both tumor location and tumor shape, projection images from multiple directions might be required to provide benchmark comparisons of tumor shape in different projection angles.
Furthermore, more complicated circumstances (e.g. phase shifts, baseline drifts and hysteresis, etc.) need to be evaluated for the proposed model. In this study, to evaluate the model's responses to general respiration motion, primary parameters changes (period, amplitude and tumor shape) were simulated in the synthesized motion in the NCAT phantom study. We did not synthesize more complicated motion, where the reasons are two folds: (1) there are too many possibilities of respiration pattern to be evaluated if more parameters are added in motion simulation. Even if all these possibilities were enumerated, it is still difficult, if not impossible, to realistically synthesize respiration since breathing pattern is far more complicated and diversified in patients; (2) the tracking accuracies responding to above practical issues are closely related to the corresponding solution scheme.  Table 4. Quantitative comparisons (Mean ± STD(p-value)) between the RoiMod, SurphaMod and the SurMod in intra-fraction tracking and inter-fraction tracking for all the clinical cases. For example, aligning the patient anatomy on the treatment day with that defined on modeling is an alternative solution to address inter-fractional baseline drift. The influence in tracking accuracy from the baseline drift is directly related to the registration uncertainties of the adopted registration algorithm, which is beyond the scope of the proposed mathematical model. However, more comprehensive model evaluation and adaption need to be addressed in future studies before successful clinical applications. Finally, in clinical practice, the external surface is obtained by optical surface imaging devices (e.g. ToF camera and AlighRT etc.), which has certain practical issues such as, image acquisition latency. Latency is a common issue in many commercial surface monitoring systems, and model forward prediction might be a practical way to accommodate this issue. For the proposed model, which was built on external DVFs estimated from DIR of external surface mesh, a forward predicted surface mesh might be needed. As an alternative, the external DVFs can be extrapolated using known phasic external DVFs, or breathing velocity/accelerated velocity calculated from directional external DVFs.

Conclusion
In this work, we proposed and validated a novel motion-tracking model in lung cancer radiotherapy, which is constructed based on the correlation between the phasic DVFs on the internal organs' surface and phasic and inter-phasic DVFs on the patient's external surface. Experimental evaluations were conducted on a synthetic 4D NCAT phantom and 4DCT images from five lung cancer patients. The evaluation results have demonstrated the capability of the proposed model in internal organs tracking with the variance of breath frequency, amplitude, inter-fraction anatomical baseline shift and deformation. The good performances on the evaluation data demonstrated the effectiveness, reliability and accuracy of the proposed model. However, more comprehensive model evaluation and model adaption are needed before adopting the proposed model in a clinical setting.