Predictions for Central Lymph Node Metastasis of Papillary Thyroid Carcinoma via CNN-Based Fusion Modeling of Ultrasound Images

2.1 Data pre-processing

First, the PTC in the US image is manually segmented by a clinician with 10 years of experience using ITK-SNAP software, and the region of interest (ROI) of the PTC lesion is obtained automatically by programming. Then, the ROI is transformed into the Fourier frequency domain to obtain the low-frequency and high-frequency components [10]. After transformation, most of the energy is concentrated in the low-frequency component, while the detail is in the high-frequency component. The image is sharpened using a Gaussian high-pass filter [11], the transfer is used as shown in formula (1).

$H(u, v)=1-e^{-D^{2}(u, v) / 2 D_{0}^{2}}$     (1)

where, D₀ is the specified non-negative number, D(u, v) is the center distance from the point (u, v) to the filter.

To enhance the contrast and highlight the detailed information, the low-frequency and high-frequency components are weighted at the same time. However, the weighting coefficient of the low-frequency component must be less than that of the high-frequency one to ensure that the low-frequency component will not be lost, and the edge details of the high-frequency component are significantly enhanced. The transfer function of the high-frequency emphasis filter is shown as formula (2).

$H_{h f e}(u, v)=a+b H_{h p}(u, v)$     (2)

Finally, the inverse Fourier transformation is carried out in the spatial domain, and then histogram equalization is used to improve the brightness and contrast and highlight the features of the thyroid ROI. The tumor image patch is generated as a squared patch circumscribing the tumor as the conventional image patch for the tumor classifying RTF inputs. The generated tumor image patch can reflect the intra tumoral texture and tumor shape information to the RTF training.where, a is the offset, b is the multiplier, and H_hp(u, v) is the Gaussian high-pass filter transfer function.

The procedure of the pre-processing is shown in Figure 2. (a) is the original figure of PTC, and (b) is a map of the focus manually segmented by the clinician, (c) is the focus map after high frequency emphasis filtering and histogram equalization, and (d) shows the internal outline of the focus.

2.png

Figure 2. Data pre-processing

2.2 Feature extraction

For the image segmentation, this paper uses the gold standard of manual segmentation for the input of this method. The descriptions are summarized as 577 features. The features of this method are divided into nine categories: clinical statistical information, first-order features, gray level co-occurrence matrix features, gray level size zone matrix features, gray level run length matrix features, gray level dependence matrix features, neighboring gray-tone difference matrix features, shape features, and wavelet transforms. All imaging features can be divided into three categories [6, 12, 13].

2.2.1 First-order statistical characteristics

The first-order statistical feature is mainly the basic statistics describing the numerical distribution of pixels. As shown in formula (3) and (4) , X is used to represent the image, N is the number of pixels, X_i is the discrete gray intensity of i position, and $\bar{X}$ is the mean value of X, the gray histogram is the description of the surface layer of the image, and the features extracted from the histogram are invariant, such as scale and displacement and rotation, the histogram of each PTC is calculated from the ultrasound image, and the first-order statistical feature (mean, standard deviation (SD), skewness, and kurtosis) are calculated easily.

Moreover, the percentile-mean and percentile standard deviation (percentile-SD) are also derived from the first 50%, 25%, and 10% of the histogram curve.

Where, α represents the percentile (50%, 25% and 10%) and M is the number of pixels in the unselected (1-α) part.

$\operatorname{mean}_{\alpha}=\frac{1}{N-M} \sum_{i=M}^{N}\left(X_{i}-\bar{X}\right)^{2}$     (3)

$S D_{\alpha}=\sqrt{\frac{1}{N-M} \sum_{i=1}^{N}\left(X_{i}-\bar{X}\right)^{2}}$    (4)

2.2.2 Texture feature

Table 1 mainly defines the following general expressions, and then calculates the statistics.

The texture feature is mainly for extracting the texture feature of the gray level cooccurrence matrix (GLCM), defining the GLCM as a matrix with N_g×N_g and recorded as P(i, j, δ, θ). This matrix represents the second-order joint probability function of the image. Here, (i, j) is the number of times that a combination of i and j appear. These pixels are all with a distance of $\delta \text { in the } \theta$ direction.

Table 1. Expressions and their meanings

(1) Energy is the square sum of the elements in the gray co-occurrence matrix, which can reflect the uniformity of texture grayscale, reflect the texture thickness of the image distribution, and quantify the uniformity of the grayscale. It is shown in formula (5) .

$\text { energy }=\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{g}} P(i, j)^{2}$     (5)

(2) Contrast is a statistic of texture thickness which reflects the uniform gray distribution of the image. It can describe the numerical distribution and local variation of the matrix. It can also reflect the clarity and groove depth, and count the change of local strength. It is shown in formula (6).

$\text { contrast }=\sum_{i=1}^{N_{B}} \sum_{j=1}^{N_{g}}|i-j|^{2} P(i, j)$    (6)

(3) Correlation is used to calculate the gray correlation of a specified location, and the correlation is mainly to calculate the degree of similarity between elements in rows or columns. It is shown in formula (7).

$\text { correlation }=\frac{\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{g}} i j P(i, j)-\mu_{x}(i) \mu_{y}(j)}{\sigma_{x}(i) \sigma_{y}(j)}$     (7)

(4) Entropy describes the uniformity of the image and quantifies the randomness of the texture. The more uniform the image, the greater the entropy. It is shown in formula (8).

$\text { entropy }=\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{g}} P(i, j) \log _{2} P(i, j)$     (8)

(5) Homegeneity is a statistic that describes local changes and reflects the homogeneity of texture. It is inversely proportional to the contrast weight, and the larger the value is, the smaller the local change is, that is, the more uniform the local is. It is shown in formula (9).

$\text { homegenity }=\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{g}} \frac{P(i, j)}{1+|i-j|^{2}}$    (9)

2.2.3 Shape features

Shape features are expected to describe the shape-related features of the lesions, including the area, perimeter, and maximum diameter of the lesions. These features can get the objective description of the focus, which is the feature extraction of the surface information of the image.

N_pixels is the number of 1 in the ROI, and L_PS denotes the number of pixels of the focus. The calculating area is $A=N_{\text {pixels }} L_{P S}^{2}$, the circumference is $C=N_{\text {pixels }} L_{P S}$. Max diameter refers to the maximum Euclidean distance in pairs of pixels in the ROI region. The circumference to area ratio is the ratio of the perimeter to the area, $R_{C / A}=\frac{C}{A}$, the circular disproportion is $C D=\frac{C}{2 \pi R}$, and the roundness is $\text { roundness }=\frac{4 \pi A}{C^{2}}$.

All features are extracted from the original US-negative PTC images [14]. The relationship between lymph node status and imaging features is explored by machine learning.

Table 2 shows a part of the radiomics features extracted by machine learning.

Table 2. Radiomics features extracted by machine learning

2.3 Model reduction and the RF-based classifier

There are 577 groups of clinical features and image features. Because of the repeatability and redundancy of image group features, dealing with high-dimensional features requires more memory and processing power. Principal component analysis (PCA) can synthesize high-dimensional variables that may have correlation to linearly independent low-dimensional variables to alleviate dimensional disaster problems. At the same time, dimensionality reduction can minimize information loss while compressing data.

The PCA is employed to reduce the dimensions of the data, and can keep only the terms corresponding to the K largest eigenvalues. Hence, this obtains a new feature vector consisting of eigenvectors of principal components. The final data computed using this feature vector and the mean adjusted original input data is obtained as shown in formula (10).

$F_{\text {data }}=R F_{\text {vector }} \times R D_{\text {adjust }}$    (10)

where, RF_vector is the matrix in which eigenvectors in the columns are transposed and RD_adjust is the mean adjusted input data. F_data is the characteristic matrix after PCA dimension reduction. The obtained subspace is spanned by the orthogonal set of eigenvectors which reveal the maximum variance in the data space.

Using PCA mapping high-dimensional data into low dimensional data reduces the calculation cost of NIDS and improves the efficiency of the analysis. Here, PCA has been used for dimensionality reduction of the forty-two dimensions and the output of the PCA method provides a set of features that are the linear combination of the original set of features. PCA provides RF with the features that allow for efficient classification from an informative view point. Then, the RFclassifier [15, 16] is used in this method and 100 trees are employed in each RF. Data enhancement is performed on the training set to reduce deviation and overfitting.

Different features according to the type of input image can be obtained because of the different training methods described above. The model trained with tumor image patches mainly learns the correlation between texture information and clinical information, while the VGG16 model trained with context image patches mainly learns the relationship of context information.

4.1 Data selection

This retrospective research protocol has been approved by the Ethics Review Board of the Department of Ultrasound at Sun Yat-sen University Cancer Center. Written informed consent from all participants has been obtained. All of the procedures are performed in accordance with the Declaration of Helsinki and relevant policies in China.

The participants include 906 consecutive patients who underwent near, sub-, or total thyroidectomy, and have been pathologically confirmed with PTC by the Department of Head and Neck in the hospital between January 2016 and May 2019.The inclusion criteria for the thyroidectomy are: (a) patients ≥18 years old, and (b) no clear CLN identified from US images. The exclusion criteria are: (a) patients who previously underwent a neck operation, (b) patients diagnosed with other types of thyroid tumor at the same time, (c) patients with a previous history of radiation therapy, and (d) patients who are unavailable to be completely evaluated by US. Patient enrolment details are shown in Figure 4. In the center, surgical procedures are performed on patients based on the recommendation of the National Comprehensive Cancer Network and the American Thyroid Association. Prophylactic CLND is commonly performed in the VI compartment (neck level VI) according to the institutional protocol, regardless of the clinical evidence of CLNM. Prophylactic CLND is typically perforated.

Clinical data, such as age, gender, and other basic patient information, are collected from the electronic medical records in the registry. Pathologic information is gathered from postoperative pathology reports in the electronic medical record. All diagnoses are rendered and reported by pathologists with 7–20 years of experience.

A comprehensive neck US examination is preoperatively performed on all the patients in the supine position, with their neck extended, using a 5–18 MHz linear array transducer machine (iU22, Philips Medical Systems; Acuson Sequoia 512, Siemens Medical Solutions; LOGIQ S8 and E9, GE Medical Systems). These tests are performed by board-certified doctors with 7–20 years of experience, specializing in head and neck imaging. The doctor who performed the US examination prospectively records the US features of the thyroid nodule and CLN status. If more than one nodule under suspicion of malignancy has been found in the thyroid gland, the maximum diameter of the most suspicious lesion is recorded and included in the data analysis.

4.png

Figure 4. Patient enrolment details

4.2 Results analysis

4.2.1 Nodule characteristics

All 906 patients are confirmed to have PTC by surgical pathology in the Sun Yat-sen University Cancer Center from January 2016 to May 2019, among which 322 cases are CLNM positive and 584 cases involved metastasis. The PTC diagnosis is based on the final pathology results. The 634 cases (about 70%) in the training set have a median age of 42.76 years (in the range of 18-74), and include 191 males (30.13%) and 443 females (69.87%), where 228 cases (35.96%) are CLNM positive. Further, 196 cases in the validation set (about 30%) have a median age of 44.30 years (in the range of 18–79), and include 93 males (34.19%) and 179 females (65.81%), where 94 cases (34.56%) are CLNM positive. No significant difference is found in demographic profiles between the two groups.

Eight features are listed based on clinical and US features. Univariate analysis shows that CLNM is related to sex, age, tumor size, microcalcification, extra-thyroid invasion, shape, vascular invasion, and aspect ratio. These eight features are identified as the key factors and are discussed considering the manual extraction of features (Table 4; P<0.001).

Table 4. Relationship between clinical and US characteristics

4.2.2 Comparison of the predict models

In the validation set, 272 cases are input into the RTF model, the fine-tuned VGG16 model, and the combined model. Table 5 shows the CLN prediction results of the three comparison methods. In the validation set, the VGG16 model shows slightly lower accuracy than the RTF model (78.31% and 81.92%, respectively), although the sensitivity of VGG16 is significantly higher than that of the RTF model. This may be because of the fact that the texture and shape of the tumor are important criteria for predicting CLNM, in addition to the tissue around the tumor. However, compared with the RTF and VGG16 models, the combined model shows an increase sensitivity of13.83% and 3.19%, respectively, and an increase in ACC of 1.17% and 4.78%, respectively.

Table 5. Performance comparison of machine-learning prediction models based on US features

Figure 5 shows the receiver operating characteristic curves(ROC) obtained using the three different models. The AUC of the VGG16 model is slightly lower than that of the RTF model. However, as context and intra-tumor information are complementary for the prediction, the combined model that includes both of these data shows an increase in overall performance, with an AUC value of 0.9228. Hence, the combination of the context VGG16 model and RTF model has a high potential for predicting CLNM in US-negative patients with PTC. The results show that the fusion of the model has the potential to assist CLNM detection.

5.png

Figure 5. Receiver operating characteristic curves obtained using the three different models

4.2.3 Visualization effect

As shown in Figure 6, (b) map is the visualization of (a) map. It can be seen from the figure that the deep convolution network has noticed the focus area and marked it clearly. This shows that the deep convolution network has known the abnormal characteristics of malignant tumors and used them as an important basis for the prediction of central lymph node metastasis. (b) map pays attention not only to the location of the lesion, but also to the edge of the thyroid gland. This shows that the surrounding tissue information is also the basis for the prediction of central lymph node metastasis, and validates that the model studied in this paper not only pays attention to the texture and shape of the focus, but also is sensitive to context information and surrounding tissue information.

6.png

Figure 6. Feature visualization

4.3 Discussion

PTC is one of the most common cancers diagnosed in the world, and its main life-threatening manifestations are lymph node metastasis and extra-thyroid spread [4, 23, 24]. The early diagnosis and treatment of CLNM is the most significant means of reducing the mortality and recurrence rate of patients with PTC [25, 26]. US is an effective tool for precise localization of PTC lesions, however, most radiologists do not have enough training or experience to read US images [4]. Additionally, most of the previous studies are limited by the relatively small number of patients [27], lacking post-operative diagnosis results as a standard [25, 28]. Furthermore, most of the studies use human-defined image features to build a diagnostic model to predict CLNM, which inevitably introduces subjective and empirical bias [19, 29-31]. Hence, the machine learning-based combined model is used to assist radiologists in segmenting thyroid nodules and predicting CLNM of PTC.

This paper performs fine-tuning based on the concept of transfer learning to overcome the problem of insufficient training data. After transferring parameters, the model continued to fine-tune by training with the US images to achieve better performance. As with previous studies, it is helpful to perform the fine-tuning of annotated nature image datasets (ImageNet) on a large scale [17]. The reasons for this might include: (1) Training a VGG16 model on a small dataset from scratch will reduce VGG16 generalizability, often resulting in overfitting; (2) The fine-tuned VGG16 with initialized weights is more robust than training VGG16 from scratch with random initialization; (3) Fine-tuning will reduce training time and get a better result.

In general, for small data medical images, the architecture of the selected model is too deep and it is easy to lose low-level information. The model employs an optimized feature extraction technique by hierarchically integrating multiple feature layers from VGG16 to incorporate both low and high-level information from US images.

Moreover, according to the experimental results shown in Table 5, the accuracy of the RTF model is higher than that of the VGG16 but lower than the sensitivity of the fine-tuned VGG16. The primary reason might be that direct training with limited training data is not enough to optimize the parameters of the VGG16. Therefore, the automatic feature extractor of the VGG16 could not produce features more effectively than the handcrafted method, since the small dataset lacks diversity and representative ability, resulting in the VGG16 model overfitting and performance degradation [8]. Later research can collect more data (ex. a thousand cases of each type) with varying diversity and representative ability that might enhance the VGG16 method and outperform the handcrafted method.

Besides evaluating the performance of the VGG16 and RTF models, the research also evaluates the performance of a combined model. According to the experimental results, the performance of the combined model is better than VGG16 and RTF models. This is because of the following main reasons. (1) RTF model with tumor image patches primarily learns the correlation between texture information and clinical information, while VGG16 model trained with context image patches primarily learns the relationship of context information. (2) VGG16 model extracts more detailed features, including low-level to high-level features from the images. (3) The handcrafted RTF model effectively extracts images and clinical features, and the RF classifier reduces the risk of overfitting and improved the robustness of the model by averaging decision trees.

A validation set of 272 cases are used to prove that the combined model diagnostic system has a satisfactory diagnostic performance. Compared with other models, the predictive ability of the combined model is superior. The VGG16 uses the visual attention map to infer the important areas related to the transfer of the US image. In addition, during the development of the US images, various sources of US noise, such as blood vessels, the trachea, the esophagus, fat, and muscle tissue, interfere with the detection of the real PTC nodule. Therefore, it is most effective to extract only the PTC nodule region for the training set, and that significantly improves the performance of the diagnosis model. The study develops and validates a model based on DL that is designed to predict CLNM in US-negative patients with PTC. The research expects that the model could be implemented as part of an improved diagnostic system in the future that could provide a simple solution for the early diagnosis of CLNM.

Although this is a single-center retrospective research with a small sample size, the diagnosis and treatment processes are standardized and unified. However, the images in the dataset only come from the same type of machine. Therefore, the robustness of the proposed system using images from different types of machines is yetun confirmed. The acquisition of images from different types of machines is necessary for fine-tuning the model to confirm its robustness. Additionally, the model currently has good prediction performance, however, the research has not performed external validation. In the future, the model requires more prospective and multicenter studies to confirm the accuracy of the predictions and improve the system. Despite these shortcomings, the machine learning diagnosis system is an effective predictor of PTC in patients using preoperative data. Hence, the research proposes that this simple and reliable method could be useful in a clinical setting to improve diagnostic sensitivity and assist in the diagnosis and treatment of PTC.