Computer Simulation and Analysis of Three-Dimensional Tumor Geometry in Radiotherapy

Radiotherapy plays an important role in the treatment for patients with solid tumors. Recently, the advantages in high-precision radiotherapy enable focusing of higher radiation energy (dose) to the tumor region while minimizing unwanted radiation exposure to surrounding normal tissue to avoid radiation injury. Intensity-modulated radiation therapy (IMRT) varies the intensities and profiles of beams from various directions to fit the tumor size and shape. This technique greatly improves dose concentration on target region and normal tissue sparing. Image-guided radiation therapy (IGRT) uses advanced imaging technology such as on-board imaging system to achieve precise and accurate dose delivery. Many studies have reported inter-fractional organ motions and efficacy of IGRT in reducing targeting errors using daily CT images (Den et al. 2009, Wang et al. 2009, Houghton et al. 2009, Pawlowski et al. 2010, Varadhan et al. 2009, Greene et al. 2009). Owing to these techniques, errors in patient set-up and dose delivery can be reduced to some extent. However, as radiotherapy typically takes several weeks, tumor and normal tissues may deform due to therapeutic effect or loss of body weight during treatment period. Shapes and locations of the tumor and the surrounding organs would be quite different from those when the treatment was planned. This results in overdosage of surrounding normal tissue or underdosage of target region. To overcome this issue, it would be useful to precisely analyse and predict the changes in three-dimensional (3D) geometries of tumors and normal tissues through the treatment.


Introduction
Radiotherapy plays an important role in the treatment for patients with solid tumors. Recently, the advantages in high-precision radiotherapy enable focusing of higher radiation energy (dose) to the tumor region while minimizing unwanted radiation exposure to surrounding normal tissue to avoid radiation injury. Intensity-modulated radiation therapy (IMRT) varies the intensities and profiles of beams from various directions to fit the tumor size and shape. This technique greatly improves dose concentration on target region and normal tissue sparing. Image-guided radiation therapy (IGRT) uses advanced imaging technology such as on-board imaging system to achieve precise and accurate dose delivery. Many studies have reported inter-fractional organ motions and efficacy of IGRT in reducing targeting errors using daily CT images (Den et al. 2009, Wang et al. 2009, Houghton et al. 2009, Pawlowski et al. 2010, Varadhan et al. 2009, Greene et al. 2009). Owing to these techniques, errors in patient set-up and dose delivery can be reduced to some extent. However, as radiotherapy typically takes several weeks, tumor and normal tissues may deform due to therapeutic effect or loss of body weight during treatment period. Shapes and locations of the tumor and the surrounding organs would be quite different from those when the treatment was planned. This results in overdosage of surrounding normal tissue or underdosage of target region. To overcome this issue, it would be useful to precisely analyse and predict the changes in three-dimensional (3D) geometries of tumors and normal tissues through the treatment.
However, even methodology to evaluate 3D tumor shapes has not been established yet. At present, tumor diameter is commonly used as an indicator to evaluate therapeutic response in cancer patients. Since the 1970s, the World Health Organization (WHO) had suggested to assess tumor response by measurement of maximum diameter and largest perpendicular diameter (World Health Organization, 1979). The response evaluation criteria in solid tumor (RECIST) guideline proposed to use only maximum diameter in categorizing therapeutic response (Therasse et al. 2000, Werner-Wasik et al. 2001. Some researchers have reported that this difference in diagnostic criteria often resulted in different categorization of therapeutic response. Rohbe et al. stated that volumetric measurement with CT might help to evaluate therapeutic response more accurately (Rohbe et al. 2007). For more precise assessment, it would be useful to examine 3D tumor morphological features.
For the purpose of prediction of therapeutic effect, various computational methods have been proposed for simulation of tumor growth and shrinkage as the effects of radiotherapy. Almost all of these method consider tumor as concentration of huge number of cancer cells and calculate the death or birth of these cells pursuant to some rules based on the cell biology (Dionysiou et al. 2004, Stamatakos et al. 2010, Kolokotroni et al. 2011. These approaches look attractive, however it seems to be difficult to implement them into clinical situation because therapeutic effect cannot be evaluated accurately only by the number of cancer cells. Usable method for prediction and evaluation of therapeutic response of tumors has been required strongly.
Authors have proposed novel numerical simulation method to predict changes in tumor volume in radiotherapy (Takao et al. 2009a). We have also established the methodology to evaluate 3D tumor shape visually using color map called surface geometry map (Takao et al. 2009b and. From these knowledge, we proposed simulation method to predict changes in 3D tumor shape through the treatment duration and evaluated tumor geometry visually and quantitatively.

Analysis of three-dimensional tumor geometry
It would be valuable if changes in 3D tumor shapes during treatment can be evaluated visually and quantitatively. Quantitative evaluation could be performed using geometric factors such as volume, surface area, and radius. For visual evaluation, we introduce some color map called surface geometry map to represent 3D tumor geometry (shape) into twodimensional plane. The volume of tumor is also represented as size of the map. Therefore, this surface geometry map is useful to understand how the tumor shrinks during treatment.

Geometric factor
For quantitative (or numerical) analysis of tumor geometry, volume V, surface area S were measured and evaluated based on the 3D surface models of the tumors constructed from CT image sets taken once a week during treatment. From the values of V and S, spherical shape factor (SSF) which represents degree of sphericity was calculated by following equation (Choia HJ and Choi HK 2007).
If SSF=1, the tumor has a perfect sphere. As surface roughness increase, the value of SSF decreases. Tumor radius R(θ, φ) was also measured in all directions from the gravitational canter of the tumor.

Surface geometry map
We introduced surface geometry map like a global map to visualyze changes in 3D tumor shape (Takao et al. 2009a). Distances from the center of the gravity to its surface, i.e. radius R(, φ) were measured in the spherical coordinate system O-Rφ, with the origin is set at the tumor center. The angle  represents the azimuthal angle (-180° ≤  ≤ 180°), and the angle φ www.intechopen.com represents the polar angle (-90° ≤ φ ≤ 90°). Radius R(θ, φ) was sampled at 10° intervals in both the polar and azimuthal directions. To enable a visual understanding of the features of the 3D geometry of the tumor, the values of radius R(θ, φ) were represented in a color scale and plotted in the - plane: red indicates the maximum radius and blue the minimum within the tumor. The warm colors (red and yellow) represent convex region and cool colors (blue) represents concave areas (Fig. 1).
Further, to evaluate changes in tumor geometry quantitatively, the image correlation was analyzed. Surface geometry maps were converted into grayscale, therefore the intensities in the maps represented the tumor radius. The intensities at every position in the maps created from CT images taken during the treatment were compared with those of the corresponding position in the map before the treatment, and correlation coefficient was calculated. Similar analysis was performed to evaluate our simulation method, especially calculation considering tumor heterogeneity (discussed later).

Equations for tumor deformation
Equations for describing changes in 3D tumor geometry were established using analogy with deformation of continuous body in solid mechanics (Takao et al. 2009b) (Fig. 2). The relationship between radiation dose D and cell mortality fraction M is described as follows, Radiation dose and mortality fraction are tentatively treated as tensor quantities in formulation process so that direction of shrinkage can be considered. C ijkl is the parameter of radioresistance and varies as the function of radiation dose. Radiation dose D should obey following equations in accordance with the rule in solid mechanics.
Assuming that cancer cells killed by irradiation are removed from the tumor, decrease in the number of cancer cells, i.e., cell mortality directly relates to the changes in tumor volume and shape. The relationship between cell mortality M and the displacement of tumor boundary R is formulated as follows, The boundary condition that prescribes the amount of radiation dose transmitted through the tumor boundary is given by The geometrical boundary condition is given by Equations (2)-(6) express tumor response to irradiation in clinical situation. These equations can be solved numerically by means of discretization method such as finite element method (FEM). Assuming that tumors locally (in a sub-millimeter scale) have an isotropic property, C ijkl parameter in Eq.
(2) for radioresistance is represented by two parameters; reduction resistance E and reduction ratio . The values of E and  should be determined based on the radiobiological model describing the relationship between radiation dose and therapeutic effect. Here we use the linear-quadratic (LQ) model, which is widely accepted and used in the field of radiobiology. According to the LQ model, cell survival fraction S is expressed by following equation, where  and  are parameters of radiosensitivity. In standard radiotherapy, fractionated irradiation, which gives small radiation dose in many times, is performed for reducing radiation damage to normal tissue. Cell survival after n times irradiation (total dose of D, administered in fractionated dose of d) is denoted as follows, This formula is rewritten as the relationship between radiation dose and mortality fraction and in incremental form as, Equation (9) denotes mortality fraction as the function of radiation dose in terms of radiobiology. Equally equation (2) represents same quantity based on the analogy with solid mechanics. Therefore, right-hand sides of both two equations can be equated and following relationship can be derived.
Equation (10) gives the condition that the parameters E and  must satisfy. When the values of E and  are determined, therapeutic displacement Ri can be calculated by means of numerical discretization method such as finite element method.

Heterogeneity of tumor radiosensitivity
Tumor shows internal heterogeneity of radiosensitivity, and it may result in uneven tumor shrinkage. For example, cancer cells under hypoxic environment are likely to be radioresistant, therefore tumor tissue in this region usually does not shrink obviously.
Tumor radiosensitivity also varies depend on many other factors such as cell cycle phase or growth rate of each cell. Because of its complexity, it still seems difficult to measure or estimate the distribution of radiosensitivity in tumor tissue. In this study, intratumoral distribution of radiosensitivity is estimated from the degree of tumor shrinkage in the early stage of treatment. The site where tumor radius decreases greatly is considered to be more radiosensitive, and the site where tumor radius decrease slightly is more radioresistant. Therefore tumor radiosensitivity can be represented as the function of position within the tumor. As radiosensitivity is represented by E parameter in this model, heterogeneity of radiosensitivity can be expressed by the distribution of the values of E.

Simulation procedure
Simulation flow-chart for calculation of tumor geometry is shown in Fig. 3. First, finite element (FE) models of the tumor at the start of the treatment was constructed from the CT images taken for the purpose of treatment planning. Next the values of parameters for tumor radiosensitivity in LQ model,  and , were set as  = 0.1 and  = 0.01 (/=10.0) as suggested by Thames et al. (Thames et al. 1990). Other simulations parameters were consecutively determined The reduction resistance  was tentatively set to 0 assuming that the relationship between tumor response and radiation dose would entirely obey the LQ www.intechopen.com model. After then, the reduction resistance E was calculated from Eq. (10) using , , cumulative dose D, and daily dose d. The E parameter was initially considered to be uniform throughout the tumor. Using these parameters, changes in tumor volume was calculated by means of FE analysis software. Calculated tumor volumes were compared with corresponding tumor volumes measured from CT image sets for treatment follow-up.
Followed by this process, the value of reduction ratio , which represents interpatient variation in radiosensitivity, was adjusted for each case so that the discrepancy between the calculated and the actual tumor volumes obtained from follow-up CT would decrease. If the calculated volume was less than the actual volume, the value of  was incremented and then the tumor volume was recalculated. This iterative process was continued till the root mean square (RMS) error of the calculated and actual tumor volumes reached a minimum.  After that, similar process was performed to determine the distribution of the parameter E (reduction resistance) to represent intratumoral heterogeneity of radiosensitivity. The E parameter was so far considered to be uniform throughout the tumor. In following calculation, the value of E parameter was set for each element, which constitutes the whole 3D tumor model for FEA calculation, as the function of azimuthal angle and polar angle φ. Depending on the distance from the gravity center of the tumor to surface, the value of E at (θ, φ) was determined by following equation.
where, superscript represents the number of iteration, R(θ, φ) M is measured tumor radius at (θ, φ), and R(θ, φ) C is calculated tumor radius at the corresponding point.

Clinical cases 4.1 Patient characteristics
The subjects of this study were three clinical cases (case A, B, and C) of metastatic cervical lymph nodes in three patients with nasopharyngeal cancer treated at the Hokkaido University Hospital, Sapporo, Japan between February 2007 and August 2007. Of the patients, case A and B were undifferentiated carcinoma and case C was squamous cell carcinoma. Initial volumes of lymph nodes were 3.5, 55.1, and 10.4 cm 3 , respectively. Case A and C were treated with IMRT; patient B received conventional radiotherapy. The dose distribution before radiotherapy intended each node in this study to be homogeneously irradiated with a dose of 66 Gy (case A) to 70 Gy (case B, C) in 2.0 Gy fractions delivered five times a week.
Pre-treatment CT images (CT0) were taken for the treatment planning. The slice thickness of the pre-treatment CT images was 2 mm. After the start of treatment, follow-up CT images were taken at weekly intervals (CT1, CT2, CT3, etc). The slice thickness of the follow-up CT images was 5 mm. All patients were immobilized by thermoplastic masks during the CT scanning and treatment. Additionally, in the head and neck IMRT treatments in our hospital, A mouthpiece with three fiducial markers (2 mm diameter gold pellets) was used for the fluoroscopic verification of the patient set-up by means of real-time tumor-tracking radiotherapy (RTRT) system. This study was conducted with written informed consent obtained from all patients and was approved by the institutional ethical committee at Hokkaido University Hospital.

Finite element modeling of tumors
Three-dimensional finite element (FE) models of tumors were constructed based on the CT images taken before and during treatment (Fig. 4). One radiation oncologist determined and contoured the boundary of metastatic cervical lymph nodes on the CT images by means of a treatment planning system (Xio). A group of sequential cross-sectional profiles of the tumor was then loaded into biomedical imaging software and interpolated to 1 mm intervals. After file format conversion with in-house software, the data was imported into the finite element analysis software (ANSYS 11.0), and the 3D FE models of the tumors were constructed.

Results
This study first aimed to evaluate changes in 3D tumor shapes during treatment visually and quantitatively. We constructed surface model of tumors and then analysed the 3D shapes by geometric factors for quantitative evaluation. Two-dimensional surface geometry map was also proposed for visual evaluation of 3D morphological features of tumors. Other main aims of this study was to propose simulation method to predict changes in 3D tumor shape during radiotherapy.

Geometric factors
This study investigated geometric factors to accurately evaluate therapeutic response in radiotherapy. Tumor volume, surface area, radius, and spherical shape factor (SSF) were used to quantitatively evaluate tumor geometry. Fig. 5-7 show changes in geometric factors of tumors through treatment duration for quantitative analysis of 3D tumor morphology. Changes in tumor volume and surface area are shown in Fig. 5 and 6. Tumor volume at the end of the treatment period was 8.7% of its initial volume in case A, 15% in www.intechopen.com case B, and 23% in case C. Changes in tumor surface area showed similar tendency with changes in tumor volume. Tumor surface areas at the end of the treatment were 20%, 28%, and 37% for case A, B, and C, respectively. Changes in SSF are shown in Fig. 7. The values of SSF were about 0.8 through the treatment duration and did not vary widely in all three cases. Fig. 8 shows changes in average, maximum, and minimum tumor radius through the treatment period. At the start of treatment, each average radius was 9.8 mm, 24.6 mm, and 14.4 mm in tumor A, B, and C. The tumor radius ranged from 7.3 mm to 13.3 mm (75 % to 137 % of average radius) in tumor A, 17.4 mm to 34.4 mm (71 % to 140 %) in tumor B, 9.2 mm to 20.3 mm (64 % to 142 %) in tumor C, respectively. At the end of the treatment, average radius decreased to 4.5 mm, 13.5 mm, and 8.9 mm, respectively. The ranges of radius were 3.3 mm to 5.8 mm (73 % to 129 % of the average radius at the end of the treatment) in tumor A, 9.4 mm to 18.2 mm (70 % to 135 %) in tumor B, 6.3 mm to 12.5 mm (71 % to 140 %) in tumor C, respectively. In Fig. 8, all values are represented as percentages of each initial average radius.

Surface geometry map
shows the changes in tumor geometry through treatment duration using surface geometry map. The red end of the color scale is convex region and blue is concave region. Therefore surface geometry map can provide information about 3D geometrical feature of the tumor. The patterns of color distribution are found to be similar in each case. This result means that tumor shrunk uniformly keeping their original morphological features through the treatment duration.

Simulation of changes in tumor geometry
A computational simulation method to predict changes in 3D tumor geometry during radiotherapy was proposed. Simulation parameters were determined based on the changes in tumor geometry in the early stage of the treatment (first 19 days in case A and B, 22 days in case C). Using these parameters, tumor geometries in the latter half of the treatment period were calculated sequentially.
Simulation results are shown in Fig. 10

Discussion
Precise assessment of therapeutic response in radiotherapy has been an important issue in the field of radiation oncology. To understand how tumor geometries change during the treatment would be useful for not only determination of prognosis but for treatment plans as well. However, there has been no research which visualized 3D tumor geometries and evaluated the therapeutic response based on the changes of tumor geometries. In this study, we proposed a method to represent 3D tumor geometry in 2D color map, and evaluated therapeutic response through the treatment period, as well as geometric factors representing therapeutic response quantitatively.
Surface geometry map introduced in this study could indicate 3D morphological features of the tumors in color scale. These figures show that tumors shrank evenly maintaining their original shape. This would be valuable information for determining the optimal radiation field in the latter half of the treatment. The degree of tumor shrinkage, i.e. decrease in tumor radius (shown in Fig. 8), varied approximately plus or minus 20 % depending on tumor region. This variation was considered to represent intratumoral heterogeneity. These findings cannot be obtained from commonly-used geometric factor, i. e., tumor volume or maximum diameter measured on CT images.  A simulation method proposed here enabled to predict changes in 3D tumor geometry considering intratumoral heterogeneity of radiosensitivity therefore would be more valuable for treatment planning or diagnosis of prognosis. The simulation method could predict tumor shapes at the end of the treatment within 2.1 mm discrepancy in radius (absolute average) by considering tumor heterogeneity. As the limitation of present method, these calculations were performed on the assumption that tumors were irradiated evenly. Although dose distribution in IMRT is considered to be uniform to a certain extent, actual dose distribution should be taken into consideration for more precise calculation.
Calculated tumor geometries in the latter half of the treatment strongly depended on changes in tumor geometries in the early stage of the treatment, because values of parameter E were determined only from early tumor geometries. Our present method cannot handle the cases that tumor geometry changes drastically during treatment. We should investigate factors affecting tumor radiosensitivity e.g. angiogenesis, hypoxia, or cell cycle, and appropriately determine the values of E parameter considering these factors.

Conclusion
In this work, tumor geometries were analysed visually and quantitatively from several perspectives. The effectiveness of surface geometry map proposed here was confirmed as the tool for precise assessment of therapeutic response based on the 3D tumor geometry. It was revealed that tumors shrunk uniformly keeping their initial morphological features during radiotherapy. This finding cannot be obtained from traditional evaluation by measuring diameters of the tumors on CT images. This study also proposed a novel simulation method to predict changes in 3D tumor geometries during radiotherapy considering intratumoral heterogeneity of radiosensitivity. The simulation results were found to conform to actual changes in tumor geometry. Although some difficulties still remain to be solved, tumor geometries could be predicted using our approach. It would lead to the idea of "Computer associated radiotherapy (CART)", which is the highly-advanced integration of computational technology and radiotherapy to achieve more precise and safety treatment.

Future directions
The role of computer technology in recent advanced radiotherapy has rapidly increased. Calculation of 3D dose distribution in patient is indispensable for precise treatment. Imaging techniques, especially image registration technique, are expected as tools for adaptive radiotherapy, which modifies original treatment plan to fit the actual patient condition. Through this research, great contribution of our computational analysis and simulation method for changes in 3D tumor geometry has been confirmed. These methods are expected to play a great important role in putting adaptive radiotherapy into practice. It would be a first step for achievement of computer associated radiotherapy (CART), and attainment of more effective and safety treatment.

Acknowledgment
This research is partly granted by the Japan Society for the Promotion of Science (JSPS) through the "Funding Program for World-Leading Innovative R&D on Science and