Analysis of Deformation Fields in Spatio-temporal CBCT images of lungs for radiotherapy patients
Bijju Kranthi Veduruparthi Jayanta Mukherjee Partha Pratim Das Mandira Saha Raj Kumar Shrimali Sanjoy Chatterjee Soumendranath Ray Sriram Prasath
stat.AP









bijjuair@gmail.com


Department of Computer Science and Engineering, Indian Institute of Technology, Kharagpur









Department of Radiation Oncology, Tata Medical Center, Kolkata


Department of Nuclear Medicine, Tata Medical Center, Kolkata


Department of Medical Physics, Tata Medical Center, Kolkata

[corra]Corresponding author



Deformable registration of spatiotemporal Cone-Beam Computed Tomography (CBCT) images taken sequentially during the radiation treatment course yields a deformation field for a pair of images.
The Jacobian of this field at any voxel provides a measure of the expansion or contraction of a unit volume.
We analyze the Jacobian at different sections of the tumor volumes obtained from delineation done by radiation oncologists for lung cancer patients.
The delineations across the temporal sequence are compared post registration to compute tumor areas namely, unchanged (U), newly grown (G), and reduced (R) that have undergone changes.
These three regions of the tumor are considered for statistical analysis. In addition, statistics of non-tumor (N) regions are taken into consideration.
Sequential CBCT images of 29 patients were used in studying the distribution of Jacobian in these four different regions, along with a test set of 16 patients.
Statistical tests performed over the dataset consisting of first three weeks of treatment suggest that, means of the Jacobian in the regions follow a particular order. 
Although, this observation is apparent when applied to the distribution over the whole population, it is found that the ordering deviates for many individual cases.
We propose a hypothesis to classify patients who have had partial response (PR). Early prediction of the response was studied using only three weeks of data.
The early prediction of response of treatment was supported by a Fisher's test with odds ratio of 5.13 and a p-value of 0.043.



Deformable Image Registration Jacobian Cone-Beam Computed Tomography(CBCT) Statistical analysis








§ INTRODUCTION

Lung cancer patients undergo a six week fractionated course of radiotherapy as part of  their curative treatment. Radiation therapy planning requires identifying the tumor based on a pre-treatment planning computed tomography(CT) taken in treatment position followed by localization of the tumor by radiation oncologists. During the process nearby “organs at risk” are also identified to ensure that radiation induced damage is avoided. This localized volume is then planned to receive a specified (often 60Gy) dose of radiation in thirty fractions over a period of six weeks. Modern radiation therapy machines can acquire a cone beam CT scan(CBCT) image on a daily or more commonly, a weekly basis prior to delivery of the radiation to negate any positional changes in the identified tumor. Such image guided radiation therapy is often the standard of care for the curative therapy of Lung cancers. Adaptive Radiation Therapy (ART) requires oncologists to adapt the radiation volume based on spatio-temporal changes in the tumor, and has the potential to personalize radiation therapy based on such changes. Registration of such images therefore remain a key area of concern and could not only help identify the target accurately, but could also potentially prognosticate the outcomes in such patients.

From a big dataset of radically treated lung cancer patients in a tertiary care center, we randomly selected their pre-treatment CBCT and on-treatment serial CBCT images. Weekly CBCTs of the patients are delineated by experienced radiologists for approval, and treatment. Each of the CBCTs are first delineated by an expert radiologist, followed by correction and approval of a second senior radiologist. Finally, the approved delineations are used by the radiation planning team to position the patient for radiotherapy.

The work studies the statistical analysis of deformation fields obtained by registering pre-treatment CBCT and serial on-treatment CBCT images of lung cancer patients who underwent radiotherapy. Using only the first three weeks of treatment as data, we aim to predict the response of the patient to treatment before the therapy completes. To ensure clinical relevance, we used the actual target volumes delineated by the clinical oncologist on pre-treatment CBCT for treatment and thereafter on CBCT to assess the changes in the targets. The delineations were used to categorize the tumor into different regions. We study the Jacobian statistic on the deformation fields in these regions and study how it varies across patients. We observe that patients satisfying the proposed hypothesis have better response to treatment than those who do not satisfy the hypothesis. Using this hypothesis, we are able to classify patients into PR and non-PR categories, and support the hypothesis using Fisher's association test.



 §.§ Related work


Computational Anatomy <cit.> is a study of anatomy using mathematical and computational models. The deformable template that exhibits diffeomorphism, a bijection between anatomical coordinates is sometimes used in such a study. In <cit.>, landmarks around the area of interest were used to track changes, where the Jacobian maps showed the direction of growth in ventral and dorsal parts of the hippocampus of a mouse. Deformation fields obtained after deformable image registration were earlier used in <cit.> for characterizing different regions of the lung, based on the motion properties like directional change, volume change and nature of change. In another study <cit.>, the deformation fields were used to track changes in tissue volumes. Different regions of the lung were segregated into blocks, and average Jacobian in these regions was observed to report lung activity and for tracking change in volume of tissues.
In <cit.>, anatomical changes in the brain were observed to be different across spatial scales using the Helmholtz decomposition of the deformation field. A difference of Gaussians (DOG) operator was applied on the irrotational component of the decomposition to identify the areas of maximal volume change. However, this method was applied only to brain MR images. 
In another work <cit.>, an average template was generated using registration of brain MR images. Statistical analysis on the deformation tensors with respect to this average template was performed to find anomalies in brain. Longitudinal analysis of spatiotemporal data to measure the growth of hippocampus in the brain is reported in <cit.>. Here, differences between growth trajectories are used to estimate a mean growth behavior of a population. This is compared with subjects to identify delay in growth for children diagnosed with autism.
Biological growth using deformation fields is mathematical modeled in <cit.>. The model is called a Growth by Random Iterated Diffeomorphisms (GRID). This model uses random seeds with radial deformations around it, to capture growth based deformations. In <cit.> and <cit.>, the growth variables in GRID model are directly captured from image data. All the aforementioned approaches use Jacobian for measuring growth and decay from images using nonlinear registration. 



 §.§ Contribution

This work presents a statistical analysis of deformation fields obtained by registering CBCT images of lung cancer patients who underwent radiotherapy. In our work, we also make use of delineation of tumor in CBCT images by radiologists. This is the first work of its kind on a dataset annotated by radiologists, and studying the behavior across population. The delineations were used to categorize the tumor into different regions. We study the Jacobian statistic on the deformation fields in these regions, and study how it varies across patients. Patients satisfying our proposed hypothesis were found to have better response to treatment. The results of the classification suggest that early prediction is feasible.




§ METHOD



 §.§ Image Registration

Given two images S, T ∈ℝ^3, the goal of image registration is to find a transformation
g⃗ : ℝ^3 ↦ℝ^3 that maps/aligns S onto T.
The approach of computing both the forward (registering S to T) and reverse (registering T to S) transformations is termed as bidirectional symmetric registration. This process involves computation of inverse of the deformation field.

A symmetric log-domain based nonlinear image registration technique was proposed in <cit.>. The deformation field obtained using this method is diffeomorphic and the true inverse can be computed efficiently at a very low cost. Hence this method is suitable for computational anatomy. This technique, however, is based on the sum of squared differences (SSD) in the intensities of the image. This technique is not quite applicable in our scenario, on account of the high amount of noise present in CBCT images. The noise in CBCT images is due to several artifacts arising from usage of low energy beams to produce the image. Further improvement in the results of the registration technique is reported  in <cit.>, where the symmetric Local Correlation Coefficient (LCC) provided a robust similarity measure. The technique was found to give smooth deformation fields owing to the regularization imposed on the total energy. We use this technique to perform non-linear image registration.


The deformation field that warps S to T is
denoted as ϕ⃗ = z⃗ - g⃗(z⃗), where z⃗ is the set of points in S that are mapped to corresponding points z⃗ - g⃗(z⃗) in T with displacement field g⃗(z⃗). Such a deformation ϕ⃗ yields proper alignment of the two images.
Trilinear interpolation was used to compute the warping of deformation field at non-integer coordinates.


The deformation field obtained using this method is diffeomorphic and the inverse can be computed efficiently at a very low cost. Moreover, this method guarantees invertibility of the deformation and is invertible. Hence this method is suitable for computational anatomy.



 §.§ Analysis of deformation fields


The deformation index for interpreting the information in a deformation field is the determinant of Jacobian, commonly referred in the literature as simply Jacobian.
The Jacobian of a deformation field ϕ is defined as:

    J(ϕ⃗(z⃗)) = [ ∂ϕ_1(z⃗)/∂ z_1 ∂ϕ_1(z⃗)/∂ z_2 ∂ϕ_1(z⃗)/∂ z_3               ; ∂ϕ_2(z⃗)/∂ z_1 ∂ϕ_2(z⃗)/∂ z_2 ∂ϕ_2(z⃗)/∂ z_3               ; ∂ϕ_3(z⃗)/∂ z_1 ∂ϕ_3(z⃗)/∂ z_2 ∂ϕ_3(z⃗)/∂ z_3 ]


The determinant of J(ϕ⃗(z⃗)), denoted as |J(ϕ⃗(z⃗))| or in short J , can also be expressed in terms of its eigen values. The Jacobian (J) measures the local volume change with respect to that of a unit cube.
This can be visualized in Fig. <ref>.
J varies between 0 and ∞. When J = 1, it means there is no change in volume.
A value of J < 1 denotes, net contraction; and J > 1 denotes, net expansion.




We computed the Jacobian in different areas of the 3D image and calculated the variation in the Jacobian statistic. The categorization can be visualized in Fig. <ref>. The Jacobian is computed individually in the following regions:


  * Deformed Tumor Region in week i (T_i): Set of tumor voxels when the image is deformed.

  * Tumor Region in week i (T^'_i): Set of tumor voxels when the image is undeformed.

  * Deformed Non-Tumor Region in week i (N_i): Set of non-tumor voxels computed as negation of the set T_i.

  * Non-Tumor Region in week i (N^'_i): Set of non-tumor voxels computed as negation of the set T^'_i.

  * Unchanged Region in week i (U_i): The set of tumor voxels that is the intersection of T_i and T_i+1 i.e, U_i = T_i ∩ T^'_i+1.

  * Reduced Region in week i (R_i): The set of tumor voxels that is defined as, R_i = T_i ∖ T^'_i+1.

  * Newly grown Region in week i (G_i): The set of tumor voxels that is defined as, G_i = T^'_i+1∖ T_i.


Applying the discussed model on real images, several slices are shown in Fig. <ref>. We study the deformation fields in these regions.

































 §.§ Distribution of Jacobian

Statistical tests on the Jacobian are performed for the entire dataset on observing the distribution of J in regions of different categories as described in Section. <ref>.
We assume normality of the dataset in all the regions due to the very large number of samples. The number of samples in R, G, U, and N regions are approximately, 1.8, 1.5, 3.5 and 169 million, respectively.

The confidence intervals around the mean can be estimated as [x̅± 1.96 σ/√n], for a confidence of 95%. We also computed the confidence intervals for the different regions using bootstrap method. In bootstrapping, the samples are resampled to compute the confidence intervals. The bootstrap confidence intervals and the mean for a population of 29 patients are shown in the Table. <ref>. It can be observed that the intervals are very narrow.



 §.§ Two-sided t-test

The two-sided t-test was performed under the null hypothesis that two independent distributions have identical expected values. The t-test was done with an assumption that the two populations have unknown identical variances. The test gives a p-value that explains the level of significance achieved. The p-value gives the probability of achieving a result equal to or more extreme than what was observed, assuming the null hypothesis is true. A very low p-value leads to the rejection of null hypothesis. The p-value was found to be close to zero for all the pair of regions. The t-statistic reveals which of the two samples have higher or lower expected values.

Samples were collected from the different regions of the tumor and the test was conducted on the population of 29 patients.
Let us denote the mean(μ) and standard deviation(σ) in the regions with the respective subscripts, where T corresponds to tumor regions, N corresponds to non-tumor regions, U corresponds to unchanged regions, R corresponds to reduced regions, and G corresponds to newly grown regions, respectively.


The t-statistic values are interpreted such that, when the alternative hypothesis is μ_X ≥μ_Y, then large positive values of t lead to rejection of the null hypothesis, where X and Y are defined as any of the four regions i.e, X,Y ∈R,U,G,N| X ≠ Y. Alternatively, when the t-statistic values are large and negative, then μ_X ≤μ_Y. 
Fig. <ref> shows the box-plot for the population in different regions.




The t-statistic values are shown in Table. <ref>. Here, X corresponds to the first column and Y corresponds to the first row. It can be seen that μ_R ≤μ_G, μ_R ≤μ_U, and μ_G ≤μ_U, as the t-statistic values are large and negative.
Similarly, μ_R ≥μ_N, μ_G ≥μ_N, and μ_U ≥μ_N, respectively, due to large positive values of the t-statistic. 
By arranging the t-statistic values on a real line, we observe that μ_N ≤μ_R ≤μ_G ≤μ_U for the population. These results corroborate with the intuition also, as the R region must be more aggressive than the N region. R region has recently reduced to non-tumor. Further, the G region has recently manifested as tumor, hence we expect it to be aggressive; however, lower than the U region which has currently not undergone any change. The above observation also is in conformation with the fact that G and U are relatively more aggressive than R and N regions. 



 §.§ Ordering of Jacobian for individual patients

Radiology Response Evaluation Criteria in Solid Tumors(RECIST) <cit.> is used to assess the imaging of each patient and analyze response to radiation treatment.
Patients who have undergone treatment are categorized according to this criteria according to their response. The response of a patient is checked by radiation oncologists typically three months after completion of radiotherapy. A follow up CT scan is performed to check for progression. Based on the tumor size and metastasis, the patient is classified into one the categories as follows:

  
  * Complete Response (CR): Disappearance of all target lesions.
  
  * Partial Response (PR): At least a 30% decrease in the sum of the longest diameter(LD) of target lesions, taking as reference the baseline sum LD.
  
  * Stable Disease (SD): Neither sufficient shrinkage to qualify for PR nor sufficient increase to qualify for PD, taking as reference the smallest sum LD since the treatment started (% change between 30% decrease and 20% increase).
  
  * Progressive Disease (PD): At least a 20% increase in the sum of the LD of target lesions, taking as reference the smallest sum LD recorded since the treatment started or the appearance of one or more new lesions.
  
  * Distant Progression (DP): The cancer is spreading from the original (primary) tumor in the prostate to lymph nodes or distant organs such as the bones, liver and lungs.







  §.§.§ Ordering hypothesis

We perform the t-test for each of the patients to record the ordering of the means in different regions.


The hypothesis to classify each patient is explained in Algorithm. <ref>.
The aim is to correctly classify partial response(PR) patients using the hypothesis.
It is to be noted that, by definition, a CR qualifies as a PR but a PR does not qualify as a CR.
Therefore, in our analysis PR category of classification refers to patients with response as either PR or CR, while the non-PR category refers to patients with response neither as PR nor as CR.
We perform the t-test for each patient in the population and test set.
The ordering from t-test is modified in the hypothesis for classifying patients into PR and non-PR categories. It is designed based on the following logic: By definition, we know that the expansion in G and U should be greater than that of R, because of more tumorous activity in these regions. Therefore, we hypothesize that μR≤{μ_G,μ_U}. This hypothesis does not consider the ordering between G and U. Additionally we add another condition based on the tumor activity in the R region. Intuitively, since R region should result in net reduction of the volume, the mean of Jacobian μ_R should be less than 1.0.
It can be observed from Fig. <ref>, that the mean of patients for PR is less than 1.0, while that of NPR is close to 1.0.





Fig. <ref> shows the block diagram of the classification process.
Our data consists of the radiation response of each patient in terms of the RECIST criteria. Table. <ref> shows the corresponding radiation response(RX_Response) for each patient. Patients whose response to radiation could not be determined by doctors are labeled as “NA". 

Based on this hypothesis, we label each patient in the “Classification" column of Table. <ref> in  <ref>. We compare the predicted and the actual response(RX_Response) to summarize into Table. <ref> and Table. <ref>, while excluding the patients whose response was “NA".
Overall, we find that 12 out of 21 patients are correctly classified as PR, and 13 out of 17 are correctly classified as non-PR patients.






Fisher's exact test for association between PR and the hypothesis was performed on the contingency table Table. <ref> and Table. <ref>, which yields an odds ratio(OR) and a p-value, to measure the association between two categories PR and non-PR. The null hypothesis is that the occurrence of PR and non-PR are equally likely. The odds ratio(OR) measures the ratio of odds of occurrence of an event to the odds of an event not occurring. The contingency table Table. <ref> corresponding to complete six weeks data yielded an OR and p-value of 4.33 and 0.051, respectively.

Using only the first three weeks data, we find that 11 out of 21 patients are correctly classified as PR, and 14 out of 17 are correctly classified as non-PR patients. From Table. <ref>, for the first three weeks data, the obtained result for OR and p-value were 5.13 and 0.043, respectively. Assuming a critical p-value of 0.10 as threshold, we can reject the null hypothesis that there is similarity in occurrence between PR and non-PR categories. The OR of 5.13 represents that it is 5.13 times more likely that PR happens when the hypothesis is satisfied than non-PR when the hypothesis is satisfied. This analysis suggests that the hypothesis is satisfied for patients with response as PR during the first three weeks of the treatment also. This early prediction can help doctors take necessary actions for improving the response of the patient to treatment. 

The accuracy, precision and recall were found to nearly the same on the three weeks data as shown in Table. <ref>. The level of significance is close to being significant considering the full six weeks and first three weeks of the dataset.





§ CONCLUSION

In this work we discussed analysis of Jacobian obtained from the deformation fields of registering spatiotemporal data of CBCT images for patients who underwent radiotherapy. Using clinical radiation oncologists delineation over the first three weeks of treatment duration, we analyzed the behavior of Jacobian in four regions of each delineated tumor. We observed very narrow variation of the Jacobian in the population in each of the regions. Two-sided t-test was performed to identify the ordering of mean Jacobian for the population in the four regions as described above. Based on the ordering obtained from the population, we proposed a hypothesis for classification of each patient within the population and test data into PR and non-PR classes. This hypothesis was used to segregate patients into those satisfying the hypothesis and not satisfying the hypothesis. We observed that patients satisfying the proposed hypothesis had better RECIST response to radiation treatment. Significant association between the proposed hypothesis and better response was confirmed by using Fisher's test. 
This early prediction has significant impact in providing radiologists necessary feedback regarding the response of patient to the treatment. Early prediction before the stipulated radiation course of 6 weeks gives enough time for doctors to take alternate actions like increased radiation dose, or alternative treatment approaches like chemotherapy, surgery. Using this result, the images of new patients who have undergone treatment can be processed using the discussed work flow, to qualitatively predict the response of a treatment. Further analysis can be performed to predict quantitative measures that can indicate prognosis for the patient.



§ ACKNOWLEDGMENT

This work is carried out under the MHRD sponsored project entitled as
“Predicting Cancer Treatment outcomes of lung and colo-rectal cancer by modeling and analysis of anatomic and metabolic images". We would additionally like to thank, Partha Sen and Gaurav Goswami for helping us procure the data and validating the image registrations.



§ REFERENCES






§ DETAILS OF THE DATASET





