Multi-center evaluation of artificial intelligent imaging and clinical models for predicting neoadjuvant chemotherapy response in breast cancer

Neoadjuvant chemotherapy (NAC) plays an important role in the management of locally advanced breast cancer. It allows for downstaging of tumors, potentially allowing for breast conservation. NAC also allows for in-vivo testing of the tumors’ response to chemotherapy and provides important prognostic information. There are currently no clearly defined clinical models that incorporate imaging with clinical data to predict response to NAC. Thus, the aim of this work is to develop a predictive AI model based on routine CT imaging and clinical parameters to predict response to NAC. The CT scans of 324 patients with NAC from multiple centers in Singapore were used in this study. Four different radiomics models were built for predicting pathological complete response (pCR): first two were based on textural features extracted from peri-tumoral and tumoral regions, the third model based on novel space-resolved radiomics which extract feature maps using voxel-based radiomics and the fourth model based on deep learning (DL). Clinical parameters were included to build a final prognostic model. The best performing models were based on space-resolved and DL approaches. Space-resolved radiomics improves the clinical AUCs of pCR prediction from 0.743 (0.650 to 0.831) to 0.775 (0.685 to 0.860) and our DL model improved it from 0.743 (0.650 to 0.831) to 0.772 (0.685 to 0.853). The tumoral radiomics model performs the worst with no improvement of the AUC from the clinical model. The peri-tumoral combined model gives moderate performance with an AUC of 0.765 (0.671 to 0.855). Radiomics features extracted from diagnostic CT augment the predictive ability of pCR when combined with clinical features. The novel space-resolved radiomics and DL radiomics approaches outperformed conventional radiomics techniques.


Introduction
Neoadjuvant chemotherapy (NAC) is currently the standard of care for locally advanced breast cancer [1]. NAC allows for downstaging of bulky tumors, allowing previously inoperable tumors to become operable and allowing for breastconserving surgery (BCS) [2,3] where mastectomy was only possible. In addition, it potentially spares axillary clearance in biopsy-proven node-positive disease that resolves clinically [4]. Relatively high pathological complete response Tan Hong Qi and Ong Hiok Hian have contributed equally to this study.
The members of JBCR and Ai3 have been provided in Acknowledgements section.
Extended author information available on the last page of the article 1 3 (pCR) rates have been observed in HER2+ and triplenegative breast cancer (TNBC) [5]. Response to NAC also reflects tumor chemosensitivity and prognosis with complete primary and nodal pathological response associated with improved event-free as well as overall survival [6] and conversely, poor response or progression during NAC can result in delayed surgery, inoperability, greater risk of local progression and distant metastasis.
Aside from selection of breast cancer subtypes more sensitive to chemotherapy and/or targeted therapy, there are no good clinical models to predict NAC response. Use of functional imaging [7] and incorporation of biomarkers [8] have been studied to improve response predictions, but validated and accurate clinical models are still lacking.
Radiomics, or the extraction of quantitative radiomic features from medical imaging coupled with machine learning [9][10][11] or deep learning (DL) [12], may offer a solution. It has been used for other image-based classification problems such as diagnosing pulmonary tuberculosis from chest X-radiographs [13] and diabetic macular edema from optical coherence tomography images [14]. In predicting response to NAC in BC, promising results have been reported using intra-and peri-tumoral textural analysis of multi-parametric MRI (mpMRI) and CT images [15][16][17][18]. The authors thus aim to develop on AI model based on CT breast images to predict NAC response in BC.

Patients
The Centralised Institutional Review Board approved this retrospective study (IRB 2019/4219 Outcomes Research in Breast Cancer Care). 361 patients treated with NAC from January 2014 to December 2017 were included in this study. 37 patients were excluded due to missing clinical information or tumors that were difficult to define on CT scans, leaving 324 patients. Of these, 173 patients had their scans performed in National Cancer Centre Singapore (NCCS) and Singapore General Hospital (SGH), while the remaining 151 patients had their scans done in Kandang Kerbau Women and Children Hospital (KKWCH) and Changi General Hospital (CGH) in Singapore. The former and latter cohort was assigned as the training dataset and testing dataset, respectively. The binary endpoint of interest in this study is pCR which is defined as achieving ypT0N0 after NAC.
All patients were staged according to the American Joint Committee on Cancer 7th edition. Histological subtypes of cancers were approximated from the tumor estrogen receptor (ER), progesterone receptor (PR) or collectively as hormone receptors (HR) and HER2 receptor positivity: HR + HER2− (ER and/or PR+ and HER2−), HR + HER2+ (ER and/or PR+ and HER2+), HR − HER2+ (ER and PR−, and HER2+) and HR − HER2− (ER and PR−, and HER2−). ER and PR staining intensity for nuclear positivity of the cells were graded as 0, 1, 2, and 3 for none, mild, moderate, and strong intensity, respectively. HER2 positivity is reported as per American Pathologists guideline recommendations.

Analysis pipelines
Four different radiomics algorithms were applied on the training datasets. Two of the algorithms rely on conventional radiomics approach applied to different segmentations-tumoral (entire gross tumor volume, GTV) and peri-tumoral regions. This approach can only extract a single radiomics feature vector representative of the entire segmentation and the meaning of this vector can be nebulous when distinct heterogeneous regions exist in it. The third approach was developed by us to overcome this shortcoming. This novel technique makes use of voxel-based radiomics which extracts radiomics feature from small kernel in the segmentation. The outputs are feature vectors defined for every voxel in the segmentations. Further processing is carried out to condense the high-dimensional features into practical feature vector for modeling purpose. This method is known in the manuscript as space-resolved radiomics to distinguish it from the conventional radiomics. The last method used DL approaches on a fix-sized volume around the centroid of the GTV.
Feature reductions' techniques were applied to the radiomics feature vectors extracted from the first three algorithms before being input into the LASSO algorithm [19] for model building. The first feature reduction technique involved removing correlated features with Spearman correlation greater than 0.80. The second feature reduction technique only applied to the first two algorithms and it removes non-robust features according to a modified method by Zwanenburg et al. [20]. A 10-fold cross-validation was used to determine the optimal lambda parameter in the LASSO model. The DL radiomics approach directly learnt the model and output the predicted probability from the input image and no features preprocessing are required. These constitute the radiomics model in this manuscript. Clinical model was constructed using stepping logistic regression with statistically significant covariates under univariate analysis. Lastly, clinico-radiomics models were obtained using similar stepping logistic regression method on the similar set of clinical covariates and an additional variable known as radiomics score. The radiomics score is calculated by applying a logit function to the predicted probability from the radiomics approach.

Image acquisition, segmentations and radiomics feature extraction
The 3D volumetric images in the entire dataset were acquired and reconstructed with 11 different CT scanners. The images with contrast were captured with two different CT scanners located in the Centre's Radiotherapy Department. The first was the GE LightSpeed VCT and the second was Siemens SOMATOM Definition. The CT scanners from external hospitals were made up of primarily GE Revolution and Canon Aquilion model. The images were reconstructed with 90 to 140 kVp X-ray with variable slice thickness of 2.0 mm, 2.5 mm, 3.0 mm and 5.0 mm. The focal spot of the scanner ranges from 0.7 to 1.2 mm. The in-slice resolution was 512 by 512 for all images and the in-slice pixel widths are sub-millimeters.
The GTV was contoured by two experienced radiation oncologists (F. Y. Wong and W. L. Nei) for all the training and testing datasets without knowledge of the patient's pathologic outcome. These segmentations were known as the Intensity Mask under the IBSI naming convention.
Radiomics features were extracted using pyradiomics v3 [21]. It is comprised of shape, first-order and secondorder textural features (GLRLM, GLSZM, NGDTM and GLCM). The features are compliant with the Image Biomarker Standardization Initiative (IBSI) [22]. CT image was first interpolated with 1 mm isotropic voxel before feature extraction [23,24]. The radiomics features were extracted from the original CT image and the filtered CT image. LoG (Laplacian of Gaussian) filter with sigma 1.0, 2.0 and 3.0 mm and wavelet filters were applied on the image. This gives a total of 1130 radiomics features per CTV and GTV segmentation. A constant bin width of 10 HU was used for textural calculation.
We designed a procedure to select a subset of features from the original 1130 radiomics features which are robust to CT scanner variation and inter-rater variation in GTV contouring. The contouring variability is simulated by performing morphological dilatation and erosion operations of up to 2.0 mm on the segmentation (this generated four additional structures per patient). The CT scanner variability is simulated by adding an independent Gaussian, Poisson and Uniform noise to the CT image. The mean of the Gaussian, Poisson and Uniform noise distributions is 20 HU and the standard deviation of the Gaussian and Uniform distribution is 20 HU as well. These parameters were estimated from a phantom study on the two CT machines in NCCS. The process of adding the noise to the original CT image was repeated five times for each distribution, resulting in fifteen different CT image for each patient. Overall, the morphological operations and noise addition generated 19 additional set radiomics features per patient.
Details are shown in Fig. 3. Intraclass correlation coefficients (ICC) with a threshold of 0.600 are used to select the robust features.

Tumoral and peri-tumoral radiomics
Tumoral radiomics features are features extracted directly from the GTV contoured by the physician, and represent the textural, shape and intensity descriptors of the entire tumor core. The peri-tumoral regions are obtained by extracting an annulus ring segment from the Boolean subtraction of GTV with 1 mm dilatation and 1 mm erosion. These features describe the gray level characteristic of the peripheral of the tumor.

Space-resolved radiomics
Voxel-based radiomics are first used to extract features heat map from the GTV segmentation. Voxel-based extraction involves defining a kernel size (5 voxels isotropically in this work) and calculating the radiomics features (firstorder and textural) for each kernel centered on every voxel within the GTV. The choice of the kernel size can be interpreted as the degree of locality the texture calculation is sensitive to; a small kernel size could result in a noisy texture map, while a large kernel size will lose the local information and become increasingly smoother. A total of 88 first-order, GLCM, GLRLM, GLSZM and NGDTM feature maps are calculated for each GTV.
It is impractical to use all the feature maps directly for modeling; thus, the second step involves condensing all the information in the 88 feature maps into a fixed size feature vector for each GTV. The following procedure is applied consistently to all the 88 feature maps. The feature values for all the voxels in the entire training dataset are combined together and are used to generate the bin edges values by constructing a histogram with 20 equal-sized bins. Feature values beyond the 2 and 98 percentiles are excluded while calculating the bin edges of the histogram to prevent the influence of outliers. The bin edges of the histogram are then used to bin all the feature values of all the voxels in each GTV and the percentage number of voxels in the 20 bins are used as the feature vectors to represent the particular feature map.
Space-resolved radiomics has another advantage over conventional radiomics approach in that it can show which are the voxels within the segmentations that are driving the prediction model. The significant features in the model are first identified and voxels within the segmentation with feature value lying within the bin can be highlighted.

Deep learning radiomics
A fixed size box of dimensions 64 × 64 × 32 voxels centered at the centroid of the GTV is extracted and serves as the input image to the DL algorithm. Data augmentations consisting of a random translational shift of 1 to 5 voxels (uniform distribution) is applied to the box in all three directions, random 90° rotations about the z axis and 180° rotations about the x and y axes were applied to the box. MultiResUnet3D architecture [25] was selected for this work and the network was trained from scratch using Tensorflow framework. Adam optimizer was used with a learning rate of 1.0 × 10 −4 with exponential decay. Several modifications were made to the open source architecture. Firstly, He normal weight initialization was used and a L2 weight regularization of 1.0 × 10 −7 was applied to all Conv3D layer. Secondly, the last Conv3D was removed and replaced with a Flatten layer follow by 128 units fully connected layer with reLu activation function and a 1-unit fully connected layer with sigmoid activation unit. A 0.50 and 0.20 dropout were applied to the first and second fully connected layer, respectively. This architecture modification is necessary because the original architecture was meant for segmentation task with probability output per voxel while we are working on binary classification task with a single probability output. The 173 patients are divided into 120 training data and 53 validation data. The model with the highest validation AUC is applied to the testing data. The testing datasets are similar across all the different algorithms. Gradient-weighted class activation map (Grad-CAM) [26] is calculated at the encoded layer of the network to highlight the important regions in the image contributing to the pCR prediction.

Assessing performance of models
The performance of the four algorithms was evaluated on the test datasets. For each algorithm, three different models comprising of clinical, radiomics and clinico-radiomics are generated for comparison. Two performance metrics are used to compare between algorithms and determine if the radiomics model value-add to the clinical model. The first performance metric is the AUC which is calculated from the ROC curve and serves as a measure of the discriminatory ability of the algorithm. The second performance metric is the calibration curve which measures how closely the predicted probability matches the observed one, and Brier's score is calculated to quantify the degree of departure from logistic regression model. The observed probabilities are partitioned into ten groups for the calculation of calibration curve. A lower Brier's score indicates better agreement between the predicted and observed probability.

Statistical analysis
All statistical tests were performed using R statistical software (version 3.4.2; R Foundation for Statistical Computing). A two-sided P-value < 0.05 was considered significant.

Patient clinical characteristics and analysis pipelines
A total of 361 patients met the inclusion criteria and were included in this study. 37 Patients were excluded from the study, leaving 324. Of the 37 patients, 23 had missing clinical information, while 13 patients had tumors that were difficult to define on CT scans (such as diffuse tumors or multifocal tumors). The baseline characteristics of the patients in the training and testing datasets are shown in Table 1. The patient characteristics distributions were similar across the two datasets. The number of pCR events are 21.4% and 24.5% in the training and testing dataset, respectively. The schematic of the training and testing dataset split is shown in Fig. 1A. Univariate analysis was also performed, and the P-values are shown in the table. HER2, ER, PR and histological subtypes were statistically significant under the analysis. Information on the chemotherapy drug usage in this study cohort is shown in Supplementary Table 1. The chemotherapy usage was evenly distributed across the two datasets with frequent usage of cyclophosphamide, taxanes and anthracyclines. No precise chemotherapy regimen information was known for 7 (0.9%) patients. More than 90% of patients received at least third-generation chemotherapy, namely taxanes in addition to an anthracycline-based regimen. Platinum-based chemotherapy was given to 52.5% of TNBC patients. In patients with HER2+ breast cancer, 94.2% received trastuzumab. Of these, 24.5% also received pertuzumab.
The overview of the analysis pipeline is shown in Fig. 1B. Three of the four radiomics approaches require significant feature preprocessing (scaling and selection) before being input into the model building pipeline, while the DL approach directly output the predicted probability of pCR. The combined model is then compared primarily with the clinical model to determine if radiomics features augment the predictive power clinical model.

Image acquisition, segmentations and radiomics feature extraction
Siemen SOMATOM Definition and GE LightSpeed VCT make up 68.2 and 31.8% of the CT machines in the training dataset (from our center). Canon Aquilion and GE Revolution make up 35.7 and 19.9% of the CT machine in the testing dataset. These two represent the majority of the CT machines used in this study. The remaining CT machines in the testing dataset comprise Philips iCT, Philips Mx8000 Quad, Siemens SOMATOM Definition and Siemens Sensation Cardiac.
Examples of the segmentations of the tumoral and peritumoral regions are shown in Fig. 2. The purple-shaded regions in the images in the first column are the GTVs contoured by our radiation oncologists. The purple ring in the second column is the peri-tumoral regions extracted using morphological operation. A single feature heatmap obtained from voxel-based radiomics is shown in the third column for comparison to show that our method is inherently different from the conventional radiomics approach.

Tumoral and peri-tumoral radiomics model
The overview of the robust feature selection procedure is shown in Fig. 3. Noise addition to the CT image and morphological erosion and dilatation operations on the segmentations are shown in the figure. ICCs are calculated for 19 sets of radiomics features extracted from the GTVs (15 from random noise addition and 4 sets from morphological operations) and a robust feature is defined by having ICC value greater than 0.600. The plots of the right show the robust features as being yellow color coded, while the non-robust features are blue color coded. First-order, GLCM, GLRLM and GLDM robust features are shown in the figure as an example. In practice, the ICCs are calculated for all 1130 radiomics features. A total of 367 features are determined to be robust by this procedure. These robust features are used for modeling (after removing correlated features) in tumoral and peri-tumoral radiomics.

Space-resolved radiomics model
The comparison between the conventional radiomics approach and our novel space-resolved radiomics is illustrated in Fig. 4. The main difference lies in the radiomics feature extraction process; conventional radiomics extract only a single feature vector for each segmentation (shown as ⃗ f in the left blue box) whereas our approach relies on voxel-based radiomics where the feature vector, ⃗ f , is extract for each kernel. The output of voxel-based radiomics is   Fig. 4. The last step of the space-resolved radiomics pipeline involves binning the feature values in each feature heatmap into histogram; all the histograms from the 88 feature heatmaps are concatenated together to form a single vector before being input into the model. The significant radiomics features and its associated coefficients in the radiomics model are shown in Table 2. The suffix represents the bin number of the histogram, i.e., Original_firstorder_RootMeanSquared_B9 means the ninth bin in the orginal_firstorder_RootMeanSquared histogram. The final clinico-radiomics space-resolved radiomics model consists of ER, HER2 and radiomics score with coefficients − 1.70 (− 3.85 to 0.44), 1.13 (0.43 to 2.64) and 1.96 (0.75 to 2.62), respectively.
The highest coefficient features in the model are selected and the voxels represented by each feature is shown in Fig. 5. The left and right images show six different true pCR and true non-pCR, respectively. The top, middle and bottom rows show voxels with original_glrlm_ShortRunLowGrayLevelEmphasis lying within the second, eighth and ninth bin, original_gldm_ DependenceVariance lying within the second and seventh bins and original_firstorder_RootMeanSquared lying within the ninth bin, respectively. The voxels driving the prediction for each feature are shown in purple.

Deep learning radiomics model
An overview of the DL radiomics pipeline is shown in Fig. 6A. A fixed size bounding volume around the centroid of the tumor is used as an input into the DL architecture. The model was trained for 120 epochs and the highest validation AUC model was used as the final DL radiomics model. The class activation maps of a true-positive and true-negative example calculated in the axial slice of the tumor centroid are shown in Fig. 6B. The map shows that the tumor contributes significantly to the explanation of pCR in the DL model. The highest achieved validation AUC is 0.746 and it happens at the fifth epoch.

Assessing performance of models
The results of the ROC curve comparisons between the four algorithms and between clinical, radiomics and clinico-radiomics models are shown in Fig. 7 (for train data) and Fig. 8  (for test data). The blue, yellow and green lines represent the clinical, radiomics and clinico-radiomics model, respectively. A higher AUCs are obtained for all the combined models compared to the clinical models in the train data. The greatest increase is observed in the space-resolved model followed by the peri-tumoral model. Looking at the test data results in Fig. 8, only the peri-tumoral, space-resolved and DL radiomics approach are able to enhance the AUC of clinical model. Similar to the train data, the highest enhancement is achieved in space-resolved radiomics where the clinical model AUC increases from 0.743 (0.650 to 0.831) to 0.775 (0.685 to 0.860) by including the radiomics score. The DL model achieves similar magnitude enhancement of 0.772 (0.685 to 0.853). The ROC curves comparison of all the four clinico-radiomics models is shown in Fig. 9. Tumoral radiomics model clearly underperforms compared to the rest of the models.
ROC curves are also calculated for all the histological subtypes in the test data using the four radiomics models. The results are shown in Supplementary Table 2. The performances of the models are generally non-uniform across the subtypes. Most of the models are able to enhance the AUCs of the clinical models except for the tumoral radiomics model which results in a lower combined AUCs than clinical models in HR + HER2−, HR + /HER2+ and HR − HER2− subtypes. Space-resolved and DL radiomics models performed the best in the TNBC (HR − HER2−) subtype (with highest combined AUCs), while the tumoral radiomics model performed the best in HR − HER2 + subtype. The compilation of the calibration curve results of all the four models is shown in Fig. 10. The blue, yellow and green lines represent the clinical, radiomics and clinico-radiomics model, respectively. The Brier score of the individual model is shown in the parenthesis in the legend. Similar to the ROC curve results, only peri-tumoral, space-resolved and DL models show improvement in the Brier's score from pure clinical model alone.

Discussion
Our results suggest that the addition of radiomic features to clinical data can improve the prediction of pCR after NAC in BC. A key advantage of this approach is that it uses CT imaging which is commonly available and we have tested this in an external validation cohort. Apart from the tumoral radiomics model, the peri-tumoral, space-resolved and DL models show improvement on the clinical model. This suggests that heterogeneity within certain parts of the tumor, in particular, the peri-tumoral ring, may yield Fig. 4 Comparison of conventional and space-resolved radiomics pipelines. This figure compares the conventional and space-resolved radiomics pipeline. The left blue box shows that a single radiomics features is extracted from each segment, while the right red box shows that radiomics feature maps are extracted from the spaceresolved radiomics pipeline. A histogram is generated from each radiomics feature map to characterize the tumor heterogeneity more information about the tumor than analyzing the whole tumor alone. This is consistent with other studies [17,26]. Braman et al. [17] showed that a combination of intra and peri-tumoral radiomics was able to give a good prediction pCR in MRI images. Features from the 0-to 3-mm peritumoral region were significantly associated with the density of tumor-infiltrating lymphocytes. Similar conclusion on the importance of peri-tumoral region is found in a separate radiomics study conducted by Sun et al. [27]. The author found that CT radiomic signatures extracted from 4-mm peri-tumoral ring were correlated with CD8 + T-cell infiltration and response to immunotherapy.
Our novel spatial radiomics method improves on the peritumoral ring method which extracts only a single feature vector for each segmentation (shown as ⃗ f in the left blue box in Fig. 4). Our approach relies on voxel-based radiomics where the feature vector is extracted for each kernel (as shown in the red box of Fig. 4). The improved discriminatory ability could be attributed to the fact that the proposed method extracts more information from the lesion and better characterizes the tumoral heterogeneity. Furthermore, this method is capable of outputting an equivalence of "class activation map" [28] which shows the region within the segment driving the prediction (shown in Fig. 5). In fact, from Fig. 5, one can observe the difference in patterns and regions driving the prediction between pCR and non-pCR patients.
Currently, to the authors' knowledge, there are no DL models developed yet that can predict NAC response in BC. While traditional radiomics depend on thousands of user-defined features, DL radiomics is not constrained with  Comparison of ROC curves of the four different algorithms using clinical, radiomics and combined of the training data. A-D ROC curves of the clinical (in blue), radiomics (in yellow) and com-bined (in green) models for peri-tumoral, tumoral, space-resolved and deep learning algorithms, respectively. The 95% CI of the AUC are shown in parenthesis Fig. 8 Comparison of ROC curves of the four different algorithms using clinical, radiomics and combined of the test data. A-D ROC curves of the clinical (in blue), radiomics (in yellow) and combined (in green) models for peri-tumoral, tumoral, space-resolved and deep learning algorithms, respectively. The 95% CI of the AUC are shown in parenthesis this pre-determined set of features but uses a convolutional neural network to learn relevant radiomics features for its task. Our result shows the DL radiomics approach outperforms conventional radiomics approach (both peri-tumoral and tumoral), and interestingly, matched our space-resolved radiomics method in terms of AUC performance. Diamant et al. [29] has similarly applied this DL radiomics approach to predict outcomes in head and neck cancer and found that it outperformed the traditional radiomics prediction.
A common difficulty in the application of these radiomics model has been its generalizability. Different scanners and scanning parameters can give vastly different results and there is a potential of overfitting on the training data and generalizing to out-of-class distribution [30,31]. We mitigated the problem of overfitting by using two approaches. First, our model has been trained and tested on CT images retrieved from various machines and hospitals. Second, we included a robustness pipeline in our conventional radiomics approach to select our features that are stable to noise and contour perturbations.
An advantage of using CT imaging is that they are usually done as part of the staging workup for breast cancer patients (c ≥ T2 or cN+) planned for NAC. A limitation, however, is that as the paradigm shifts toward offering neoadjuvant therapy to earlier stages of breast cancer CT imaging may not be routinely performed in all centers.
Our models could potentially be easily integrated into current workflows. CT scanning protocols and parameters are also relatively more uniform across centers than compared to MRI. Furthermore, the CT voxel values (in Hounsfield units) are standardized across different scanners or protocol. This confers an additional robustness against scanner heterogeneity especially when there are many different brands and models of CT scanners. Following the definition of radiomics quality score (RQS) defined by Lambin et al. [32], our assessment of this work has a RQS of 18 out of 36 [33].
Despite the encouraging results from our study, the main limitation of our study is the retrospective nature of the study and the relatively small cohort size. Further validation of our model in a larger cohort and in prospective studies is being planned.

Conclusion
Peri-tumoral, space-resolved and DL CT radiomics features when combined with clinical features improve the prediction of pCR after NAC in BC. Conventional radiomics underperform our proposed space-resolved radiomics and DL radiomics approaches. The new proposed radiomics technique is a promising way to characterize tumor heterogeneity compared to conventional textural radiomics and could be used to develop better predictive and prognostic models.