Abstract
Background/Aim: We established a data-driven method for extracting spatial patterns of dose distribution associated with radiation injuries, based on patients with prostate cancer who underwent iodine-125 (I-125) seed implantation. Patients and Methods: Seventy-five patients underwent I-125 seed implantation for prostate cancer. We modeled the severity of lower urinary tract symptoms (LUTS) to be estimated using a linear model, which is formulated as an inner product between the dose distribution D and voxel-wise radiosensitivity B inside the prostate. For the estimation, tensor regression based on a low-rank decomposition with generalized fused lasso penalty was applied. Results: The spatial distribution of B was visually assessed. Positive parameters appeared dominantly in the region close to the urethra and the prostate base. Conclusion: Our tensor regression-based model can predict intra-organ radiosensitivity in a data-driven manner, providing a compelling parameter distribution associated with the development of LUTS after I-125 seed implantation for prostate cancer.
- Prostate cancer
- brachytherapy
- heterogeneous intra-organ radiosensitivity
- tensor regression
- anatomical standardization
To ensure good quality of life (QOL) for patients treated with radiotherapy, optimizing dose distributions to avoid susceptible sites in organs at risk (OARs) is quite essential. To evaluate radiosensitivity of OARs, conventional techniques rely on dose-volume histograms (DVHs). However, a DVH-based approach basically assumes that the radiosensitivity within a region of interest (ROI) is homogeneous. Moreover, it abstracts the three-dimensional (3D) spatial information of dose distributions into a twodimensional (2D) relationship between the volume and the dose. Therefore, when using DVHs, it is difficult to determine if a particular part of an organ is more sensitive to radiation than other parts. These challenges that come with conventional approaches have limited our knowledge of the heterogeneous sensitivity of OARs to radiation.
The example considered here is the prostate, and the adverse effects after iodine-125 (I-125) seed implantation for prostate cancer. Brachytherapy is widely used as definitive monotherapy or as a boost after external beam radiotherapy (EBRT), because it provides excellent tumor control (1-3). However, in some cases, irritable or obstructive lower urinary tract symptoms (LUTS) occur as an adverse effect of the treatment (4). Generally, the urinary symptoms develop rapidly within 3 months and resolve gradually within approximately 12 months after the procedure (5-7). Several risk factors for LUTS have been reported, including higher biological effective dose, prostate volume, the number of needles implanted, pre-treatment International Prostate Symptom Score (IPSS), and the use of hormonal therapy (8-10). However, there has not been a study on the precise dose-response relationship that considers a heterogeneous intra-prostatic radiosensitivity for the development of LUTS.
Anatomically, the prostate exists below the bladder. The urethra passes through the center of the prostate. The prostate has a base and an apex, where the base, which constitutes the upper surface of the organ, is fused with the bladder neck. Because the internal sphincter muscle located in the bladder neck is important for urinary continence, it is hypothesized that urinary-related toxicities are worse when radiation causes damage to the prostate base.
In order to determine vulnerable sites inside the prostate, several researchers have exploited a DVH-based approach. In the majority of these studies, the prostate was manually segmented into several compartments according to the hypothesis regarding susceptible sites, which were determined beforehand. Then, the dose delivered to each compartment was examined for any adverse effects. Based on these studies, an intra-prostatic region around the prostate base and the bladder neck appears to be particularly sensitive to radiation (11-14). However, there are some inconsistencies among studies (15, 16). These inconsistencies may be because when using the DVH-based approach, advanced insight is needed into the candidate region so that manual segmentation can take place, both of which would be arbitrary due to inter-observer variation. Therefore, we believe that a more data-driven approach exploiting state-of-the-art computational techniques is necessary.
Recently, we proposed an approach to reveal intra-prostatic heterogeneous radiosensitivity after brachytherapy based on principal component analysis (PCA) (17). In this research, each patient’s dose grid inside the prostate was mapped to a common coordinate space using a technique called anatomical standardization, which performs contour-based non-rigid image registration to a template structure. Then, standardized dose grids were flattened into a one-dimensional (1D) dose vector for the subsequent PCA analysis. The study was successful in statistically identifying the responsible eigenvectors for the severity of LUTS without any hypothesis for the vulnerable sites in advance; however, it was still difficult to interpret the spatial distribution of parameters represented as the weighted sum of eigenvectors. Because vectorizing 3D dose information into a 1D array to apply a PCA can cause a loss of spatial information, we extended the previous framework to handle the spatial information of dose distribution as it is.
Here, we demonstrate a tensor regression-based model to establish a data-driven modeling for the heterogeneous intra-organ radiosensitivity of the prostate. Based on previous research using contour-based image registration (17-19), we applied the same methodology for the anatomical standardization of intra-prostatic dose distributions. Then, spatial dose distributions were treated as 3D tensors. For the dose-response relationship, a simple model was prepared, where a summation between voxel-wise radiosensitivity and irradiated dose to each voxel should have a linear relationship for the development of LUTS. In machine learning models using multidimensional array, such as 3D tensor, the dimensionality of input can exceed the number of patients, causing the model to overfit the training data. Therefore, we further applied tensor decomposition in combination with generalized fused lasso penalty. The latter is based on the assumption that adjacent voxels in the prostate should have similar sensitivity to radiation. Resultant parameter distribution representing the radiosensitivity inside the prostate was visually evaluated.
Patients and Methods
Patients and treatment. Because this was an extensive research based on our previous technique, the patient population and treatment details in the present study were the same as in our previous report (17). Briefly, we retrospectively identified 75 patients with prostate cancer from our institution’s records. They received I-125 seed implantation with a prescription dose of 160 Gy from May 2009 to December 2013 at our institution. The details on our treatment protocol and brachytherapy techniques are described elsewhere (20). Patient characteristics are summarized in Table I.
Patient characteristics.
Scoring of LUTS. For the quantitative evaluation of LUTS, a composite score of the results from the IPSS questionnaire was used. The IPSS is based on a patient’s response to seven questions concerning urinary symptoms and one question on QOL. The total score of IPSS can range from 0 to 35. In our institutional protocol, a patient’s IPSS score is obtained before brachytherapy. Then, patients respond to the IPSS questionnaire at each follow-up visit every 3 months for the first year after treatment. In this study, IPSS scores were collected retrospectively from the database. For the calculation, the pre-treatment score was a baseline and the maximum increase in IPSS score during the first 12 months after treatment was obtained for each patient. Then, these values representing an individual increase in IPSS scores were normalized with a mean of 0 and a standard deviation of 1 among patients (Z-score normalization). Finally, the normalized values were concatenated to be a 1D array, denoted as y.
Pre-processing for the dose distribution. Data on intra-prostatic dose distribution delivered in brachytherapy were obtained based on computed tomography scanning at 1 month after treatment. At the scanning procedure, a foley catheter was inserted to visualize the urethra. Because the intensity of radiation is inversely proportional to the square of the distance from the source, the dose values were converted to the natural logarithm for smoothing the value distribution. Based on the assumption that small doses less than 1.0 Gy have little effect on organ function, we rounded values less than 1.0 Gy to 0.0 Gy for simplicity. Then, as in the previous report (17), anatomical standardization to map each patient’s intra-prostatic dose distribution to a common coordinate with a template prostate was applied. As a result, patient radiation doses were standardized in the coordinate with a size of 46×35×36. By excluding the voxels outside the template prostate, we considered a total of 31,018 effective covariates for the prediction model of LUTS severity. See Figure 1 for examples of standardized dose distributions.
Example results of anatomical standardization. Intra-prostatic dose distributions, which are different between patients, were mapped to a common coordinate space with a template prostate through a contour-based non-rigid registration. After anatomical standardization, quantitative comparison of dose distributions can be enabled. Note that the global topologies of dose distributions were preserved before and after the standardization.
Tensor regression model for intra-prostatic radiosensitivity. Generally, normal tissue complication of parallel organs can be modeled as a summation of local radiation damage (21). Here, we considered the prostate as a parallel functioning organ regarding the development of LUTS, where a summation between voxel-wise radiosensitivity and irradiated dose to each voxel should have a linear relationship for the development of LUTS. Then, we tried to model this linear relationship based on 3D tensors for the preference of spatial information.
Tensors provide a natural representation for multidimensional data. A first-order tensor is a 1D vector, a second-order tensor is a 2D matrix, and tensors of order three or higher are called higher-order tensors. We define a parameter distribution representing intra-prostatic radiosensitivity as a tensor B∈ℝp1×…×pD and an intra-prostatic dose distribution as X∈ℝp1×…×pD. Here, these tensors are considered as 3D (D=3). Then, based on the assumption of the prostate as a parallel organ, the severity of LUTS, y, is formulated as follows: ŷ =<B,X>, where an inner product between two tensors is calculated as <B,X>=∑i1,…,iDβi1…iDxi1,…iD.
When maximum likelihood estimation is applied to the above equation, it will be severely compromised by the huge number of regression parameters of B because the number of patients in this study was relatively limited. To mitigate the ultrahigh dimensionality of the tensor in the regression problem, we exploited a computational approach proposed by Zhou et al. (22). Briefly, the approach starts with rank-R tensor decomposition, which approximates original signals as a linear combination of eigenvalues and the outer product of their corresponding eigentensors. Here, R denotes the number of eigentensors needed for the approximation. Based on this technique, a tensor B∈ℝp1×…×pD admits a rank-R decomposition as follows:
where 
indicates 1D vectors along with each direction. An outer product of 1D vectors 
constitutes an r-th rank-one tensor, and the sum of R rank-one tensors can approximate the original tensor B. Subsequently, when a tensor is decomposable as a rank-R tensor, based on a mode-d matricization and a vec operator (23), the tensor can be rewritten as:
denotes Khatri-Rao product between tensors (24), and 1R is a vector of R ones. Based on these equations, our formulation of the predictive model can be written as follows:
This simplifies the estimation of all parameters of B as a regression problem of 1D vectors. Zhou et al. (22) solved this regression problem by using an algorithm called block relaxation, which alternatively updates Bd (d=1,…,D) while keeping other components fixed. This enables us to break the simultaneous estimation of all parameters into a sequence of low-dimensional parameter optimization. As a result, the update rule can reduce the dimensionality from the order of ∏dpd to the order of ∑dpd for regression.
In addition to the tensor decomposition described above, we further considered regularization for the parameter estimation. Generally speaking, regularization is necessary in the computation of tensor regression for stabilizing the estimates and avoiding overfitting. The selection of the kind of constrains is usually based on insights into the subject the model addresses. To make the matter simple, we consider that intrinsic radiosensitivity in the prostate should be contiguous, where units adjacent to each other may have similar sensitivity to the radiation injury. Thus, a penalty on the pairwise difference to make successive variables similar was thought to be reasonable. This type of regularization is known as generalized fused lasso, which can be incorporated into the prediction model as follows:
where β∈ℝd and λ1,λ2≥0. λ1 indicated the sparcity of critical regions, and λ2 regulated the similarity between adjacent voxels. L1 norms penalize both the variables and their pairwise differences. The calculation was performed using SPArse Modeling Software (SPAM) (25, 26).
Finally, we applied these computational techniques to minimize a mean squared error (MSE) between the recorded severity of LUTS y and the estimate ŷ. When the best model is acquired in terms of prediction accuracy, the parameter distribution of B is visually assessed to reveal the vulnerable sites in the prostate.
Results
Hyperparameter optimization. Our model has several hyperparameters, such as the number of rank in the tensor decomposition R and the balancing terms in the generalized fused lasso penalty, λ1 and λ2. These hyperparameters were optimized via k-fold cross-validation (k=15), splitting the original dataset (n=75) into a validation (n=15) and a test dataset (n=5). Hyperparameter values were ranged as follows: R={1,2,…,5}, λ1={0,1.0×10–6),1.0×10–5),…, 1.0×10–1)}, λ2={0,1.0×10–6),1.0×10–5),…,1.0×10–1)}. The best combination of hyperparameters was obtained based on grid search in terms of the minimization of MSE between the estimation of the model ŷ and the actual values y in the test dataset. Table II shows the three best sets of hyperparameters obtained. Figure 2 demonstrates the parameter distribution of B from each model configuration.
Three best parameter combinations minimizing MSE.
The parameter distributions representing intra-prostatic radiosensitivity of the three models. The magnitudes of parameters were highlighted and overlaid on the template prostates. Note that the parameters of greater value are clustered around the urethra and the prostate base.
Comparison of model accuracy with a baseline. The prediction accuracy of the obtained models was compared with a constant model, which was defined as y=0. This means that the estimates of the constant model always correspond to the mean value of the severity of LUTS in the dataset, which is normalized to be 0. Using this constant model as a baseline, in model 1, the maximum improvement in terms of MSE is obtained by 17.8%.
Evaluation of the effect of generalized fused lasso penalty. Figure 3 presents an example effect from the generalized fused lasso penalty based on the visualized parameter distribution of B in model 3 (see Table II for the configuration of the model), where the regionality of parameter distribution showed up more clearly in the regularized model. We also confirmed that MSE of the model was improved from 1.06 to 0.74 by imposing the penalty.
An example/result of the regularization effect by imposing generalized fused lasso. The region associated with the development of lower urinary tract symptoms (LUTS) appeared more clearly in the regularized model (rank=3, λ1=1.0×10−5, λ2=1.0×10−1) than in the unregularized model (rank=3). We also confirmed that the model accuracy was improved in terms of mean squared error (MSE), from 1.06 to 0.74.
Visual assessment of the parameter distribution. Although the three models differed in terms of signal coarseness, there seemed to be a shared pattern, where large parameters gathered dominantly in the region close to the urethra and the prostate base. To confirm this preference, we created an averaged model from the three models. Figure 4 shows overlapped regions with large parameter values around the urethra and the prostate base.
The result of the averaged model parameter distribution. Parameter distributions of the three models were averaged. Then, the mean distribution was demonstrated on the template prostate. There was the same tendency of distribution when large parameters were preferentially located around the urethra and the prostate base.
Discussion
Considering the differences in radiosensitivity within tumors or OARs should be the next step for refining state-of-the-art radiotherapy, in order to optimize dose distribution based on the balance between tumor control and complications. Current DVH-based models presume that all parts of an OAR have the same intrinsic properties. However, there has been accumulating evidence showing that healthy organs can be heterogeneous in function as well as in sensitivities for toxicities (27).
Some attempts have been made to introduce spatial consideration into the model of normal tissue complication in radiotherapy. For example, by exploiting a technique called dose-surface maps, Buettner et al. (28) revealed that late complications in the rectum after prostate radiotherapy are related to the dose distribution shape. They also proposed a quantitative prediction model based on a 3D dose parameterization in the rectal wall, which outperformed the standard DVH-based model (29). More recently, Liang et al. (30) leveraged image registration techniques to map each patient’s dose grid into a standardized coordinate space and identified specific dose patterns associated with acute hematologic toxicity based on PCA regression. In line with this research, we applied anatomical standardization using contour-based image registration in combination with PCA regression to identify the vulnerable sites in the prostate after brachytherapy (17). Because the previous study was limited in terms of the loss of spatial information when applying PCA to the dose array, we tried to mitigate this point by extending the framework into a tensor regression problem.
Although the tensor-based approach can introduce spatial consideration into the dosimetric analysis, it has only been used in a limited number of prior studies. For example, Ospina et al. (31) provided tensor-based population value decomposition to extract characteristics of spatial dose patterns associated with particular clinical outcomes. More recently, Fargeas et al. (32) classified spatial dose patterns based on the distances along with subspaces spanned by basis vectors derived from tensor decomposition, demonstrating that the tensor-based model can outperform the classic DVH-based model in the prediction of clinical outcome.
Based on the assumption that radiosensitivity inside OARs should be heterogeneous, we proposed a novel framework leveraging tensor regression in a data-driven manner to parameterize the radiosensitivity for an adverse effect. Notably, several state-of-the-art techniques, including block relaxation algorithm based on rank-R tensor decomposition and generalized fused lasso, effectively eased the tensor regression problem. We consider that these techniques were essential because the number of patients was too small to estimate the huge number of parameters contained in tensors in a naïve manner. This substantial imbalance between the number of samples and the number of voxels in medical imaging, including dose distribution of radiotherapy, may be quite common because it is generally difficult to obtain a large amount of patient data from a hospital. Therefore, we consider that our proposal is a versatile approach for extracting spatial patterns of covariates that are related to some clinical outcomes.
The results showed that the region around the urethra and the prostate base has large values of parameters representing radiosensitivity for the development of LUTS. These results are convincing because the internal sphincter muscle located near the prostate base may play an essential role in urinary continence. Moreover, the implication that high doses to the urethra may be associated with urinary symptoms is also compelling. Note that we obtained these results in an intuitive manner without any hypothesis for candidate intra-prostatic sites beforehand. In addition to that, our model also provides a quantitative relationship between the spatial dose distribution and the severity of an adverse effect as a normal tissue complication model. Thus, by revealing sensitive areas within OARs, our model can enable treatment planning of radiotherapy of I-125 seed implantation for prostate cancer to better meet the patient’s QOL. Additionally, our simple linear model is extensible. Therefore, in the future, the model should incorporate other clinical factors that also affect the development of LUTS to establish precision brachytherapy for better QOL for patients.
Conclusion
Our tensor regression-based model can predict intra-organ radiosensitivity in a data-driven manner, providing a convincing parameter distribution associated with the development of LUTS after I-125 seed implantation for prostate cancer. Because our technique is versatile for other OARs, it will widen the opportunity for refining state-of-the-art radiotherapy by optimizing dose distribution in a more precise manner.
Acknowledgements
This work has been financially supported by MEXT KAKENHI Grant Number 17K16497, and JST CREST Grand Number JPMJCR1689.
Footnotes
Authors’ Contributions
All Authors made substantial contributions to conception and design, acquisition of data, and analysis and interpretation of data: K.K., N.M., K.T. and K.I. collected cases. K.K. completed all data. K.K. designed the study and analyzed data. K.K. wrote the article. H.I, R.H. and J.I. critically revised the article.
This article is freely accessible online.
Conflicts of Interest
The Authors have no conflicts of interest to declare regarding this study. The Authors alone are responsible for the content and writing of the paper.
- Received October 9, 2020.
- Revision received October 19, 2020.
- Accepted October 20, 2020.
- Copyright© 2021, International Institute of Anticancer Research (Dr. George J. Delinasios), All rights reserved
References
In this issue
Jump to section
Related Articles
Cited By...
- No citing articles found.