Review on Multispectral Photoacoustic Analysis of Cancer: Thyroid and Breast

In recent decades, photoacoustic imaging has been used widely in biomedical research, providing molecular and functional information from biological tissues in vivo. In addition to being used for research in small animals, photoacoustic imaging has also been utilized for in vivo human studies, achieving a multispectral photoacoustic response in deep tissue. There have been several clinical trials for screening cancer patients by analyzing multispectral responses, which in turn provide metabolomic information about the underlying biological tissues. This review summarizes the methods and results of clinical photoacoustic trials available in the literature to date to classify cancerous tissues, specifically of the thyroid and breast. From the review, we can conclude that a great potential exists for photoacoustic imaging to be used as a complementary modality to improve diagnostic accuracy for suspicious tumors, thus significantly benefitting patients’ healthcare.


Introduction
In biomedical studies, the characterization of molecular and functional information about the underlying tissue can significantly improve our intuition in analyzing the morphology, treatment efficacy, and metabolites of target tissues. Among many biomedical imaging techniques, optical imaging methods have been widely applied for small animal studies due to their cost-efficiency, ease of implementation, non-ionizing radiation, and real-time imaging capability [1]. More importantly, optical imaging techniques can also provide molecular and functional information by tuning the wavelength of the light source. While advantageous, the strong optical diffusion of the pure optical imaging techniques in biological tissues leads to a reduced penetration depth, thus limiting clinical translation.
Photoacoustic imaging (PAI), one of the branches of optical imaging, provides the added advantage of increased imaging depth [2]. Compared to other optical imaging techniques, PAI inherits ultrasound imaging characteristics (USI), which increases its ability to visualize structural information in deep tissue. The signal generation in PAI is based on the photoacoustic (PA) effect, which is energy transduction from light to ultrasound (US) [3]. In brief, PA images can be achieved through the following procedure: (i) pulse laser illumination, (ii) light absorption by chromophores, (iii) momentary heat generation, (iv) acoustic wave (i.e., PA wave) generation through thermoelastic expansion, (v) signal detection by US transducer, and (vi) image generation. The resulting PA images, formed from the acoustic wave, include the optical absorption characteristics of the underlying biological tissue. Thus, PA images can provide molecular and functional information with a good ultrasound resolution in deep tissue [4,5]. In addition to the endogenous chromophores, contrast-enhanced PAI [6,7] has been studied by developing exogenous

Principles of Multispectral Photoacoustic Analysis
Similar to other optical imaging techniques, PA investigates the molecular composition of chromophores by applying spectral unmixing techniques to multispectral data [54]. From the images acquired with various excitation wavelengths, the spectral responses from different chromophores are delineated by mathematical calculations. The typical application of the spectral unmixing technique in PAI is in calculating hemoglobin oxygen saturation (sO 2 ) levels in blood vessels [55][56][57]. The sO 2 level in tumors indirectly represents metabolomic information about the tissue. Thus, PA-guided tumor assessment typically measures the sO 2 levels to differentiate benign from malignant tumors [58][59][60].
The typical unmixing method for multispectral PA data is linear unmixing, which can be simplified as a matrix multiplication as shown below: where P i (i = 1, 2, · · · , K) is the i-th pixel in images and λ j (j = 1, 2, · · · , N) is the j-th excitation wavelength, V P i , λ j is the measured value, µ HbO λ j and µ HbR λ j are optical absorption coefficients, and C HbO (P i ) and C HbR (P i ) are concentrations of oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR), respectively. The concentrations of hemoglobins are calculated by taking the pseudo-inverse as shown below.
Consequently, the concentration of the total hemoglobin (HbT) and the sO 2 level in each pixel are derived using the following equations.
In addition to the linear unmixing technique, the model-based unmixing [61] and the deep-learning technique [62] have also been used for multispectral PA analysis.

Multispectral Photoacoustic Analysis of Thyroid Gland
Thyroid cancer is one of the most common cancers, with an increasing global incidence rate in men and women [63][64][65]. The gold standard for diagnosing thyroid nodules is fineneedle aspiration biopsy (FNAB) [66]. The triaging for FNAB of the nodule is determined by the characteristics of nodules in USI [67][68][69]. Although the sensitivity of US-guided triaging is greater than 90%, the lack of functional metabolomics results in a low specificity of 20-50% [70]. The high false-positive rate leads to unnecessary FNAB, which results in the over-diagnosing of the tumor. Thus, clinical trials have been conducted to enhance the accuracy of triaging thyroid nodules using PAI due to its molecular and functional imaging capability.
Dogra et al. analyzed 88 resected tissues (13 malignant nodules, 30 benign nodules, 13 colloid accumulations, and 32 normal tissues) from 50 patients (11 malignant and 39 benign) [71]. Four different wavelengths (760, 850, 930, and 970 nm) were used for the spectral unmixing of HbO, HbR, lipid, and water components from multispectral PA data. Statistically significant differences were found in HbO and HbR between malignant and other types of tissues ( Figure 1a). In particular, HbR components were significantly different between malignant and normal tissue, with a p-value of 0.003 in the student t-test. The results showed the promising feasibility of PA-guided classification with a sensitivity of 69.2% and a specificity of 96.9%, but this study was limited to ex vivo environments only. Thus, for clinical translation, further in vivo validation is needed.
Dima et al. demonstrated the in vivo imaging capability of their PA and US system for the human thyroid [72]. They recruited two healthy volunteers to acquire PA images with a single excitation wavelength of 800 nm (Figure 1b). US Doppler images were also acquired in the same region to verify the blood vessel's position (yellow box in Figure 1b). By comparing the PA images with the US Doppler images, surrounding blood vessels extending from the isthmus and carotid artery to the anterior of the thyroid gland were identified. The results showed the feasibility of in vivo PAI using the arc array US transducer by confirming the matched positions of blood vessels (white arrows in Figure 1b). However, the spectral analysis of cancerous nodules was not available in this study. Yang et al. compared in vivo PA images between papillary thyroid cancer (PTC) patients and healthy volunteers (Figure 1c) [73]. Although they achieved PA responses from cancerous nodules, the number of patients (10 PTC and 3 normal) included in the study was insufficient for statistical analysis. In addition, multispectral analysis was also not available in this study because they used a single excitation wavelength of 1064 nm. transducer by confirming the matched positions of blood vessels (white arrows in Figure  1b). However, the spectral analysis of cancerous nodules was not available in this study. Yang et al. compared in vivo PA images between papillary thyroid cancer (PTC) patients and healthy volunteers (Figure 1c) [73]. Although they achieved PA responses from cancerous nodules, the number of patients (10 PTC and 3 normal) included in the study was insufficient for statistical analysis. In addition, multispectral analysis was also not available in this study because they used a single excitation wavelength of 1064 nm. Roll et al. presented multispectral PA analyses for differentiating tissue disorders in the thyroid gland [74]. The composition of HbO, HbR, fat, and water were spectrally unmixed from the in vivo PA images of the enrolled patients (6 Graves' disease, 3 malignant, 13 benign, and 8 healthy), obtained using eight excitation wavelengths (700, 730, 760, 800, 850, 900, 920, and 950 nm). The sO2 levels of the thyroid were also visualized and investigated ( Figure 2a). The contours of the thyroid glands were determined by the corresponding US images. Statistical analyses demonstrated significant differences between diseased and normal thyroid tissues (Figure 2b). This study was also unable to calculate the classification accuracy due to a small number of samples, thus limiting the applicability of the results. Roll et al. presented multispectral PA analyses for differentiating tissue disorders in the thyroid gland [74]. The composition of HbO, HbR, fat, and water were spectrally unmixed from the in vivo PA images of the enrolled patients (6 Graves' disease, 3 malignant, 13 benign, and 8 healthy), obtained using eight excitation wavelengths (700, 730, 760, 800, 850, 900, 920, and 950 nm). The sO 2 levels of the thyroid were also visualized and investigated ( Figure 2a). The contours of the thyroid glands were determined by the corresponding US images. Statistical analyses demonstrated significant differences between diseased and normal thyroid tissues (Figure 2b). This study was also unable to calculate the classification accuracy due to a small number of samples, thus limiting the applicability of the results. Recently, Kim et al. presented a multispectral PA analysis with a statistically sufficient number of samples (23 PTC and 29 benign), the largest number of patients in a clinical thyroid study to date [75]. They achieved multispectral PAI using five wavelengths (700, 756, 796, 866, and 900 nm). The corresponding US data were also acquired simultaneously for delineating the boundary of nodules. Similar to the previous studies, the sO2 levels in nodules were acquired through the spectral unmixing of HbO and HbR ( Figure  2c). Three parameters were quantified and used for training the decision function ( Figure  2d): (i) PA spectral gradient: the slope of a first-order polynomial fitted line to the average value of the top 50% of PA signals within the nodule boundary at each wavelength; (ii) relative sO2 level: the average value of the top 50% sO2 values within the nodule; (iii) skewed angle of sO2 distribution: the skewed angle of the Gaussian-fitted distribution for the top 50% of sO2 values within the nodule. With the values of the three parameters scattered in a 3D plane, a support vector machine was trained to determine the 3D decision boundary, which showed a good classification accuracy with a sensitivity of 78% and a specificity of 93%. The classification accuracy was further enhanced using a novel scoring method (ATAP score), which combined a conventional USI-based scoring method (i.e., US, ultrasound; sO 2 , hemoglobin oxygen saturation; HbR, deoxy-hemoglobin; HbT, total hemoglobin; PTC, papillary thyroid cancer; SVM, support vector machine; ATAP, the American Thyroid Association and the PA probability of PTC; Se, sensitivity; Sp, specificity; CA, carotid artery; TH, thyroid; ND, nodule; **, p ≤ 0.01; *, p ≤ 0.05. The images are reproduced with permission from [74,75]. Recently, Kim et al. presented a multispectral PA analysis with a statistically sufficient number of samples (23 PTC and 29 benign), the largest number of patients in a clinical thyroid study to date [75]. They achieved multispectral PAI using five wavelengths (700, 756, 796, 866, and 900 nm). The corresponding US data were also acquired simultaneously for delineating the boundary of nodules. Similar to the previous studies, the sO 2 levels in nodules were acquired through the spectral unmixing of HbO and HbR (Figure 2c). Three parameters were quantified and used for training the decision function ( Figure 2d): (i) PA spectral gradient: the slope of a first-order polynomial fitted line to the average value of the top 50% of PA signals within the nodule boundary at each wavelength; (ii) relative sO 2 level: the average value of the top 50% sO 2 values within the nodule; (iii) skewed angle of sO 2 distribution: the skewed angle of the Gaussian-fitted distribution for the top 50% of sO 2 values within the nodule. With the values of the three parameters scattered in a 3D plane, a support vector machine was trained to determine the 3D decision boundary, which showed a good classification accuracy with a sensitivity of 78% and a specificity of 93%. The classification accuracy was further enhanced using a novel scoring method (ATAP score), which combined a conventional USI-based scoring method (i.e., ATA guideline score) and the photoacoustic probability of malignancy. The novel scoring method improved the sensitivity to 83% and the specificity to 93%. Thus, the results showed a great potential for enhancing the triaging accuracy of thyroid nodules using a multiparametric analysis of multispectral PA data as a complementary method to the conventional triaging method.
While PA analyses of thyroid nodules have been conducted by various groups worldwide (Table 1), the validation of multispectral PA analysis is still at the initial stage of evaluation. Further studies are required to address the following issues for successful clinical translation. (i) Larger number of patients are needed to enhance the reliability of this technique. (ii) In addition to PTC, the classification of other types of thyroid cancers such as follicular, medullary, and anaplastic thyroid cancers would expand the application area. (iii) Quantitative analyses of PA responses in skin color are needed. (iv) System improvement with a deeper imaging depth, faster frame rate, and smaller size would enhance the image quality for multispectral analyses.

Multispectral Photoacoustic Analysis of Breast
Breast cancer is the most common cancer for women worldwide [76,77], with more than 110,000 new cases estimated to be diagnosed in the United States alone in 2022 [78]. Similar to other cancers, early diagnosis of breast cancer increases the survival rates of affected patients [79]. Diagnostic imaging modalities such as X-ray digital mammography, magnetic resonance imaging (MRI), and USI are currently used for screening breast cancers. While X-ray mammography is the most commonly used screening device and has been found to reduce mortality, its sensitivity and specificity are low, especially in dense breast patients and small tumors [80][81][82]. MRI has a high sensitivity in detecting cancer but its specificity is low, resulting in unnecessary biopsies. In addition, MRI is cost-prohibitive and is not available for routine examination [83]. The low-cost, commonly available USI can overcome the disadvantages of mammography, especially in dense breasts, but the lack of functional information degrades its diagnostic accuracy [84]. Thus, there have been efforts to use PAI to increase diagnostic accuracy in breast cancer by providing molecular and functional information with multispectral analysis. By combining this with structural information (shape of nodules) obtained by USI, the PAI showed a great potential to distinguish benign from malignant tumors in vivo.
The feasibility of PAI for in vivo human breast has been evaluated by several researchers [85][86][87]. Hemispherical array transducers have typically been used for the 3D visualization of the breast. Kruger et al. successfully visualized 3D volumes of breasts from four healthy volunteers (Figure 3a) [88]. The resulting image clearly showed the blood vessels, with a sufficient contrast-to-noise ratio (CNR) of~200 and an imaging depth of 30 mm. However, since they used a single-wavelength (756 nm) laser source, the results lacked functional information. Schoustra et al. acquired the 3D PA images of two healthy breasts using arc-shaped detector arrays and a laser with two wavelengths (755 and 1064 nm) [89]. Although they used multiple wavelengths, the study's goal was focused on penetration depth and optical absorption (Figure 3b).
Metabolites 2022, 12, x FOR PEER REVIEW 7 of 14 molecular and functional information with multispectral analysis. By combining this with structural information (shape of nodules) obtained by USI, the PAI showed a great potential to distinguish benign from malignant tumors in vivo.
The feasibility of PAI for in vivo human breast has been evaluated by several researchers [85][86][87]. Hemispherical array transducers have typically been used for the 3D visualization of the breast. Kruger et al. successfully visualized 3D volumes of breasts from four healthy volunteers (Figure 3a) [88]. The resulting image clearly showed the blood vessels, with a sufficient contrast-to-noise ratio (CNR) of ~200 and an imaging depth of 30 mm. However, since they used a single-wavelength (756 nm) laser source, the results lacked functional information. Schoustra et al. acquired the 3D PA images of two healthy breasts using arc-shaped detector arrays and a laser with two wavelengths (755 and 1064 nm) [89]. Although they used multiple wavelengths, the study's goal was focused on penetration depth and optical absorption (Figure 3b). Lin et al. acquired tomographic PA images of the breast using a full-ring array transducer with a single wavelength of 1064 nm [90]. They used four data acquisition modules to increase the imaging speed of the system significantly. A volumetric PA image of the entire breast was achieved within a single breath-hold (~15 s), thus minimizing motionrelated artifacts in the images. In seven patients with malignant tumors, they quantified the vessel density of the tumors and compared it with the normal tissue (Figure 3c). The results clearly showed the differences between malignant and normal tissues, but the small number of samples limited statistical verification. In addition to the volumetric geometry, linear array transducers were also used for breast imaging. Wang et al. demonstrated the feasibility of PAI using a linear array transducer on a healthy volunteer [87]. From the same group, Nyayapathi et al. improved their system to achieve mammogramlike PA images by using linear array transducers located both superior and inferior to the Lin et al. acquired tomographic PA images of the breast using a full-ring array transducer with a single wavelength of 1064 nm [90]. They used four data acquisition modules to increase the imaging speed of the system significantly. A volumetric PA image of the entire breast was achieved within a single breath-hold (~15 s), thus minimizing motionrelated artifacts in the images. In seven patients with malignant tumors, they quantified the vessel density of the tumors and compared it with the normal tissue (Figure 3c). The results clearly showed the differences between malignant and normal tissues, but the small number of samples limited statistical verification. In addition to the volumetric geometry, linear array transducers were also used for breast imaging. Wang et al. demonstrated the feasibility of PAI using a linear array transducer on a healthy volunteer [87]. From the same group, Nyayapathi et al. improved their system to achieve mammogram-like PA images by using linear array transducers located both superior and inferior to the breast [91]. They recruited 38 patients with malignant tumors, and acquired PA images of both diseased and healthy breasts with a wavelength of 1064 nm [92]. They defined four quantitative parameters which were calculated as the ratio of values in cancerous and normal tissues. (i) Vessel mean ratio: the ratio of the mean values of PA signals in the vessels. (ii) Contrast ratio: the ratio of the contrast that was calculated as the vessel values divided by the background signal. (iii) Standard deviation ratio of vessels: the ratio of the standard deviations of PA signals in the vessels. (iv) Standard deviation ratio of background: the ratio of the standard deviations of PA signals in the background. They used these to compare the PAI results of malignant and healthy breasts. The results showed statistical differences with p-values of~0.05, but the diagnostic accuracy was not evaluated. Thus, the imaging results showed the feasibility of visualizing the vascular structure in human breasts, but multispectral analysis to assess the functional information was still missing.
Toi et al. investigated multispectral PA analysis using two wavelengths (755 and 795 nm) [93]. They recruited 22 patients with malignant tumors, and successfully visualized tumor-related blood vessels using a hemispherical array transducer (Figure 4a). They defined a novel imaging parameter named the S-factor, which depends upon the sO 2 level and HbO composition in blood vessels. The comparison of vessels before and after chemotherapy showed the degradation of the S-factor, which indicates hypoxia. Although they could acquire high-quality images for a relatively large number of subjects, no statistical analysis was performed for differentiating malignant masses from benign or normal tissues. Diot et al. evaluated multispectral PAI with 28 wavelengths for the human breast in the range of 700-970 nm at 10 nm intervals [94]. Using the spectral unmixing of the multispectral data, they obtained the distributions of HbO, HbR, lipid, and water within the tissue. The corresponding total blood volume, which was related to the HbT level, was also calculated. They compared the findings between 10 patients with malignant tumors and 3 healthy volunteers. Compared to the healthy tissues, the malignant tumors showed significantly high values of total blood volume, indicating angiogenesis in the tumorous region. Statistical analysis was still unavailable due to the small number of patients. Thus, further statistical studies with a sufficient number of subjects are necessary to indicate the reliability of the results.
Neuschler et al. reported classification results using multispectral PA data on a large number of subjects (1079 benign and 678 malignant cases) [95]. They also presented a statistical analysis of their results. They used a laser with two wavelengths (757 and 1064 nm) to visualize HbO, HbR, and HbT distributions in human breasts. From the acquired images, the authors quantified five scores based on the corresponding features ( Figure 4b). (i) Vessel score: the combination of the number of individually resolvable vessels and their relative degree of deoxygenation. (ii) Blush score: the average volume of vessels that are too small to distinguish. (iii) Hemoglobin score: the amount of hemoglobin. (iv) Boundary zone score: vascular morphology and deoxygenation in the boundary of nodules. (v) Peripheral zone score: the number of radiating vessels in the peripheral region. They combined the PA scores with the BI-RADS (breast imaging reporting and data system) grades, the conventional image-based risk stratification method for breast cancers. In benign cases, PA scores downgraded 34.5% of high-grade BI-RADS masses, while only 6.0% were upgraded from low BI-RADS grades. In contrast, PA scores upgraded 30.6% of malignant masses from low BI-RADS grades, and downgraded 16.5% of high-graded malignant masses. These results showed that the PA scores could correct the BI-RADS grades to improve diagnostic accuracy. The classification result yielded a sensitivity of 98.6% and a specificity of 43.0%. This is the largest clinical study with patients recruited from multiple sites, showing the clinical feasibility of multispectral PA analysis. The subsequent study validated the PA-based BI-RADS correction with 209 patients [96]. For 47.9% of benign masses, the proposed PA analysis method successfully downgraded the BI-RADS scores to two or three. The original scores of those masses were 4a or 4b, which are classified as highly suspicious for cancer in BI-RADS. The result showed the great potential of the PA-guided scoring method to decrease unnecessary biopsies or surgical operations for breast cancer patients. fined a novel imaging parameter named the S-factor, which depends upon the sO2 level and HbO composition in blood vessels. The comparison of vessels before and after chemotherapy showed the degradation of the S-factor, which indicates hypoxia. Although they could acquire high-quality images for a relatively large number of subjects, no statistical analysis was performed for differentiating malignant masses from benign or normal tissues. Diot et al. evaluated multispectral PAI with 28 wavelengths for the human breast in the range of 700-970 nm at 10 nm intervals [94]. Using the spectral unmixing of the multispectral data, they obtained the distributions of HbO, HbR, lipid, and water within the tissue. The corresponding total blood volume, which was related to the HbT level, was also calculated. They compared the findings between 10 patients with malignant tumors and 3 healthy volunteers. Compared to the healthy tissues, the malignant tumors showed significantly high values of total blood volume, indicating angiogenesis in the tumorous region. Statistical analysis was still unavailable due to the small number of patients. Thus, further statistical studies with a sufficient number of subjects are necessary to indicate the reliability of the results.  The breast is one of the accessible organs which can be imaged using PAI. In addition to being superficial with protruding geometry, the homogeneous structure of the breast is optically transparent compared to other organs, lending itself as an ideal organ to be imaged using PAI. Thus, many clinical trials have been conducted to identify breasts cancer using multispectral PA images (Table 2). These studies have demonstrated the feasibility of PAI with a relatively large number of subjects. However, further validation is required for its successful translation to the clinic, as outlined here. (i) Improved specificity is needed for the reliability of this method. (ii) The study of additional excitation wavelengths and multiparametric analysis will allow the opportunity to define additional quantification methods, increasing the classification accuracy further. (iii) Imaging with contrast agents will assess the lymph node metastasis, which is one of the key factors for treatment planning [97][98][99].

Discussion and Outlook
Multispectral PAI is a promising method for investigating the metabolism of cancer cells and the molecular composition of biological tissue without injecting any contrast agents. In most cases, the quantification of sO 2 level has been used for detecting tumor hypoxia, which is one of the well-known characteristics of cancerous tissue [100]. Using the functional information of PAI, clinical trials have been conducted for visualizing, diagnosing, and analyzing tissue in various clinical applications including cancers [36,101], melanoma [52,53], and vascular diseases [102].
This review summarized multispectral PAI for analyzing thyroid and breast cancers in humans. The multispectral PA analysis results showed the feasibility of this technique for classifying cancerous nodules from healthy tissues, but some limitations remain to be overcome for successful clinical translation. (i) The reproducibility of initial results should be verified and validated with larger-scale clinical studies. (ii) In addition to a simple cancer classification, various types of diseased nodules should be differentiable. (iii) A method to improve both sensitivity and specificity should be developed and verified. To solve the issues above, clinical PAI systems should be also improved with the following points: (i) real-time quantification for multispectral PA features, (ii) enhanced image quality, (iii) improved imaging depth, (iv) compact systems for better mobility, and (v) a userfriendly interface.
Although PAI is still in the initial stage of clinical translation, many studies have already established, with statistical analyses, the feasibility of distinguishing cancerous masses from benign tissues. With additional verification using a larger number of patients from multiple sites, a PAI-based modality could become an excellent method for screening cancer patients. In the near future, with sufficient data collected, machine learning techniques could be used to further increase classification accuracy.