A subregional prediction model for radiation-induced hypothyroidism

Abstract

Background

Considering the potential association between radiation-induced hypothyroidism (RHT) and the thyroid subregions as well as the received radiation dose in each subregion, this study aims to develop a subregional prediction model for RHT.

Methods

CT images and dose images of 128 patients with nasopharyngeal carcinoma were collected retrospectively. The thyroid subregion was obtained by clustering thyroid voxels and voxel entropy. After extracting 1781 radiomics features and 1767 dosiomics features, a subregional RHT prediction model was established, and its performance was compared with that of the whole thyroid model. The phenotype and dosimetry parameters of each subregion were analyzed by AUC, T test and Delong test.

Results

Three subregions (S1, S2, S3) were identified. The subregional prediction model was constructed based on 34 radiomics and dosiomics features. According to the Delong test, the prediction performance of the subregional model was significantly superior than that of the whole thyroid model (0.813 VS 0.624, p = 0.038). Subregional analysis suggests that S1 and S3 regions may have higher radiosensitivity than S2 regions.

Conclusions

In this study, a subregional model for predicting RHT was established and the radiosensitivity of the relevant subregions was evaluated. The subregion-based RHT prediction model may help to improve radiotherapy plan design for better thyroid function protection.

Background

Radiotherapy is one of the main treatment methods for nasopharyngeal carcinoma (NPC) [1]. Currently, the five-year survival rate of NPC patients who received radiotherapy had reached 80%, resulting from the continual advancement of radiotherapy technology [2, 3]. Radiation-induced hypothyroidism (RHT) is a prevalent complication observed in patients undergoing radiotherapy for the treatment of head and neck tumors, with an incidence of approximately 40% and a peak occurrence ranging from several months to years after the completion of radiotherapy [4,5,6]. RHT exerts a sustained adverse impact on patients’ quality of life and patients with RHT have to receive lifelong thyroid hormone replacement therapy [7, 8]. Therefore, establishing a predictive model for RHT is of great significance for personalized planning design and radiotherapy application of NPC patients.

RHT is associated with multiple factors, among which radiation dose has always been considered one of the primary contributors. However, there is currently no consensus on the optimal dose limits for the thyroid tissue. Several dosimetric parameters, including Dmean (average dose), V50 (the proportion of thyroid volume exposed to a dose higher than 50 Gy), or VS45 (the absolute thyroid volume spared from 45 Gy or less), have been consistently demonstrated in numerous studies as reliable predictors for RHT [9,10,11]. Nevertheless, these dosimetric parameters just evaluate thyroid tissue as a comprehensive entity, overlooking the spatial heterogeneity inside thyroid gland. The latest research suggests the potential existence of regional disparities in the radiosensitivity of normal tissues. It is reported that late radiation damage was related with the location of dose delivery in rat parotid glands, which could be potentially applicable in human parotid as well [12]. Furthermore, the concept of thyroid functional subregions was raised up in a study using ¹⁸F-FDG-PET to effectively predict the occurrence of RHT, suggesting a potential relationship between radiotherapy complications and the dose received in each internal functional subregions [13].

Radiomics is a research methodology that was capable of extracting quantitative features from CT, MRI, PET, or ultrasound images through a high-throughput manner to help realize tissue or organ characterization at pathophysiological level. Currently, RHT prediction based on radiomics method has drawn great interest. By extracting radiomics features from the whole thyroid region, a RHT prediction model was developed for head and neck tumor patients who received radiotherapy [14]. By combining radiomics features derived from contrast-enhanced CT images and dosimetric parameters, a predictive model for RHT with superior predictive performance was established, exhibiting the advantage of the integrated model [15].

The dosiomics method, inspired by the radiomics approach, enables the extraction of features from three-dimensional dose distribution images, providing a more comprehensive and multidimensional information compared to conventional dosimetric parameters. Presently, we have extracted dosiomics features from dose distribution images of 145 NPC patients and established a predictive model for RHT using various machine learning algorithms [16]. Our findings indicate that the prediction models based on dosiomics features demonstrate superior predictive performance than the models based on dosimetric parameters. Furthermore, dosiomics can be integrated with radiomics for bi-omics analysis. Existing study has reported a dual omics model for RHT prediction by integrating radiomics and dosiomics [17]. According to their research findings, the dual-omics model not only outperforms the normal tissue complication probability (NTCP) model, but also provides guidance in designing radiotherapy plans for patients with NPC. However, all these studies focused on feature extraction of the whole thyroid gland without considering the impact of subregional phenotypic distribution within the thyroid gland on RHT occurrence.

Currently, research on predicting radiotherapy complications based on internal subregions is constantly emerging. In patients with oropharyngeal squamous cell carcinoma, a predictive model for mandibular radiation necrosis was developed by delineating the mandible into sixteen subregions based on dental anatomy [18]. The findings of aforementioned study indicate that the manifestation of mandibular necrosis varies among individuals and is intricately associated with the subregional dose distribution. Another study also investigated the potential of subregional radiomics models based on parotid glands to accurately predict late xerostomia following radiotherapy in patients with head and neck cancer, demonstrating exceptional precision and earlier identification using subregional radiomics models [19]. Furthermore, functional subregions identified based on a threshold ranging from 20 to 80% of the maximum perfusion count on single-photon emission computed tomography (SPECT) images also exhibited superior performance in the prediction of radiation-induced pneumonia. The does volume parameters derived from these functional lung subregions outperformed conventional dose-volume histogram (DVH) metrics [20]. Similarly, the association between specific subregions and acute dysphagia resulting from radiotherapy was also proved in another study by employing voxel-wise analysis to identify swallowing-related tissues [21].

If variations in radiosensitivity exist within the thyroid and can be identified, selectively protecting more susceptible subregions of the thyroid may aid in mitigating the risk of developing RHT. Nevertheless, to the best of our knowledge, the prediction of risk for RHT based on subregions remains unexplored. Therefore, the goal of this research is to explore and compare the varying levels of accuracy among predictive models for RHT that utilize subregions and whole thyroid gland volume, while analyzing the subregion response relationship.

Methods

Data collection

This retrospective study collected clinical data, CT images, and dose images of 128 patients with NPC who underwent radiation therapy between January 2012 to January 2015. The treatment plan for each patient was designed based on the TomoTherapy Planning System (Accuray Inc., Madison, WI) or Pinnacle³ (v9.0) treatment planning system, with a prescribed dose of 70–74 Gy administered into 33 sessions. The inclusion criteria for this study were as follows: 1) No history of previous head and neck radiation therapy, thyroid surgery, or any thyroid-, hypothalamus-, or pituitary-related diseases; 2) Serum thyroid function tests (TFTs) were conducted prior to treatment; 3) The comprehensive follow-up results of serum TFTs data for determining the functional status of the thyroid gland.

The clinical endpoint of this study is the occurrence of either clinical or subclinical hypothyroidism. Clinical hypothyroidism is defined as an elevation in serum thyroid stimulating hormone (TSH) levels (> 5.5μU/mL) during follow-up, accompanied by a decrease in serum free thyroxine (FT4) levels (< 12.0 pmol/L). Subclinical hypothyroidism is characterized by an elevation in serum TSH levels with normal FT4/T4 levels observed at least twice during follow-up monitoring [22]. In this study, enrolled patients underwent regular assessments of free triiodothyronine (FT3), FT4, and TSH tests every 3 or 6 months within the first two years following radiotherapy, followed by annual evaluations starting from the third year onwards. The maximum follow-up duration was five years since the completion of initial radiotherapy.

A total of 128 enrolled patients were assigned numbers and subjected to stratified random sampling to create a training set (n = 76) and a testing set (n = 52) at a ratio of 3: 2. A comparable proportions of RHT were realized between training set and testing set by the stratified random sampling. The training set was utilized for feature selection and modeling, while the test set was used for independent validation and performance comparison between different models. Additionally, to further validate the subregional model's performance, this study also developed a model based on the whole thyroid gland. The primary methodology employed in this study is illustrated in Fig. 1.

Fig. 1

Workflow in this study

Image acquisition

The CT images utilized in this study are radiotherapy localization images acquired by SIEMENS Definition AS and Philips Brilliance Big Bore scanners. The tube voltage is set at 120 kVp, the pixel resolution ranges from 0.8 to 2.5 mm, and the slice thickness is 3 mm. The patient's dose image is exported from the TPS system, with different radiotherapy techniques resulting in dose space rates (grid size) of 2.73 × 2.73 × 3 mm³ for TOMO and 4 × 4 × 4 mm³ for IMRT.

Subregion segementation and feature extraction

The region of interest (ROI) in this study pertains to the thyroid region derived from the radiotherapy plan. Subregions are obtained through clustering of the ROI as follows: Firstly, taking each voxel within the ROI as the center, calculate the entropy of 9 × 9x9 patches individually and consider it as the local entropy of that voxel. Next, compile a two-dimensional matrix by combining the intensity and local entropy values of voxels within the ROIs from all patients in the training set, followed by applying horizontal clustering using K-means algorithm. Horizontal clustering refers to the clustering patterns performed among patients, with selection of its K value based on Calinski-Harabasz (CH) score trends.

After obtaining subregions, Pyradiomics was utilized to extract features from both the whole thyroid and its subregions separately. A total of 1781 radiomics features and 1767 dosiomics features were extracted from each patient region to enhance model accuracy. Detailed information regarding the features and extraction parameters can be found in Supplementary Table S1 and Supplementary Table S2.

Feature selection and modeling

The feature selection and modeling process were conducted separately based on either the whole thyroid gland or their corresponding subregions. A two-step method was employed for feature selection on the training set. Firstly, t-tests were utilized to analyze the radiomics and dosiomics features between the RHT and the Non- Radiation-induced hypothyroidism (NRHT) patients in the training set. For clinical features, chi-square tests, Mann Whitney U-tests, and Student's t-tests were selected based on specific circumstances for analysis, with significant differences (p < 0.05) being retained for further investigation. Subsequently, the least absolute shrinkage and selection operator (LASSO) algorithm was applied additional feature selection. LASSO is a linear regression technique that incorporates regularization to control the number of features by adjusting regularization parameters accordingly. In this step, a threefold cross-validation strategy combined with grid search was employed using the reserved feature matrix from previous steps to determine optimal parameters. Once the most suitable subset of features was determined through LASSO analysis, logistic regression was used for modeling.

Subregional analysis

The voxel intensity and entropy between subregions were analyzed using a t-test. To investigate the phenotypes and radiosensitivity of the subregions obtained in this study, the differences in V50, VS45, VS60, and Dmean between these subregions were analyzed [9,10,11]. The calculation method for dosimetric parameters in the subregions is as follows:

$$\begin{array}{c}{V}_{a}\left(Sn\right)=\frac{\sum_{i\in {Sn}_{a}}{v}_{i}}{\sum_{i\in Sn}{v}_{i}}\\ {VS}_{b}\left(Sn\right)={\sum }_{i\in {Sn}_{b}}{v}_{i}\\ Dmean\left(Sn\right)=\frac{\sum_{i\in Sn}{v}_{i}{d}_{i}}{\sum_{i\in Sn}{v}_{i}}\end{array}$$

Among them, $Sn$ represents subregion $n$, ${Sn}_{a}$ demotes the areas within subregion $n$ with doses greater than or equal to $a$, and ${Sn}_{b}$ refers to the areas with subregion $n$ with doses less than $b$. The variables ${v}_{i}$ represents the voxel size, while ${d}_{i}$ represents the dose received by each voxel.

After completing the calculation of dosimetric parameters, t-tests were employed to analyze the dosimetric parameters across different subregions in both RHT and NRHT cohorts. Subsequently, the diagnostic performance of dosimetric parameters for RHT was evaluated using the area under the receiver operating characteristic (ROC) curve (AUC), while any potential differences in performance among subregions were assessed using the Delong test.

Statistical analysis

All data in this study were analyzed using R software (version 4.1.2) and Python (version 3.7), with a significance level set at 0.05. Depending on the situation, the differences in clinical factors and dosimetric parameters between RHT and NRHT were evaluated using the Student's t-test and Mann–Whitney U test, respectively. The performance evaluation of the model is based on the AUC curve of the subjects, and the Delong test is used to compare the performance differences between different models.

Results

Patient characteristics

The main clinical factors and dosimetric parameters of the patients were presented in Table 1. There were no significant differences observed between the training and testing sets with respect to gender, age, T-stage, N-stage, radiotherapy technique, treatment mode, and serum thyroid function outcome (RHT or NRHT).

Table 1 Characteristics of patients in the training and test sets

Subregional clustering and feature selection

The process of determining the optimal value for K in K-means clustering is illustrated in Fig. 2. When K = 3, the CH score reaches its maximum value. Consequently, based on this analysis, the number of thyroid subregions in this study was set to three, denoted as S1 region, S2 region, and S3 region respectively.

Fig. 2

Subregion segmentation. a is a CT image, (b) is a thyroid gland image, (c) is a DOSE image and (d) is a clustering subregional image. e is the trend chart of the Calinski–Harabasz score in the training set under different K value

For the whole thyroid model, a total of 24 radiomics and dosiomics features were determined using a two-step methodology. For the subregional model, a comprehensive set of 34 features were obtained from three subregions. Specifically, 14 features were originated from the S1 region, while 6 and 14 features were derived from the S2 and S3 regions, respectively. The detailed outcomes of feature selection for both the whole thyroid model and subregional model can be found in Supplementary Table S3 and Supplementary Table S4.

The clinical factors, as presented in Table 1, did not exhibit any statistically significant differences between RHT and NRHT groups with regards to gender, age, T stage, N stage, radiation therapy technique, and treatment mode. The comparative analysis of clinical factors can be found in Supplementary Table S5.

Model comparison

The AUC value of the whole thyroid model on the test set is 0.624, as depicted in Fig. 3, whereas the subregional model achieves an AUC value of 0.813. Statistical analysis using the Delong test reveals a significant difference in AUC between the whole thyroid and subregional models (p = 0.038). Detailed performance metrics for these models are summarized in Table 2.

Fig. 3

ROC curves of different models on the test set. a represents the whole thyroid gland model, and (b) represents the subregional model

Table 2 The prediction performance for whole-thyroid and subregional models on test set

Subregional analysis

The average voxel intensity and voxel entropy of the S1, S2, and S3 regions are presented in Fig. 4. Table 3 displays the average dosimetric parameters for these regions. As depicted in Fig. 5, there were no statistically significant differences between RHT and NRHT groups regarding Dmean and V50 in the S1, S2, and S3 regions. However, a statistically significant difference was observed between RHT and NRHT groups for VS45 in the S1 and S3 regions but not in the S2 region. Additionally, a statistically significant difference was found between RHT and NRHT for VS60 specifically in the S3 region while no significant difference was noted between the S1 and S2 regions.

Fig. 4

Subregional phenotype analysis of thyroid gland. The ROI is divided into three subregions, as shown in (a). The red region represents the S1 region, the blue region represents the S2 region, and the yellow region represents the S3 region; (b) shows the comparison of voxel intensity in different subregions of the thyroid gland; (c) shows the comparison of voxel entropy in different subregions of the thyroid gland

Table 3 Comparison of dosimetric parameters between RHT and NRHT in different Subregions

Fig. 5

Comparison of dosimetric parameters in different subregions. a, (b), (c), and (d) represent the comparisons between Dmean, V50, VS45, and VS60 in RHT group and NRHT group for different subregions. * p < 0.05, * * p < 0.01, * * * p < 0.001, n. s represents no significant difference

The AUC values of dosimetric parameters for S1, S2, and S3 regions were presented in Table 4. Based on Delong test, the AUC values within S3 region showed a statistically significant difference with regard to VS45 and VS 60, compared to either S1 or S2 regions.

Table 4 AUC values of dosimetric parameters in different subregions

Discussion

To the best of our knowledge, this study represents the first subregions-based investigation on predicting RHT. In this study, we employed clustering methods to identify distinct subregions within the thyroid gland and compared the disparities of predictive performances between whole-thyroid-gland based model and subregions-based model. The findings from this study demonstrate that the prediction model relying on surpasses its counterpart model based on the whole thyroid gland (0.813 vs 0.624, p = 0.038). Specifically, although clinical factors such as gender and age were considered as predictive values for RHT in relevant studies [23], these clinical factors did not show statistically significant differences between the RHT and the NRHT groups in our study. Therefore, this study did not incorporate clinical factors for predictive model construction.

Previous studies have demonstrated that adult stem cells/progenitor cells exhibit greater radiation sensitivity when compared to normal tissue cells. Irradiation on the regions where these adult stem cells/progenitor cells are located is more likely to induce radiation damage [24]. The presence of distinct cell compositions may be correlated with varying tissue densities [25]. These variations in regional tissue density can be indicated by the voxel intensity, while the regional heterogeneity can be reflected by voxel entropy. In this study, subregional analysis revealed that the S1 region displayed moderate average voxel intensity and voxel entropy. The S2 region exhibited the lowest average voxel intensity and highest average voxel entropy. Conversely, the S3 region had the highest average voxel intensity but lowest average voxel entropy. In terms of VS45, significant differences were observed between RHT and NRHT groups in either S1 or S3 region. However, no statistically significant difference was found in the S2 region. For VS60, only the S3 region showed statistically significant difference between RHT and NRHT groups. According to Delong's test, the AUC values for VS45 and VS60 showed a statistically significant difference between the S3 region compared to either S1 or S2 regions. These findings suggest that cell density, radiosensitivity, and heterogeneity were moderate in the S1 region; heterogeneity was highest in the S2 region with lower cell density and radiosensitivity; whereas heterogeneity was lowest with higher cell density and radiosensitivity observed in the S3 region. Therefore, limiting dosage in both the S1 and S3 regions or transferring high-dose areas to the S2 region may help reduce the risk of developing RHT.

In this study, functional images or voxel analysis were not utilized; instead, voxel clustering methods were employed to acquire subregions. Although several studies have demonstrated the value of functional image-based subregions in predicting radiotherapy complications prediction [20, 26, 27], there is a lack of research specifically focusing on the association between specific functional images and RHT. In an existing study that aimed to predict RHT based on ¹⁸F-FDG PET/CT images, the subregion segmentation process was subjectively determined through setting a specific threshold, which would potentially affect the objectivity of extracted radiomics/dosiomics features from these subregions [13]. For voxel analysis, the segmentation results of subregions heavily rely on the registration effect, which poses significant challenges for organs with high volume variability such as thyroid gland. In contrast to these two methods, our study employs clustering algorithms to derive subregions. This data-driven approach not only ensures relative objectivity but also demonstrates excellent predictive performance, as mentioned in numerous radiotherapy prognosis studies.

Radiotherapy is a primary treatment modality for NPC.The clinical management of NPC requires a delicate balance between minimizing radiation-induced complications and maximizing local tumor control when designing radiotherapy plans. The continuous advancements in high-precision radiotherapy technologies, such as IMRT, are gradually making personalized dose optimization for specific subregions within normal tissues feasible to mitigate the occurrence of radiotherapy complications. For instance, Van Luijk et al. discovered that the radiation dose in a specific subregion of the salivary gland containing adult stem cells/progenitor cells is associated with salivary gland function one year after radiotherapy [28]. Reducing the radiation dose in this particular region can help decrease the incidence of Sjogren's syndrome. Similarly, adult stem cells/progenitor cells have also been confirmed to exist in the thyroid gland and may play a role during in vivo regeneration following thyroid injury [29, 30]. Furthermore, previous studies have indicated that decreased thyroid gland function is more likely to be linked to its immune volume when exposed to certain doses of radiation [6, 31, 32], suggesting that identifying the subregions where thyroid adult stem cells/progenitor cells are located and optimizing dosage can aid in protecting thyroid function. However, due to the unclear spatial distribution of adult stem cells/progenitor cells within the thyroid gland, there still exist challenges in designing functional protective radiotherapy plans for specific subregions of the thyroid. In this study, we have established a predictive model for subregional radiation-induced hypothyroidism (RHT) and assessed the radiosensitivity of relevant subregions, which may offer valuable insights for optimizing radiotherapy plans aimed at preserving thyroid function.

This study has several limitations. Firstly, this study is a retrospective single-center study. Consequently, the generalization ability of the model cannot be determined with certainty. Therefore, future research should focus on collecting a substantial amount of data from multiple centers to further evaluate and validate the application value and stability of the model. Secondly, due to the absence of sequencing information and imaging data for cell clusters in thyroid subregions, direct analysis of biological phenotypes associated with these subregions was not possible in this study. To address this limitation, future investigations could consider employing single-cell sequencing or spatial transcriptome analysis specifically targeting the thyroid gland to provide insights into these subregions. Lastly, it is worth mentioning that this study primarily focused on analyzing factors related to the thyroid gland without considering other variables such as pituitary radiation dose which may impact model accuracy. In subsequent studies, efforts will be made to incorporate an analysis of pituitary dose distribution in order to enhance overall model performance.

Conclusion

In this study, we utilized data-driven methods to identify subregions within the thyroid gland and developed a predictive model for RHT based on these subregions. Our research findings suggest that the subregional model exhibits superior predictive performance for RHT and reveals significant phenotypic variations among different subregions, which may offer valuable technical support in designing radiation therapy plans that aimed at preserving thyroid function in patients with NPC.