Tumor Prognostic Prediction of Nasopharyngeal Carcinoma Using CT-Based Radiomics in Non-Chinese Patients.

PURPOSE: We aimed to construct predictive models for the overall survival (OS), progression-free survival (PFS), and distant metastasis-free survival (DMFS) for nasopharyngeal carcinoma (NPC) patients by using CT-based radiomics.
MATERIALS AND METHODS: We collected data from 197 NPC patients. For each patient, radiomic features were extracted from the CT image acquired at pretreatment via PyRadiomics. Feature selection was performed in two steps. First, features with high inter-observer variability based on multiple tumor delineations were excluded. Then, stratified bootstrappings were performed to identify feature combinations that most frequently achieved the highest (i) area under the receiver operating characteristic curve (AUC) for predicting 3-year OS, PFS, and DMFS or (ii) Harrell's C-index for predicting time to event. Finally, regularized logistic regression and Cox proportional hazard models with the most frequently selected feature combinations as input were tuned using cross-validation. Additionally, we examined the robustness of the constructed model to variation in tumor delineation by simulating 100 realizations of radiomic feature values to mimic features extracted from different tumor boundaries.
RESULTS: The combined model that used both radiomics and clinical features yielded significantly higher AUC and Harrell's C-index than models using either feature set alone for all outcomes (p < 0.05). The AUCs and Harrell's C-indices of the clinical-only and radiomics-only models ranged from 0.758 ± 0.091 to 0.789 ± 0.082 and from 0.747 ± 0.062 to 0.767 ± 0.074, respectively. In comparison, the combined models achieved AUC of 0.801 ± 0.075 to 0.813 ± 0.078 and Harrell's C-indices of 0.779 ± 0.066 to 0.796 ± 0.069. The results showed that our models were robust to variation in tumor delineation with the coefficient of variation ranging from 4.8% to 6.4% and from 6.7% to 9.3% for AUC and Harrell's C-index, respectively.
CONCLUSION: Our results demonstrated that using CT-based radiomic features together with clinical features provided superior NPC prognostic prediction than using either clinical or radiomic features alone.

1 Introduction

Nasopharyngeal carcinoma (NPC) is one of the most malignant head and neck cancers worldwide and is endemic in Southeast Asia and Southern China (1, 2). According to the International Agency for Research on Cancer, there were about 130,000 new cases of nasopharyngeal carcinoma in 2020 (3). The mainstay treatment of NPC is intensity-modulated radiation therapy (IMRT) with chemotherapy. However, the treatment strategies, e.g., concurrent chemoradiotherapy (CCRT), CCRT followed by adjuvant chemotherapy, or induction chemotherapy followed by CCRT need personalized data for therapeutic planning (4).

In NPC, clinical staging according to the American Joint Committee on Cancer (AJCC) is conventionally used to guide the optimal treatment and determine cancer prognosis (5–8). Until recently, plasma Epstein–Barr virus (EBV) DNA has been widely used for early detection, prognostication, and monitoring of treatment response of NPC (8–12). Nevertheless, the EBV DNA level was undetectable up to 40% in non-Chinese series (13). Therefore, an additional effective biomarker is needed to improve the prognostic performance.

Medical imaging is a routine practice in oncologic patient management for tumor diagnosis and staging, treatment planning, and response monitoring. Computed tomography (CT) and magnetic resonance imaging (MRI) are common modalities in NPC imaging. Tumor appearance in radiographic images and demographic features such as patients’ age, tumor stage, and performance status were used to select an appropriate treatment. However, conventionally, the radiologic tumor description is provided qualitatively by focusing on tumor size and anatomical extension without considering intratumoral heterogeneity, which is predictive of the prognostic outcomes (14, 15). Radiomics, an emerging technique for tumor characterization, is the process of converting medical images into minable high-dimensional data (16–19). It refers to the extraction of quantitative features, so-called radiomic features, which describe detailed tumor characteristics, from whole tumor mass. Examples of radiomic features include tumor size, geometry, voxel intensity, and texture patterns. Several studies showed potential of radiomic features in head and neck cancer prognostic prediction. Aerts et al. (20) found that the radiomics model was able to capture intratumoral heterogeneity and was significantly associated with the gene-expression profile pattern. Moreover, results from various studies showed high correlation between radiomic features and prognostic outcomes, such as in head and neck cancer (21, 22), esophageal cancer (23, 24), and nasopharyngeal cancer (25, 26).

In NPC, several studies on radiomics also showed great promise for prognosis. Du et al. (25) used pretreatment contrast-enhanced T1- or T2-weighted MRI to construct predictive models for 3-year disease progression after IMRT treatment. The results showed that combining radiomics with other clinical features including age, sex, and TNM staging yielded higher prognostic performance than using clinical features alone. Similarly, Ming et al. (26) found that combining clinical features with MRI-based radiomics yielded higher performance. Some studies integrated pretreatment Epstein–Bar virus (EBV) DNA information as part of clinical features (27–29). The combined model of clinical and radiomic features yielded the highest prognostic performance in these studies as well.

In countries with limited resources, CT scans are more widely used and accessible than MRI because they are typically less expensive. However, while there were several research studies that used MRI data of non-metastasis NPC patients, none of them investigated the predictive value of CT-based radiomics in NPC. Therefore, in our study, we used CT-based radiomics to build the predictive model for NPC. We hypothesized that CT information, like MRI’s, should also improve the predictive performance. Moreover, if successful, our model would be beneficial for centers without an MRI simulator.

Specifically, we developed prognostic models for NPC based on conventional clinical features, pretreatment CT radiomic features, and the combination of the two. We compared their prognostic performances in terms of overall survival (OS), progression-free survival (PFS), and distant metastasis-free survival (DMFS) at 3 years as well as time-to-event outcomes.

2 Materials and Methods

2.1 Dataset

Data from newly diagnosed NPC patients at King Chulalongkorn Memorial Hospital, Thai Red Cross Society, Bangkok, Thailand, between October 2010 and September 2015 were collected. This group of patients was part of the patients who were enrolled in a previously reported randomized study comparing the IMRT technique between sequential versus simultaneous integrated boost technique (30). This study was approved by the Institutional Review Board of Faculty of Medicine, Chulalongkorn University, Bangkok, Thailand (IRB no. 745/61). The inclusion criteria of this study were as follows: (a) newly diagnosed with NPC patients, (b) no evidence of distant metastasis, (c) at least 3-year follow-up time, (d) underwent CT and MRI simulation, (e) received the IMRT with chemotherapy, and (f) available pretreatment plasma Epstein–Barr virus (EBV) DNA level. The clinical data included sex, age, and tumor staging. Patient were restaged according to the 8th edition of the American Joint Committee on Cancer (AJCC) TNM staging system (5). Patients with more than T1 or positive nodal disease received IMRT 70 Gy in 33–35 fractions with a concurrent chemotherapy by weekly cisplatin 40 mg/m2 for a maximum of seven cycles. Cisplatin 80 mg/m2 and 5-fluorouracil (5-FU) 1,000 mg/m2/24 h were given as adjuvant chemotherapy over a 96-h continuous infusion at 4-week intervals for three cycles.

2.2 Region of Interest Segmentation

For each patient, the gross tumor volume (GTV) of the primary tumor in the nasopharynx was contoured on CT images by the treating radiation oncologist. When a GTV was drawn, the co-registered MR image was also presented. Of 197 patients, thirty patients were randomly selected and their GTVs were additionally contoured by two radiation oncologists. This created a multiple tumor delineation dataset (with 3 GTVs for each of these 30 patients) for testing inter-observer variability.

2.3 Model Construction

2.3.1 Radiomic Feature Extraction

To ensure the spatial consistency in the radiomics analysis, all CT images were resampled into 0.5 × 0.5 × 3 mm3 voxels. A total of 842 radiomic features per GTV were calculated via Pyradiomics version 2.0.0 (31). The extracted features were classified into four feature classes (20), including the shape-based class, first-order intensity class, texture-based class, and wavelet-based class.

2.3.2 Radiomic Feature Selection

After feature extraction, the number of features was reduced via a feature selection process to prevent model overfitting. To select robust and informative features, we performed inter-observer variability test and univariate analysis as follows.

2.3.2.1 Inter-Observer Stability Test

The multiple delineation dataset was used for the inter-observer variability test. Good radiomic features should be consistent regardless of the radiation oncologist who drew the contour. The intra-class correlation (ICC) analysis was used to assess the correlation of features from multiple regions of interest (ROIs). According to a study of Koo et al. (32), features whose ICC values were less than 0.5 were excluded from further analyses.

2.3.2.2 Feature Performance Analysis

To assess the feature’s performance, a univariable analysis using logistic regression and Cox proportional hazard regression analysis for the binary outcome [having survived (OS), being disease-free (PFS), and being DMFS at 3 years] and survival outcomes, respectively, was performed. We repeatedly sampled 80% of the dataset with stratification for 100 times to calculate the area under the receiver operating characteristic (ROC) curve (AUC) and Harrell’s C-index. For each repeat, each outcome, and each radiomic feature class (20), the top single feature that yielded the highest AUC or C-index on the training set was selected. The most frequently selected set of 4 radiomic features was chosen for subsequent radiomics model construction.

2.3.3 Clinical Feature Selection

All clinical features including age, sex, T-stage, N-stage, overall stage (AJCC 8th edition), and pretreatment plasma EBV DNA level (cutoff = 2,300 copies/ml) (13) were initially included into the multivariate logistic regression or Cox proportional hazard models. Then, backward feature elimination was performed to repeatedly remove unimportant features with Wald test p-values greater than 0.1 until no feature can be removed. This process was repeated 100 times by sampling 80% of the dataset with stratification. The most frequently retained set of clinical features were selected for subsequent clinical-only model constructions.

For the combined models which used both clinical and radiomic features, backward feature elimination of clinical features was performed similarly as described, with two additional constraints: the set of 4 radiomic features selected above were always retained in the models, and at least one clinical feature must remain in each model. The most frequently retained set of clinical features were selected for subsequent combined model constructions.

2.3.4 Prognostic Model Construction

For binary outcomes, which were 3-year OS, 3-year PFS, and 3-year DMFS, logistic regression models were built. L1 (Lasso) or L2 (Ridge) regularizations with inverse strength ranging from 0.001 to 10 were considered. For time-to-event outcomes, Cox proportional hazard regression models were used. Elastic net regularizations with strength ranging from 0.001 to 10 and L1/L2 ratio ranging from 0.001 to 1 were considered. The best hyperparameters were selected based on AUC or Harrell’s C-index from 20 realizations of 5-fold cross-validation (for a total of 100 repeats). Regularized model developments were performed in Python using the LogisticRegression module of scikit-learn library and the CoxPHFitter module of lifelines library (33, 34). The overall process of model construction is shown in the .

Figure 1

Overall process in this study.

2.4 Model Comparison (Statistical Analysis)

The sign test in the scipy package with default settings was used to compare the Radiomics, Clinical, and Combined models and the model that used only tumor volume. The Benjamini/Hochberg procedure was used for multiple-comparison correction. We also constructed the Kaplan–Meier plot for 100 repeats of the test set, where patients were classified as high risk if their predicted scores were higher than the median and vice versa. All statistical analyses, except sign test and log-rank test which were analyzed in Python, were performed using STATA version 15.0 (StataCorp LLC, Texas, USA) (35).

2.5 Testing Model Robustness to Variation in Radiomic Feature Extraction

To investigate the model’s robustness to the variation in tumor delineation by different radiation oncologists, we generated 100 realizations of a normal random vector whose mean was equal to the radiomic features extracted and whose variance was equal to the sample variance calculated on the multiple-delineation dataset. For each outcome (OS/PFS/DMFS and binary/time-to-event), we computed the mean and standard deviation of the corresponding AUC or Harrell’s C-index across 100 realizations ().

Figure 2

Process of model robustness to radiomic feature variation due to segmentation.

3 Results

3.1 Patient Characteristics

A total of 197 patients had non-metastatic nasopharyngeal cancer, and approximately 80% of them were classified as undifferentiated carcinoma (). For tumor (T) classification, the proportion of patients was 22.3%, 37.1%, 24.4%, and 16.2% for T1, T2, T3, and T4, respectively. For node (N) classification, the proportion was 3.0%, 24.4%, 49.2%, and 23.4% for N0, N1, N2, and N3, respectively, while for stage grouping, the proportion was 0.5%, 14.2%, 49.2%, and 36.1% for stages I, II, III, and IVA, respectively. Fifty percent of patients had pretreatment plasma EBV DNA level ≥2,300 copies/ml. The median age was 50 years (interquartile range: 43 to 56), and 79.2% were male. The median follow-up time was 4.9 years.

Table 1

Patient characteristics.

Characteristics	n = 197 patients (%)
Median age (years) (IQR)	50 (43 to 56)
Sex
Male	156 (79.2)
Female	41 (20.8)
T classification
T1	44 (22.3)
T2	73 (37.1)
T3	48 (24.4)
T4	32 (16.2)
N classification
N0	6 (3.0)
N1	48 (24.4)
N2	97 (49.2)
N3	46 (23.4)
Stage group
I	1 (0.5)
II	28 (14.2)
III	97 (49.2)
IVA	71 (36.1)
Pretreatment plasma EBV DNA level
Undetectable or < 2300 copies/ml	100 (50.76)
≥2,300 copies/ml	97 (49.24)
Median EBV value (copies/ml) (IQR)	7,795 (3150 to 18000)
Pathologic classification
Undifferentiated carcinoma	159 (80.71)
Differentiated non keratinizing carcinoma	37 (18.78)
Poorly differentiated squamous cell carcinoma	1 (0.51)

3.2 Radiomics Model

Out of 842 features, we found that 384 features passed the ICC criterion for the inter-observer stability test. Across 100 bootstrap repeats, the set of four features (one for each class) which was most frequently selected for each type of outcome is listed in . Interestingly, “original_firstorder_Uniformity” from the first-order class was consistently selected in all binary and time-to-event outcomes. For the texture class, “original_glfm_DependenceNonUniformity” was selected for OS outcomes while “original_glrlm_GrayLevelNonUniformity” was selected for PFS and DMFS outcomes.

Table 2

Most frequently selected features of OS, PFS, and DMFS in each class on NPC patients.

Statistical analysis	Feature group	Feature name
		OS	PFS	DMFS
Logistic regression (binary outcome)	Shape	original_shape_MajorAxisLength	original_shape_MajorAxisLength	original_shape_MajorAxisLength
	First order	original_firstorder_Uniformity	original_firstorder_Uniformity	original_firstorder_Uniformity
	Texture	original_gldm_DependenceNonUniformity	original_glrlm_GrayLevelNonUniformity	original_glrlm_GrayLevelNonUniformity
	Wavelet	wavelet-LHL_gldm_LargeDependenceEmphasis	wavelet-LHL_glrlm_GrayLevelNonUniformity	wavelet-HHL_ngtdm_Busyness
	Clinical	age	age	age
		plasma EBV DNA level	plasma EBV DNA level	plasma EBV DNA level
Cox regression (time-to-event outcome)	Shape	original_shape_MajorAxisLength	original_shape_SurfaceArea	original_shape_SurfaceArea
	First order	original_firstorder_Uniformity	original_firstorder_Uniformity	original_firstorder_Uniformity
	Texture	original_gldm_DependenceNonUniformity	original_glrlm_GrayLevelNonUniformity	original_glrlm_GrayLevelNonUniformity
	Wavelet	wavelet-LHL_glrlm_RunVariance	wavelet-HLL_glrlm_RunLengthNonUniformity	wavelet-HLL_glrlm_RunLengthNonUniformity
	Clinical	age	age	age
		plasma EBV DNA level	plasma EBV DNA level	plasma EBV DNA level
				N stage (8th edition)

3.3 Clinical Model

The results from 100 bootstraps showed that the models for both binary and time-to-event outcomes selected the same feature sets. For OS and PFS outcomes, frequently selected clinical features were age, T-stage, and pretreatment EBV DNA level. For DMFS, age, T stage, N stage, and pretreatment EBV DNA level were selected.

3.4 Combined Model

We found that age and pretreatment EBV DNA level were commonly selected across all outcomes in the combined models. Furthermore, N stage was also selected in the model for the time-to-event DMFS outcome ().

3.5 Classification and Evaluation

The results of multivariate analyses showed that the combined model yielded the best performance in every outcome as compared to the clinical and the radiomics models ().

Figure 3

Boxplots of AUC values for (A) OS, (B) PFS, (C) DMFS, and C-indices for (D) OS, (E) PFS, and (F) DMFS of the clinical, radiomics, and combined models.

For binary outcome predictions, the combined model yielded AUCs of 0.817 ± 0.016 to 0.844 ± 0.018 and 0.801 ± 0.075 to 0.813 ± 0.078 on the training and test sets, respectively. In comparison, the AUCs of the clinical and radiomics models ranged from 0.769 ± 0.023 to 0.804 ± 0.024 on the training set and from 0.758 ± 0.091 to 0.789 ± 0.082 on the test set.

Harrell’s C index from the combined model in predicting time-to-event outcomes ranged from 0.791 ± 0.015 to 0.815 ± 0.017 and from 0.779 ± 0.066 to 0.796 ± 0.069 on the training and the test sets, respectively. These were higher than those from the clinical (0.767 ± 0.016–0.781 ± 0.017 and 0.756 ± 0.066–0.767 ± 0.074) and the radiomics models (0.747 ± 0.014–0.759 ± 0.018 and 0.747 ± 0.062–0.749 ± 0.077) on both the training and test sets. The performances of the combined model were significantly higher than those of the clinical or the radiomics model for all outcomes (p < 0.05).

The tumor volume yielded AUCs of 0.746 ± 0.024 to 0.773 ± 0.017 and 0.747 ± 0.097 to 0.776 ± 0.083 on the training and test sets, respectively. Harrell’s C index from the tumor volume ranged from 0.731 ± 0.022 to 0.744 ± 0.014 and from 0.734 ± 0.099 to 0.746 ± 0.064 on the training and test sets, respectively. The tumor volume yielded significantly lower AUC and Harrell’s C-index than those of the combined model for all outcomes (p < 0.001). Comparison of the tumor volume with the radiomics model in the test sets showed no significant difference between them in terms of AUC in OS, PFS, and DMFS predictions, and in terms of C-index in DMFS. For the other outcomes (OS and PFS predictions based on Harrell’s C-index), the radiomics model yielded higher performance metrics than the tumor volume with p < 0.05.

The Kaplan–Meier plot for the test sets showed a significant separation (log-rank test p < 0.001 in all 100 repeats) between the survival curves of the high-risk and low-risk groups, which were defined by the median of predictions as described in Section 2.4 ().

Figure 4

Kaplan–Meier plot for OS, PFS, and DMFS stratified by the median of the combined model’s score.

3.6 Model Robustness to Variation in Radiomic Feature Extraction

After adding uncertainty to the radiomic feature values and reconstructing the models, we found that the mean and standard deviation among 100 computed AUC values were 0.886 ± 0.043, 0.818 ± 0.052, and 0.781 ± 0.047 in OS, PFS, and DMFS, respectively. On the other hand, Harrell’s C-indices were 0.860 ± 0.064, 0.757 ± 0.070, and 0.819 ± 0.055 in OS, PFS, and DMFS, respectively. The corresponding coefficient of variation ranged from 4.8% to 6.4% and from 6.7% to 9.3% for AUC and Harrell’s C-index, respectively.

4 Discussion

Clinical data such as TNM staging (5), age, and blood test (EBV concentration level) were traditionally used for predicting prognosis in NPC patients and determining treatment strategy. Based on Lertbutsayanukul et al. (13), the optimal cutoff value for plasma EBV DNA for predicting OS, PFS, and DMFS was 2,300 copies/ml. Hence, the same cutoff was applied in our analysis. This transformed the pretreatment plasma EBV DNA level into a discrete variable. In our study, we constructed NPC tumor prognostic models using various combinations of CT-based radiomic features and clinical features. The corresponding AUC and Harrell’s C-index of OS, PFS, and DMFS outcomes for each constructed model were used to compare performances and identify the best models.

The combined model incorporating both clinical and radiomic features was superior to the clinical or radiomics model for OS, PFS, and DMFS predictions. This demonstrated the predictive value of quantitative information extracted from medical images when used in conjunction with conventional clinical features. Beyond the volumetric measurement, radiomic features also described tumor characteristics, especially intratumor heterogeneity, which has been shown to be related to tumor prognosis.

In our predictive models, we found that the most frequently selected radiomic features from the shape class were “original_shape_MajorAxisLength” and “original_shape_SurfaceArea.” Features in the shape class represent the size and shape of the ROI, independently of the gray-level intensity distribution. Both “original_shape_MajorAxisLength,” which was calculated by using the largest principal component, and “original_shape_SurfaceArea,” which indicated the total surface area of the ROI, were related to the tumor size. From the first-order class, the most frequently selected radiomic feature was “original_firstorder_Uniformity,” which indicated homogeneity of the ROI intensity. For the texture and wavelet classes, the most frequently selected radiomic features were “gldm_DependenceNonUniformity,” “glrlm_GrayLevelNonUniformity,” and “glrlm_RunLengthNonUniformity.” These features measure similarity of gray-level patterns in the ROI. For “gldm_DependenceNonUniformity” which measured the similarity in dependency throughout the image, a lower value indicated more heterogeneity in the ROI and poorer prognosis. The “glrlm_GrayLevelNonUniformity” feature measured the similarity in gray-level intensity in the image, with a lower value indicating more tumor homogeneity. Similarly, the “glrlm_RunLengthNonUniformity” feature measured the similarity in run lengths throughout the images, with a lower value indicating more heterogeneity in the tumor. Most of the selected features in texture and wavelet classes reflected the pixel intensity distribution of ROI, which can be observed with the naked eye in tumors with poor prognoses (31). Apart from radiomic features, clinical features such as age and pretreatment plasma EBV DNA level were also frequently selected by the backward feature elimination process. This finding is in concordance with those of previous studies which have explored the correlation between age and disease prognosis in patients with NPC. It has been reported that younger patients had a better prognosis than older patients (36–39), and that plasma EBV DNA levels were linked to NPC early identification, prognostication, and treatment response (8–12).

Variability among observers remains a challenging issue in manual tissue segmentation, particularly for tumor delineation. In addition to intra-class correlation (ICC) analysis, which was performed to select consistent features for model construction, we also evaluated our model’s robustness. Specifically, we simulated the variation across observers based on tumor delineations by three radiation oncologists to generate synthetic radiomic feature values. The results showed that our models were robust with coefficient of variations less than 10% for both AUC and Harrell’s C-index.

There have been radiomics studies on NPC prognosis. MRI-based radiomics showed great promise for the prognosis of non-metastasis NPC patients (25, 27, 28, 40). Peng et al. (29) employed PET/CT-based radiomics as well as several clinical features to build a model to predict disease-free survival. They demonstrated that a radiomics model based on PET/CT outperformed a clinical model. In contrast to our findings, the clinical model produced better predictive performance than CT-based radiomic features. Furthermore, to the best of our knowledge, the only CT-based radiomics study on NPC was Zhu et al. (41), which used CT-based radiomic and clinical features including plasma EBV DNA level to predict local recurrences after IMRT in a cohort of 156 NPC patients from a hospital in China. Their results were also consistent with our findings even though we studied CT-based radiomics for OS, PFS, and DMFS prediction in non-Chinese patients.

Limitations of this study were that it was a retrospective study in a single center and that tumor delineation was based on CT images with MR images shown side by side. However, our findings serve as a proof of principle that CT-based radiomics is useful and we could identify important radiomic features that would be beneficial for future research. One of our future directions would be to construct a predictive model using radiomics from ROI drawn from CT images without support from MR images. Most importantly, the proposed predictive models should be further validated in a larger dataset from multiple centers or in patients who had induction chemotherapy followed by concurrent chemoradiation. The strength of this study includes a homogeneous group of patients who had been treated with IMRT concurrent with chemotherapy followed by adjuvant therapy. To the best of our knowledge, this is the first and largest study evaluating CT-based radiomics combined with clinical features on non-Chinese nasopharyngeal cancer patients to predict OS, PFS, and DMFS.

In conclusion, our study showed that CT-based Radiomic features, when used in conjunction with conventional clinical features, were able to improve the prognostic prediction performance for OS, PFS, and DMFS in NPC patients. This has a positive impact, especially for non-MRI institutes, on the screening of high-risk patients for aggressive therapeutic treatment strategies.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

The studies involving human participants were reviewed and approved by the Institutional Review Board, Faculty of Medicine, Chulalongkorn University, Bangkok, Thailand. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

SI: methodology, software, formal analysis, investigation, writing—original draft, visualization. YC: software, formal analysis. MV: software, formal analysis, investigation. SO: investigation. SK, AP, DK, and CC: data curation. SS: methodology, software, formal analysis, investigation. CL: conceptualization, methodology, resources, data curation. YR: conceptualization, methodology, formal analysis, investigation, resources, supervision, project administration, funding acquisition. All authors contributed to manuscript revision and read and approved the submitted version.

Funding

This project was funded by the National Research Council of Thailand (NRCT) (Grant: NRCT5-RSA63001-14). Sararas Intarak was supported by the 100th Anniversary Chulalongkorn University Fund for Doctoral Scholarship.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.