Exploring the diagnostic value of photoacoustic imaging for breast cancer: the identification of regional photoacoustic signal differences of breast tumors

: We examined 14 benign and 26 malignant breast nodules by a handheld dual-modal PA/US imaging system and analyzed the data using the quantitative and semi-quantitative method. The PA signal spatial density and PA scores of different regions of the benign and malignant nodules were compared, and the diagnostic performances of two diagnostic methods based on PA parameters were evaluated. For both quantitative and semi-quantitative results, significant differences in the distributions of PA signals in different regions of benign and malignant breast lesions were identified. The PA parameters showed good performance in diagnosing breast cancer, indicating the potential of PAI in clinical utilization. The PA scores in the tumor region and central region of benign tumors were both significantly higher than that in malignant tumors. No statistically significant difference was determined between the S S and S P between benign and malignant tumors. These results were similar to differences of the PA signal spatial density between benign and malignant lesions. The mean S S − T , S S − C and S P − C of benign tumors was significantly lower than that of malignant tumors (S S − T : 0.62 ± 0.64 VS 1.25 ± 0.63, P = 0.005; S S − C : 1.51 ± 0.72 VS 2.14 ± 0.47, P = 0.02; S P − C : 0.40 ± 0.34 VS 0.83 ± 0.39, P = 0.001).


Introduction
As the most common malignancy among women globally, breast cancer caused 2.089 million new incidences and 627,000 deaths worldwide in 2018, and great importance has been attached to the early diagnosis and treatment of the disease to improve the prognosis of the patients [1]. At present, mammography and ultrasound are the routine imaging examinations for detecting breast cancers, both of which distinguish benign and malignant breast nodules based on morphological features [2,3]. Despite its high accuracy and wide acceptation in many countries, mammography is less efficient in identifying nodules of dense breast tissue. Besides, it is also radioactive and can cause uncomfortableness [4]. Breast ultrasound has gained increasing recognition by clinicians recently for its excellent performance in detecting malignant lesions, of which the use is still hindered by the relatively low specificity and operator dependence [5]. Additionally, there exists an overlap in morphological manifestations of benign and malignant breast nodules, which might be accounted for by the heterogenicity of breast tumors [6,7]. As a result, the diagnostic accuracy of mammography and breast ultrasound may be impaired only with the reference of morphology shown up on the images.
Breast Imaging Report and Data System (BI-RADS), enacted by the American College of Radiology and deemed the standard protocol for breast imaging diagnosis, proposed identifying morphological characteristics of lesions for making diagnosis [8]. Meanwhile, the role of blood flows is less emphasized by the system. However, the diagnostic role of blood vessels of breast two methods were integrated into two diagnostic models, and their diagnostic performances in differentiating benign and malignant breast lesions were investigated.

Methods and materials
The study (No. JS-1424) was designed as a prospective study and was approved by the Institutional Review Board of Peking Union Medical College Hospital. Informed consents were received from all the patients who underwent an examination.

Patient recruitment
From November 2017 to January 2018, the female patients with focal breast lesions were consecutively recruited from the breast surgery department of Peking Union Medical College Hospital. The breast lesions were detected by ultrasound before hospitalization, and the patients received ultrasound and dual-modal PA/US imaging examination after recruitment. The inclusion criteria of the breast lesions were: 1) the largest diameter of the lesion was smaller than 3 cm; 2) the upper boundary of the lesion was at least 5 mm from the skin surface; 3) the lower boundary of the lesion was within 3 cm from the skin surface due to the limited penetration depth of photons in tissue with high scattering coefficient.
After imaging examinations, the patients received ultrasound-guided core needle biopsy and/or surgery within one week. And all the lesions were then confirmed pathologically consequently.

Imaging system
The PA/US dual mode imaging was performed using a dual-modal PA/US imaging system with a handheld L9-3 linear probe (Resona 7, Mindray Bio-Medical Electronics Co., Ltd., Shenzhen, China). The horizontal and vertical resolution of the ultrasound and photoacoustic systems are less than 1mm. The probe had 192 elements, with a size of 0.2 mm for each element. The central frequency was 5.8 MHz. An OPO tunable laser (Spitlight 600-OPO, Innolas laser GmbH), which generated 700-850 nm laser pulses at 10 Hz, was also utilized. In this study, the wavelength was set at 830nm, where the OPO laser maintained at a high energy level, and the blood has a relatively strong optical absorption coefficient. The optical fluence on the skin for each patient was around 4 mJ/cm 2 , with fluctuations less than 5%. The PA signals were normalized by laser energies recorded in real time to account for the fluctuation of laser energy. In our study, for all PA images, we applied a threshold of 35% of the maximum pixel value. PA signals below the threshold were not shown in PA images and were not used for calculation.

Imaging procedure
The imaging examination was performed by a senior radiologist with 20-year experience of breast ultrasound and 3-year experience of PAI. Firstly, conventional ultrasound was performed for each patient to identify the anatomical structure of each lesion. Subsequently, a video was recorded for each tumor where the whole tumor as well as the peri-tumor region were scanned along the short axis of the tumor with dual-modal PA/US imaging. The 2D PA/US imaging result was then extracted from each frame of the video.

Quantitative analysis of PA signals and semiquantitative scoring of PA signals
For each nodule, the images containing the entire tumor and its surrounding 5mm tissue were extracted for analysis, and one frame in every 1mm was selected for analysis. Then we divided the tumor region and peri-tumor region into four parts based on their positions, listed as following.
Region S: the region surrounding the tumor but within 5 mm from the tumor boundary. Region T: the whole tumor region.
Region P: the periphery region inside the tumor and within length L from the tumor boundary, L is 1/8 length of the short axis of the tumor.
Region C: the central region inside the tumor but not belonging to Region P. Figure 1 demonstrates the four regions defined above in a straight way. These regions were identified and delineated based on the shape of each tumor, together with scanning time and the number of frames in each scanning. The region division was performed by two radiologists with five-year experience in breast ultrasound and one-year experience in PAI. One of the radiologists firstly performed the region division of the tumor. Then the other one double checked the results. When the boundary of the lesion was unclear and difficult to discern, the two radiologists would recheck the region divisions and discuss with each other to reach a consensus. The tumor is divided into the surrounding region (Region S) and the tumor region (Region T), and the tumor region is further divided into the central region (Region C) and the peripheral region (Region P). Three types of slices, including Slice 1, Slice 2 and Slice 3 are extracted from the full-length imaging data of a tumor. b. Section of Slice 1: the surrounding region. c. Section of Slice 2: the surrounding region + the peripheral region. d. Section of Slice 3: the surrounding region + the peripheral region+ the central region.

Quantitative analysis of PA signals
After sectioning regions with Image J [28], for each case, we calculated PA signal spatial density of each region averaged by slice that includes the correspondent region using MATLAB (Mathworks, Inc., USA). PA signal spatial density was defined as the quotient of the number of pixels with PA signals in a certain region and the total number of pixels in that region. In this study, we referred to the averaged PA signal spatial density of Region S, Region T, Region P, and Region C as P S , P T , P P . P C respectively. Besides, we also analyzed the relative difference of PA signal spatial density between Region S and Region C, Region P and Region C, and Region S and Region T, defined as P S−C = (P S -P C )/P C , P P−C = (P P -P C )/P c , and P S−T = (P S -P T )/P T , respectively.

Semiquantitative scoring of PA signals
In addition to the calculation of PA signal spatial density, we also performed semi-quantitative analysis by scoring each PA/US imaging result. The PA signals of each region in each image were scored according to Adler's method for Doppler scoring [29], as absent (Score 0), minimal (Score 1), moderate (Score 2), or marked (Score 3), with detailed criteria shown as follows: Score 0: no PA signals detected; Score 1: one or two bars of PA signals with a diameter less than 0.1 cm detected; Score 2: three to four bars of PA signals or a large bar of PA signals penetrating the lesion detected; Score 3: five or more bars of PA signals detected. PA score of each region in each case was averaged by slice including the correspondent region. For each case, the PA scores of Region S, Region T, Region P, Region C averaged by slice that includes the correspondent region were referred to as S S , S T , S P , and S C , respectively in the following analysis. We also analyzed the difference of PA scores between the Region S and Region C (S S−C , S S−C = S S -S C ), the difference of PA scores between Region P and Region C (S P−C , S P−C = S P -S C ), and difference of PA scores between Region S and Region T (S S−T , S S−T = S S -S T ). The PA semi-quantitative scoring for the regions was also performed by the two radiologists, who were blinded to the clinical information of the patients. If there was a discrepancy in scoring, the two radiologists would re-read the images and obtain the final scores for further statistical analysis.
The streamlined Standard Operating Procedure (SOP) for PA image acquisition and data processing was demonstrated in Fig. 2.

Statistical analysis
The software package SPSS for Windows 22.0 (SPSS Inc, Chicago, IL), R (http://www.R-project. org) and EmpowerStats software (X&Y Solutions), were used for the statistical data analysis. Continuous data were expressed as the mean ± standard deviation. Student's t-test was applied to evaluate the differences of signal intensity and scores between benign and malignant tumor regions. A P value of <0.05 was considered statistically significant.
In order to investigate the diagnostic performance of the two methods of quantitative analysis and semi-quantitative scoring of PA signals, we used a backward stepwise logistic regression method to construct two diagnostic models that combined the variables in the quantitative analysis (P S−C , P P−C , and P S−T ) and semi-quantitative methods (S S−C , S P−C , and S S−T ), respectively, based on the Akaike information criterion (AIC). And the bootstrap resampling method of random drawing of 500 samples was applied to perform the internal validation of the models. The combination of the variables of each method was presented in the form of an equation and a nomogram. The diagnostic performances of the two diagnostic models were evaluated by calculating a series of statistical indexes for diagnostic test, including sensitivity, specificity, positive likelihood ratio (PLR), negative likelihood ratio (NLR), predictive value (PPV), negative predictive value (NPV). The receiver operating curves (ROC) for the two diagnostic models were also delineated, and the area under the curve (AUC) values were calculated accordingly. We used 2 × 2 contingency tables, chi-square test for comparing sensitivity and positivity, a generalized estimating equation for comparing PPV and NPV, and the method proposed by DeLong et al. for comparing AUC values. Decision curve analysis (DCA) was also utilized to compare the diagnostic performance of the two models, which was designed to show the net benefit at a variety of threshold levels and weigh the harm of false-negative and false-positive diagnosis. The DCA was conducted according to the methods proposed by Vickers [30].

Results
From November 2017 to January 2018, a total of 40 solid breast tumors in 39 female patients whose ages ranged from 22 to 86 years (mean age, 45.6 years±14.0) were recruited from the breast surgery department of Peking Union Medical College Hospital. There were 14 benign and 26 malignant lesions, the benign lesions included 10 fibroadenomas, 2 sclerosing adenosis, one intraductal papilloma and one benign phyllodes tumor. The malignant lesions included 23 invasive ductal carcinomas and 3 ductal carcinomas in situ. (Table 1). The average depth (the distance from skin surface to the lower boundary of tumor) of the benign and malignant tumors were 1.46 ± 0.28 cm and 1.75 ± 0.38 cm, respectively. The average size (the maximum diameter) of the benign and malignant tumors were 1.88 ± 0.88 cm and 1.88 ± 0.68 cm, respectively. No statistically significant difference was found between either tumor depth or tumor size of the benign and malignant tumors.

Comparisons of the PA parameters of the benign and malignant lesions
The kappa value of the two radiologists was 0.825 (95% CI: 0.689-0.946, P < 0.001) in PA scoring, presenting almost perfect inter-observer agreement.  Fig. 3 show the statistical results of PA signal spatial density and the PA scores of the benign group and the malignant group. P T of the benign group was significantly higher than that of the malignant group (0.07 ± 0.04 VS 0.04 ± 0.03, P = 0.015). P C of the benign group was also significantly higher than that of the malignant group (0.06 ± 0.03 vs. 0.03 ± 0.02, P = 0.003). No significant difference of either P S or P P was found between the benign group and the malignant group, although both P S and P P of the benign group was slightly higher than those of the malignant group. From statistical results of P S−C and P P−C , we observed that for both the benign group and the malignant group, the PA signal spatial density in either Region S or Region P was larger than the PA signal spatial density in Region C. Interestingly, we also found that P S−C of the benign group was significantly lower than that of the malignant group (1.05 ± 1.53 VS 4.01 ± 4.23, P = 0.016). Similarly, P P−C of the benign group was significantly lower than that of the malignant group (0.46 ± 0.80 VS 1.45 ± 1.58, P = 0.036). The significant difference of P P−C as well as P S−C between the benign group and the malignant group indicates that the heterogeneity of lesion was more prominent in the malignant group than in the benign group. Similar as PA signal spatial density, S T of benign tumors significantly higher than that of the malignant ones (2.22 ± 0.68 VS 1.71 ± 0.66, P = 0.026). The mean S C of benign tumors was also significantly higher than that of the malignant ones (1.33 ± 0.77 VS 0.82 ± 0.51, P = 0.016). The PA scores in the tumor region and central region of benign tumors were both significantly higher than that in malignant tumors. No statistically significant difference was determined between the S S and S P between benign and malignant tumors. These results were similar to differences of the PA signal spatial density between benign and malignant lesions. The mean S S−T , S S−C and S P−C of benign tumors was significantly lower than that of malignant tumors (S S−T : 0.62 ± 0.64 VS 1.25 ± 0.63, P = 0.005; S S−C : 1.51 ± 0.72 VS 2.14 ± 0.47, P = 0.02; S P−C : 0.40 ± 0.34 VS 0.83 ± 0.39, P = 0.001).

Fig. 3. Boxplots of PA signal spatial density and PA scores in different regions of benign and malignant tumors.
Differences in PA signal spatial density and PA scores among different regions. P S−T (Relative difference of signal spatial density between Region S and Region T), P S−C (Relative difference of signal spatial density between Region S and Region C), P P−C (Relative difference of signal spatial density between Region P and Region C), S S−T (Difference of PA scores between Region S and Region T), S S−C (Difference of PA scores between Region S and Region C), S P−C (Difference of PA scores between Region P and Region C). The orange bars behind the boxplots represent the means of the parameters.

Diagnostic performance of the two diagnostic models
Using multivariate logistic regression analysis, two diagnostic models constructed from the combination of the variables of the quantitative analysis of PA signal spatial density and semiquantitative PA scoring were generated. The equations of the two models are listed as following: Logit(Y)=−1.52691+2.15350* P S−T −0.91390* P S−C +2.04078* P P−C ; logit(Y)=−4.37267+2.40903* S S−T −0.34354* S S−C +5.24288* S P−C . The nomograms that incorporated the variables are illustrated in Fig. 4. The comparisons of statistical indexes of diagnostic test between the two models are listed in Table 3 (Model 1 for quantitative PA analysis, Model 2 for PA semi-quantitative scoring), including sensitivity, specificity, PPV, NPV, and AUC, and the ROCs for the two models are presented in Fig. 5. The two models based on PA analysis demonstrated an excellent diagnostic performance with a best threshold of 0.567 and 1.16, respectively, between which no significant difference was detected in those statistical parameters (Sensitivity  0.865, P=0.635). The results showed that both of the models could be utilized in classifying benign and malignant breast lesions with similar diagnostic efficiency. Figure 6 illustrates the decision curves of the two models. The two curves overlapped between the threshold probabilities of 0.2-0.6, indicating both of the models were useful in differentiating malignancies with similar diagnostic capacity. Two typical cases using the diagnostic models to obtain imaging diagnosis is shown in Fig. 7.

Discussion
Although deemed as an established indicator of breast malignancy, the angiogenesis of the tumor has been proved to be highly involved with the occurrence, invasion, and metastasis of breast cancer [10]. High-quality vascular imaging can provide useful information for an accurate evaluation of the tumor and help predict treatment efficacy and prognosis, which still cannot be achieved by conventional breast imaging methods [14]. In this study, a handheld dual-modal PA/US imaging system was utilized to examine breast nodules of the 40 patients, which clearly displayed vascular distribution of the breast nodules. The diagnostic models based on quantitative and semi-quantitative PA analysis of the volume data presented good performance in differentiating benign and malignant lesions, indicating its potential clinical value for diagnosing breast cancers. In this study, volume data containing continuous cross-sectional images of breast lesions and the surrounding tissues were obtained, and the PA signals within and around the lesions were extracted for analysis. The differences in PA signal spatial density between the peripheral regions and the central regions, as well as the surrounding regions and the central regions of the malignant tumors are larger than that of the benign tumors, which reflect the heterogeneous distribution of blood vessels of the malignant tumor. On the basis of the above findings about the vascularization of breast tumors, we proposed the indicators representing the differences in PA signal spatial density of different tumor regions to illustrate the vessel distribution of the tumor, including the differences among surrounding, peripheral and central regions of the tumors (P S , P T , P P . P C , P S−C , P P−C , P S−T ; R S , R T , R P , R C , R S−C , R P−C , R S−T ). Subsequently, those quantitative and semi-quantitative indicators were integrated in two diagnostic models with the form of nomograms, both of which presented high diagnostic efficiency in detecting breast cancer with an AUC of 0.824 and 0.865, respectively, with no significant difference. Hence, the model based on semiquantitative PA scoring method also had favorable diagnostic capacity, not inferior to the quantitative analysis. With more clinical trials of newly developed PA systems launched out, the complex procedures of signal quantification and data analysis after imaging examination has posed a challenge to researchers, which might hinder the clinical promotion of PAI. Therefore, we put forward this semi-quantitative PA scoring with the reference of a widely used scoring system of Doppler imaging, to provide an easy and feasible way for the clinicians to interpret the PA images. In comparison with the quantitative analysis method, which requires further of image processing and calculation, the semi-quantitative PA scoring method is more time-saving and accessible. Also, the inter-observer agreement of the region division and PA scoring was proved to be very good in this study, thus validating the repeatability of the method. In terms that the two methods can achieve similar diagnostic efficiency in differentiating breast tumors, the semi-quantitative PA scoring method can be an alternative to quantitative PA analysis, which may also be further incorporated into the diagnostic strategies for breast cancer.
The results revealed that the signal spatial density within the benign tumors was higher than the malignancies, which was inconsistent with some of the previous studies that suggested a more abundant distribution of micro-vessels and higher blood perfusion of malignant breast lesions than that of benign ones [31]. This discrepancy could be explained by the constructions of the pathological types of the included lesions. For the nine fibroadenomas of the 14 benign cases, abundant PA signals were observed both in the central and surrounding regions with no significant difference of signal distribution. Previous studies have also suggested that rich blood flow signals might be detected in fibroadenoma, making it difficult to distinguish from malignant tumors and false-positive diagnosis [6,32]. Among the 26 breast malignant tumors in this study, 21 cases of the luminal type were identified with fewer micro-vessels within the lesions. This finding was consistent with the study of Dogan et al., who analyzed the relationship between the molecular subtypes and PA signals, and found that the internal PA signals of luminal types were sparser than the outside signals [33]. The results also showed that malignant tumors were characterized by rich PA signals in the peripheral areas tumor and sparse signals in the central areas. More signals in the peripheral areas were also detected in the benign ones, but the difference of signal distribution in malignant tumors between the central and peripheral areas is larger than that of the benign tumors. This finding is similar to several previous studies about breast PAI. Oraevsky et al. suggested that PAI's feature representing the abundant blood vessels in the peripheral regions had a positive predictive value of > 90% for diagnosing breast cancer [26]. Doit et al. also observed rich peripheral vascularization and a low level of intra-tumoral vessel distribution in the malignant tumors [27]. This might be accounted for by the fast growth and insufficient blood supply to the central area of the malignant tumor, resulting in a low level of hemoglobin and oxygenation, thus presenting low intensity of PA signals in the central areas of the tumors. On the other hand, the high intensity of the peripheral areas' PA signals represents the abundant micro-vessels and the highly invasiveness of the local tissues.
There exist several limitations in this study. Firstly, apart from delineating the distribution of blood vessels of the tumors, the system could also measure the oxygenation level at two wavelengths for more functional information, which needs to be further explored in the future study. Secondly, although high level of inter-observer agreement in PA scoring of the two radiologists was shown in this preliminary study, the stability of the PA scoring method should be validated with more participants, which requires further researches. Also, the sample size of this preliminary study was relatively small. To reduce the impact of small sample size as much as possible, we used nonparametric test in statistics. More over, in this study, we mainly focused on exploring the feasibility of using spatial characteristics of PA signals for diagnosis of breast cancer. To further validate the efficiency of the approach we proposed, the sample size should be enlarged and breast nodules with a variety of pathological and molecular subtypes should be recruited in future studies. At last, the diagnostic models in this study only incorporate PA indicators, and no conventional imaging parameters, such as the ultrasonic features of breast tumors, which are essential in the diagnosis of breast cancer, was included. Therefore, in order to construct a comprehensive diagnostic model for breast cancer, more information should be taken into consideration in further studies.

Conclusion
In this study, a dual-modal PA/US imaging system was used to examine the patients with breast lesions, and a PA scoring system based on the difference of signal intensity between the peripheral and central areas of the lesions was proposed. The difference in the distribution of PA signals in different regions of benign and malignant breast lesions was identified, which showed good performance in diagnosing breast cancer. The handheld dual-modal PA/US system possesses clinical potentials in diagnosing breast cancer by providing additional information of tumoral blood vessel distributions.