Risk assessment of thyroid cancer recurrence based on quantitative multiphoton imaging

Thyroid cancer is a common endocrine malignancy with rising incidence, and although its overall prognosis is relatively favorable, postoperative recurrence remains a clinical challenge due to limited insight into tumor microstructure and stromal remodeling from current risk assessment methods. To address this issue, this study investigated the association between collagen fiber morphological and textural characteristics and tumor recurrence using multiphoton microscopy (MPM), a high-resolution, label-free imaging technique capable of visualizing tissue microstructure and stromal architecture in detail, providing detailed quantitative information on both cellular and extracellular components. To identify key features linked to recurrence, a multi-stage pipeline-univariate testing, correlation-based feature reduction, and least absolute shrinkage and selection operator regression were employed. The distributions of the selected features were visualized using box plots and kernel density curves. Finally, based on the selected features, a logistic regression model was developed and its ability to discriminate between high- and low-risk patients was evaluated using stratified five-fold cross-validation, with performance quantified by the receiver operating characteristic curve and the area under the curve (AUC). Qualitatively MPM imaging revealed pronounced collagen linearization and aligned fiber remodeling in high-risk tissues, contrasting with the wavy, interwoven architecture of low-risk tissues. Through the feature selection process, ten key features were identified, all showing significant distributional differences between risk groups (P < 0.0001). The resulting logistic regression model achieved strong discrimination, with an AUC of 0.959 in cross-validation. These findings demonstrate that the developed model effectively stratifies recurrence risk in thyroid cancer and may serve as an objective, quantitative tool to support individualized postoperative monitoring and prognostic assessment.

Thyroid cancer, one of the most common malignancies of the endocrine system, has shown a significant global increase in incidence in recent years, which may be attributable in part to environmental and lifestyle factors [1–3]. Most thyroid cancers originate from follicular epithelial cells, with papillary thyroid carcinoma (PTC) being the predominant subtype [4, 5]. Although PTC generally exhibits low invasiveness and favorable outcomes, postoperative recurrence and lymph node metastasis remain important clinical concerns, necessitating vigilant long-term surveillance [6, 7]. Current risk stratification for PTC recurrence primarily relies on clinicopathological staging systems, such as AJCC/ATA. However, these systems do not fully capture inter-individual variability and have limited predictive performance, underscoring the need for novel microenvironment-based biomarkers [8, 9].

The tumor microenvironment (TME) is a complex ecosystem composed of tumor cells, interstitial fluids, immune cells, blood vessels, and other stromal components [10–12]. Recent studies have suggested that tumor cells not only drive tumorigenesis and progression through genetic mutations but also promote malignant transformation and invasive behavior by continuously remodeling the extracellular matrix (ECM), thereby creating a feedback loop that further supports tumor aggressiveness [13–15]. As a major structural component of the ECM, collagen fibers play a critical role in the TME, not only serving as a passive scaffold but also actively regulating cell adhesion and migration, supporting tumor angiogenesis, and guiding tissue remodeling and repair [16–19]. Alterations in collagen fiber architecture, such as changes in fiber density and alignment, can create physical tracks for tumor cell invasion, making collagen organization a key determinant of malignant behavior [20, 21]. Thus, collagen fibers represent a promising yet underexplored source of imaging biomarkers for thyroid cancer risk assessment.

Multiphoton microscopy (MPM) uses near-infrared femtosecond laser excitation to generate two-photon excited fluorescence (TPEF) and second-harmonic generation (SHG) signals, enabling high-resolution, three-dimensional visualization of cellular and ECM architecture without exogenous dyes [22, 23]. Owing to its high sensitivity to collagen and elastic fibers as well as a range of endogenous fluorophores, MPM has become a powerful tool for characterizing collagen fiber distribution, stromal remodeling, and cell-matrix interactions within the TME [24, 25]. TPEF images provide morphological context by highlighting cellular structures, glandular architecture, and necrotic areas, thereby facilitating accurate ROI selection and anatomical localization of SHG signals. In quantitative analysis, feature extraction was performed exclusively on SHG images, given the specific sensitivity of SHG to collagen fibers. Pioneering works have established that SHG can visualize collagen and other structural proteins [26, 27], with subsequent studies extending its application to cancer diagnosis and prognosis [28, 29]. In contrast to conventional histopathological stains, such as Masson’s trichrome or picrosirius red staining, MPM is inherently label-free, thereby reducing systematic variability introduced by staining batch effects, protocol differences, and time-dependent signal loss due to fading or photobleaching [30, 31]. Previous studies have applied MPM to interrogate collagen architecture in multiple tumor types, including breast cancer, and have demonstrated that collagen fiber density and alignment are closely associated with tumor invasion and disease progression [32–34]. In thyroid neoplasms, SHG-based collagen characterization has been applied to differentiate benign from malignant nodules [35, 36] and to quantify capsular collagen organization [37]. In particular, Tokarz et al. revealed that malignant thyroid tissues exhibit a higher degree of ultrastructural disorder in collagen organization using SHG microscopy [36]. These findings collectively suggest that collagen architecture in the thyroid tumor microenvironment carries diagnostically relevant information, motivating our investigation into its prognostic value for recurrence risk.

This study aimed to investigate the relationship between the morphological and textural characteristics of collagen fiber within the TME of PTC and tumor recurrence. Recently, several studies have applied SHG imaging combined with deep learning for diagnostic classification of PTC, such as distinguishing benign from malignant lesions and detecting lymph node metastasis [38, 39]. While these works demonstrated the power of SHG-based artificial intelligence, they focused on diagnostic tasks rather than prognostic assessment. In contrast, our study addresses the distinct clinical need for postoperative recurrence risk stratification, employing a regularized linear model to identify interpretable collagen-derived features with explicit contributions to risk prediction. To this end, MPM was used to image tissue specimens from PTC patients, followed by quantitative image analysis, and an initial set of 142 features was extracted from the acquired images. The features were first screened using a dual-criterion approach based on Bonferroni-corrected P values and Cohen’s d effect sizes. Next, pairwise Pearson correlation coefficients were computed, and highly correlated features (|r| > 0.95) were removed to reduce redundancy and improve feature independence. Finally, least absolute shrinkage and selection operator (LASSO) regression was applied for further feature selection, yielding a final subset of features with the strongest predictive contribution. Collectively, these results deepen understanding of collagen-related microenvironmental remodeling in PTC and highlight collagen-derived features as potential biomarkers for postoperative surveillance and prognostic assessment. In addition, this study demonstrates the utility of MPM for TME characterization, supporting its broader applicability in oncologic research.

Sample preparation

This retrospective study included nine patients diagnosed with PTC at Fuzhou University Affiliated Provincial Hospital, China. Among them, five patients were classified as high-risk for recurrence and four patients as low-risk for recurrence. The classification was based on the 2015 American Thyroid Association (ATA) response-to-therapy criteria [8] and was independently performed by two experienced endocrinologists. High-risk patients were defined as those with evidence of biochemical or structural incomplete response during a minimum follow-up of 3 years, whereas low-risk patients were defined as those with excellent response (negative imaging and undetectable or low serum Tg). Surgically resected tissue specimens were collected from each patient and sectioned into 5-µm-thick slices for multiphoton imaging. Label-free multiphoton images were acquired from selected regions of interest (ROIs) within each specimen, with simultaneous collection of TPEF and SHG signals.

As the objective of this study was to establish associations between MPM-derived features and recurrence risk rather than to perform histopathological validation, ROIs were selected primarily based on SHG image quality and collagen content from the tumor capsule and the tumor-stroma interface, where collagen remodeling is most active. To minimize the impact of capsular heterogeneity, multiple ROIs (3–5 per patient) were selected from different areas of these regions to account for intra-specimen variability, and all ROIs were confirmed by experienced researchers. Importantly, three-dimensional (3D) imaging was performed on each ROI, with 30–40 images acquired at different depths, capturing the heterogeneity of collagen architecture across the tissue thickness. In total, 1720 SHG images were acquired, each with a resolution of 1,024 × 1,024 pixels, corresponding to a physical area of 465 × 465 μm. Of these, 956 images (55.6%) were from the high-risk group (5 patients) and 764 images (44.4%) from the low-risk group (4 patients). Given the relatively large and irregular thickness of the fresh specimens, hematoxylin and eosin-stained sections from the identical imaging regions were not available.

Multiphoton microscopy imaging system

The multiphoton microscopy imaging system used in this study integrated a Leica DMi8 inverted scanning microscope with an InSight X3 femtosecond laser that produces ultrashort pulses with a pulse width of <120 fs at a repetition rate of 80 MHz and provides a tunable excitation wavelength range of 680–1,300 nm. The DMi8, configured with a 25×/0.95 water-immersion objective, enables precise focusing of the excitation beam and efficient collection of emitted signals. The spatial resolution of this system is approximately 0.5 μm laterally and 2 μm axially, as estimated from diffraction-limited formulas based on the 25×/0.95 water-immersion objective and the excitation wavelength of 810 nm. The microscope supports automated operation and flexible configuration for high-quality imaging [40]. Signal generation relies on two-photon excitation, whereby fluorophores absorb two lower-energy photons nearly simultaneously to produce fluorescence emission. In this experiment, SHG signals were detected at 395–415 nm for collagen visualization, and TPEF was collected at 600–650 nm for tissue background and cellular structures.

Feature extraction

Quantitative analysis of collagen fibers was performed on the acquired SHG images. Briefly, each SHG image (1,024 × 1,024 pixels, corresponding to 465 × 465 μm) was converted to 8-bit grayscale (intensity range: 0–255) and preprocessed using a 3 × 3 median filter to suppress noise. Collagen fibers were then segmented using a fixed intensity threshold of 30, which was determined empirically from pilot experiments and kept constant across all images to ensure consistency. Based on this mask and the original grayscale images, a total of 142 features were extracted across four categories: network topology, orientation, intensity, and texture. Specifically, these comprise 8 network topology features computed using the FIRE toolbox, including metrics such as the fiber density, mean length, and inter-junction spacing, to characterize the complexity and connectivity of the collagen architecture. One additional orientation-related feature (anisotropy index) was quantified via Fourier transform to assess fiber alignment. 6 first-order intensity histogram features (e.g., mean, variance, and skewness) were extracted to describe the grayscale distribution. Orientation-related features were quantified via Fourier transform to assess fiber anisotropy. For texture analysis, a total of 80 second-order features were derived from the gray-level co-occurrence matrix (GLCM), including contrast, correlation, energy, and homogeneity, and were calculated at four directional angles (0°, 45°, 90°, 135°) with pixel distances of 1–5. For each direction and distance combination, the four GLCM statistics were computed independently and treated as separate features rather than averaged across directions. This design choice was intentional, as different collagen fiber orientations may reflect distinct aspects of stromal remodeling, and preserving directional information has been shown to be biologically relevant in the context of tumor microenvironment characterization. In addition, 48 local texture descriptors were obtained using multi-scale, multi-orientation Gabor filters across six orientations (0°, 30°, 60°, 90°, 120°, 150°) and four scales. Detailed definitions and calculation methods of all features have been described in previous studies [41–45]. Collectively, these features provided a comprehensive and fine-grained quantitative representation of microenvironmental alterations in collagen fiber organization.

Feature selection

A three-stage feature selection strategy was applied to the 142 collagen-derived features to identify an optimal subset with the greatest statistical significance, discriminative capability, and independence for predicting thyroid cancer recurrence status. In the first stage, all features were screened using dual criteria based on statistical significance and effect size. Specifically, the Wilcoxon rank-sum test was used to assess distributional differences between recurrence-risk groups, and Bonferroni correction was applied to control the false-positive rate arising from multiple comparisons. To ensure adequate clinical discrimination, features were additionally required to exhibit at least a moderate effect size (Cohen’s d ≥ 0.5). Of note, these two criteria served as sequential independent filters rather than a combined statistical test: the Wilcoxon test assessed statistical significance, while Cohen’s d subsequently ensured biologically meaningful differences. In the second stage, correlation-based redundancy removal was performed to improve interpretability and model stability. Pearson correlation coefficients were computed among the retained features; pairs with an absolute correlation coefficient exceeding the predefined threshold (|r| > 0.95) were considered highly redundant, and the feature with the smaller P value was retained to maintain a compact and relatively independent feature set. After initial screening and de-redundancy, the third stage employed regularization-based LASSO logistic regression for final selection. Prior to LASSO regression, all features were Z-score normalized to ensure comparability across different scales. LASSO enforces sparsity by shrinking coefficients of less informative features toward zero, thereby enabling automated feature selection. The optimal regularization parameter λ was determined via five-fold stratified cross-validation, with the data partitioned at the patient level rather than at the image level and stratified by recurrence risk group to maintain a balanced representation across folds. This patient-wise splitting strategy ensured that all images from the same patient were assigned exclusively to either the training or validation set in each fold, thereby preventing data leakage. The ‘lambda.1se’ criterion was adopted to select the most parsimonious model. As a consequence of the patient-level partitioning, the exact composition of each fold varied across iterations due to the random assignment of different patient combinations to each fold and the varying number of images contributed by each patient (depending on the number of ROIs and depth layers). Accordingly, approximately four-fifths of the images were used for training and one-fifth for validation in each iteration, rather than a fixed exact number. No independent hold-out test set was reserved, given the limited sample size. All hyperparameters were set to default values in the R glmnet package (α = 1), and features with the largest absolute coefficients in the final model were selected as the most predictive biomarkers. Through this sequential procedure, 82 features were retained after the first stage, 20 after the second stage, and 10 after the final LASSO selection.

Qualitative analysis of multiphoton images

MPM, which combined TPEF and SHG, was employed to visualize the tumor microenvironment in tissue samples from patients stratified into low- and high-recurrence-risk groups (Figure 1). In the low-risk group (Figures 1A–C), SHG images revealed collagen fibers exhibiting a wavy, interwoven architecture with irregular orientation, and the extracellular matrix architecture appeared relatively preserved, maintaining a more compliant fibrous network. By contrast, tissues from the high-risk group displayed marked stromal remodeling, characterized by thickened, straightened, and densely bundled collagen fibers that exhibited pronounced co-alignment and linearization (Figures 1D–F). These collagen alterations suggest that stromal morphological features captured by label-free multiphoton imaging may serve as potential imaging biomarkers for assessing the recurrence risk in PTC. Of note, the SHG images shown in Figure 1 were acquired from selected ROIs in each specimen. Given that multiple ROIs (3–5 per patient) were selected from different regions and 30–40 images were acquired at different depths for each ROI, the single-tile images presented are representative of the broader collagen architecture within each sample and adequately illustrate the typical morphological differences between the two risk groups.

Quantitative image analysis

A total of 142 quantitative features were extracted from the SHG images. The Wilcoxon rank-sum test was first applied to assess differences between recurrence-risk groups, combined with Bonferroni correction and Cohen’s d effect size to control for multiple comparisons and ensure adequate discriminative power. Features with statistically significant adjusted P values (P < 0.05) and at least moderate effect sizes (d ≥ 0.5) were retained, resulting in 82 features. Subsequently, redundant features were eliminated based on Pearson correlation analysis. When the absolute correlation coefficient between any two features exceeded r ≥ 0.95, the feature with the smaller P value was retained, yielding 20 features without strong intercorrelations. These 20 features were then used to construct a LASSO logistic regression model, with the optimal regularization parameter (λ) determined by five-fold cross-validation. Based on the absolute magnitude of the non-zero regression coefficients, the top ten features most strongly associated with thyroid cancer recurrence were identified, as shown in Figure 2 As listed in Table 1, GLCM_Correlation_90°_1 pixel and GLCM_Correlation_45°_2 pixel exhibited the largest LASSO regression coefficients, contributing most significantly to the model. From an imaging perspective, these GLCM correlation features primarily reflect the grayscale continuity and directional organization of collagen SHG signals at specific orientations and spatial scales, capturing the degree of alignment and anisotropy in collagen fiber architecture. This indicates that texture correlation, that is, the directional order of collagen fibers at vertical and oblique orientations, is highly relevant for distinguishing recurrence risk in thyroid cancer.

In addition, the fiber density feature calculated using the FIRE toolbox also showed a relatively large coefficient in the model, indicating that the collagen fiber number density is closely associated with the recurrence risk. The remaining multi-scale Gabor mean and variance features characterize the texture intensity and heterogeneity of the collagen network across different orientations and spatial scales, providing complementary information to the model. All ten features remained statistically significant after Bonferroni correction in the univariate Wilcoxon test. Figure 3 presents the overall distributions of these ten features between the high- and low-recurrence-risk groups, where red denotes high recurrence risk and blue denotes low recurrence risk; significant differences were observed for each feature. Specifically, the high-risk group exhibited significantly lower values for GLCM_Correlation_90°_1 pixel, GLCM_Correlation_45°_2 pixel, and fiber density compared with the low-risk group. In contrast, the high-risk group showed significantly higher values for Gabor_mean_30°_1 scale, Gabor_variance_0°_1 scale, Gabor_mean_150°_3 scale, Gabor_variance_150°_2 scale, Gabor_variance_90°_1 scale, Gabor_mean_90°_4 scale, and Gabor_variance_150°_1 scale, indicating a clear separation in the overall feature patterns between the two groups.

Kernel density estimation was further applied to the ten selected key features to generate probability density curves (Figure 4). By comparing the peak location, curve shape, and degree of overlap between the low-recurrence-risk group (blue) and the high-recurrence-risk group (red), the discriminative ability of each feature for recurrence outcome was visualized and the consistency of the preceding statistical tests and effect-size analysis were assessed. Overall, all features in Figure 4 exhibited evident distributional differences between the high- and low-risk groups. For GLCM_Correlation_90°_1 pixel and GLCM_Correlation_45°_2 pixel (Figures 4A,B), the density curves of the high-risk group (red) were shifted toward higher values, with more prominent peaks than those of the low-risk group (blue). This indicates stronger grayscale correlation at specific vertical and oblique orientations and spatial scales in high-risk tissues, consistent with greater continuity and directional organization of collagen fibers and reflecting stromal linearization and enhanced anisotropy. For fiber density (Figure 4C), the high-risk group showed a peak closer to the higher-value range with a broader distribution, suggesting increased collagen fiber density and greater spatial heterogeneity, whereas the low-risk group clustered within a lower-density range, indicative of a relatively looser stromal architecture. Regarding the multi-scale Gabor mean and variance features (Figures 4D–J), most features displayed an overall rightward shift of the high-risk curves, implying stronger texture responses and more complex local variations across multiple orientations and scales. This pattern is consistent with increased multiscale textural heterogeneity of the collagen network and aligns with stromal remodeling and fiber realignment. For some features, such as Gabor_mean_90°_4 scale and Gabor_variance_150°_1 scale (Figures 4I,J), the two curves overlapped substantially, suggesting limited discriminative power when considered individually; however, these features can still provide complementary information within the multivariable LASSO model, thereby improving the overall predictive performance. Taken together, the systematic differences between the red and blue curves in terms of the location, peak characteristics, and shape indicate that SHG textures in high-risk patients are more directional, ordered, and heterogeneous. These multi-orientation, multi-scale remodeling patterns constitute an important imaging basis for distinguishing recurrence-risk groups.

A predictive model was then built based on the ten selected features. Notably, while the distributions shown in Figures 3, 4 reflect original feature values to facilitate clinical interpretability, all features were Z-score normalized prior to model construction to ensure comparability across different scales. The model’s discriminative performance was evaluated using five-fold stratified cross-validation, and the predicted probabilities from the test sets of all folds were pooled to generate a receiver operating characteristic (ROC) curve, as shown in Figure 5. The ROC curve lies well above the chance line and bends toward the upper-left corner, with an area under the curve (AUC) of 0.959, indicating strong discrimination between the high- and low-recurrence-risk groups. The model achieved an accuracy of 91.0%, sensitivity of 92.3%, and specificity of 89.4%. It should be noted that this evaluation was performed at the image level; patient-level classification performance was not assessed due to the limited cohort size.

Thyroid cancer is among the most common endocrine malignancies, and its incidence has continued to increase globally in recent years. This increase is particularly evident for differentiated thyroid carcinoma, likely driven in part by advances in imaging and screening practices [46, 47]. Although the overall prognosis for most differentiated thyroid cancers is favorable, locoregional or distant recurrence remains relatively frequent, with reported recurrence rates of approximately 7%–30%, often occurring within 5–10 years after initial treatment. Recurrence can adversely affect long-term survival, quality of life, and the intensity of follow-up required [48–50]. Current clinical risk stratification for recurrence primarily relies on pathological staging, tumor size, lymph node involvement, and imaging assessments, such as ultrasonography and radionuclide-based examinations [51]. However, these approaches largely depend on postoperative information and experience-based judgment, and they are limited in their ability to capture microstructural alterations and TME remodeling. There is therefore a clear need for more objective and quantitative adjunctive tools for recurrence risk assessment. Notably, while recent studies have successfully applied SHG imaging combined with deep learning for diagnostic classification of thyroid tissues [38, 39], the prognostic question of postoperative recurrence stratification has remained unexplored.

To address this need, quantitative feature extraction and analysis were performed on MPM images of thyroid cancer tissue to characterize differences between low- and high-recurrence-risk groups at the levels of tissue microstructure and biochemical information in this study. MPM offers several advantages, including label-free imaging, high-resolution three-dimensional visualization, and strong sensitivity to structural components, such as collagen and elastic fibers. It enables acquisition of rich optical contrast without compromising tissue integrity, providing a basis for constructing quantitative “optical phenotypes.” In this study, a multi-stage selection pipeline was applied to the initial set of 142 SHG-derived quantitative features, sequentially incorporating univariate statistical testing with effect-size assessment, correlation-based redundancy removal, and LASSO logistic regression. This process ultimately retained ten features with significant intergroup differences and low redundancy as core variables for recurrence-risk discrimination. Boxplots and kernel density curves were then generated for these key features, visually illustrating group-wise differences in medians, dispersion, and overall distribution patterns, and providing graphical support for the robustness and consistency of the statistical screening results. Building on these features, a logistic regression–based predictive model was further developed and evaluated using five-fold stratified cross-validation. The resulting model achieved an AUC-ROC of 0.959, indicating excellent discrimination between patients at high versus low risk of recurrence. Importantly, MPM-derived quantitative features provide objective, reproducible, and readily standardizable information on tissue architecture and stromal status at the microscopic scale. In addition, the logistic regression framework is inherently interpretable, facilitating understanding of the relative contribution of individual features to recurrence risk. Regarding GLCM feature extraction, we acknowledge that averaging features across the four directions is standard practice for achieving rotation invariance. In our study, we intentionally preserved directional information because collagen fiber remodeling is inherently anisotropic and direction-specific features may capture biologically meaningful aspects of stromal remodeling that rotation-invariant features would obscure. However, we recognize that without orientation normalization, feature values may depend on sample orientation relative to the scanning coordinate system, potentially affecting reproducibility across different setups. Future studies could incorporate dominant orientation normalization procedures to enhance generalizability while preserving directional biological information. This interpretability is particularly important for our goal of identifying collagen-derived biomarkers that are clinically actionable. We therefore chose logistic regression with LASSO regularization over more complex classifiers for three reasons. First, given the limited sample size (9 patients), complex models such as SVM or deep learning are prone to overfitting, whereas LASSO performs automatic feature selection and shrinkage, making it better suited for small-sample, high-dimensional data. Second, our primary objective is not merely to achieve high classification accuracy but to reveal which specific collagen features drive recurrence risk, logistic regression provides explicit coefficients for this purpose. Third, while our previous studies have demonstrated the effectiveness of deep learning for diagnostic classification using SHG images, those models address a different clinical question and are less appropriate for the prognostic task addressed here. Finally, while we chose LASSO-logistic regression for its interpretability and suitability for small-sample data, we acknowledge that a systematic comparison with other classifiers (e.g., SVM) was not performed. Given the limited patient number (n = 9), such comparative experiments would be prone to overfitting and yield unstable results. Moreover, the primary goal of this study was feature discovery and biomarker identification rather than classifier optimization. Future multicenter studies with larger patient cohorts will systematically compare multiple modeling approaches to identify the optimal classifier for clinical application. Comparison with previous SHG-based thyroid studies, previous studies have demonstrated that SHG-derived collagen features can differentiate benign from malignant thyroid nodules and detect ultrastructural alterations in malignant tissues [35–37]. However, these studies focused on diagnostic discrimination. While in this study, we extended the application of SHG imaging to the prognostic domain by extracting collagen fiber features for postoperative recurrence risk stratification and employing an interpretable model that explicitly quantifies each feature’s contribution to risk prediction. This approach differs from and complements the aforementioned previous studies. This study applies SHG-based quantitative features specifically for recurrence risk stratification in thyroid cancer, thereby extending the utility of SHG imaging from diagnosis to prognosis.

Despite these promising results, several limitations of this study should be acknowledged. First, this was a single-center study with a relatively small sample size; the case composition and treatment strategies may introduce selection bias. Specifically, the limited number of patients precluded reliable performance evaluation at the individual patient level, the current cross-validation was conducted at the image level. Consequently, external validation in larger, multicenter cohorts is required to confirm the model’s generalizability and to enable patient-level assessment in future studies. Additionally, while we selected multiple ROIs from different regions of the tumor capsule and tumor-stroma interface and performed 3D multi-layer imaging to capture depth-wise heterogeneity, we acknowledge that manual ROI selection without fully systematic heterogeneity control may introduce bias. As recently demonstrated, the collagen framework of thyroid nodule capsules is heterogeneous at both the fiber architecture level [52] and the ultrastructural level [53]. Future studies will incorporate automated or semi-automated ROI selection strategies combined with whole-capsule mapping to more comprehensively address this heterogeneity. Second, the follow-up duration was limited, and some late recurrence events may not yet have occurred, potentially leading to conservative estimates of recurrence risk. Finally, this study primarily employed logistic regression, a classical linear model that facilitates interpretability but may not fully capture nonlinear relationships within a high-dimensional feature space. To build upon this work, future studies may incorporate more advanced machine learning or deep learning approaches and integrate MPM quantitative features with clinical, pathological, and molecular biomarkers through multimodal fusion to develop recurrence prediction models with improved accuracy while maintaining interpretability.

This study demonstrated that quantitative analysis of MPM-derived features can effectively reveal systematic microstructural differences between thyroid cancer patients with low versus high recurrence risk and enables construction of predictive models with strong discriminative performance. By providing a method to quantitatively assess collagen architecture and stromal remodeling, this approach establishes a new technical avenue and methodological basis for individualized recurrence risk assessment in thyroid cancer. Consequently, quantitative MPM analysis may serve as a valuable complement to existing risk-stratification frameworks, offering an objective guidance for postoperative surveillance strategies and personalize long-term management planning for patients with thyroid cancer.