System-level optimization in spectroscopic photoacoustic imaging of prostate cancer.

This study presents a system-level optimization of spectroscopic photoacoustic (PA) imaging for prostate cancer (PCa) detection in three folds. First, we present a spectral unmixing model to segregate spectral system error (SSE). We constructed two noise models (NMs) for the laser spectrotemporal fluctuation and the ultrasound system noise. We used these NMs in linear spectral unmixing to denoise and to achieve high temporal resolution. Second, we employed a simulation-aided wavelength optimization to select the most effective subset of wavelengths. NMs again were considered so that selected wavelengths were not only robust to the collinearity of optical absorbance, but also to noise. Third, we quantified the effect of frame averaging on improving spectral unmixing accuracy through theoretical analysis and numerical validation. To validate the whole framework, we performed comprehensive studies in simulation and an in vivo experiment which evaluated prostate-specific membrane antigen (PSMA) expression in PCa on a mice model. Both simulation analysis and in vivo studies confirmed that the proposed framework significantly enhances image signal-to-noise ratio (SNR) and spectral unmixing accuracy. It enabled more sensitive and faster PCa detection. Moreover, the proposed framework can be generalized to other spectroscopic PA imaging studies for noise reduction, wavelength optimization, and higher temporal resolution.

Introduction

Prostate-specific membrane antigen (PSMA) is a type-II integral membrane protein. PSMA has been widely studied for molecular and functional imaging of prostate cancer (PCa) due to its higher expression on the cell surface of aggressive PCa than in healthy prostate tissue and indolent malignancies [1], [2], [3], [4], [5]. PSMA expression correlates well with staging and aggressiveness of PCa. Thus, it has been leveraged as an important molecular target for PCa management. Among several imaging modalities for PSMA detection, radionuclide imaging, e.g., positron emission tomography and single-photon emission computed tomography, enables whole-body imaging to detect primary PCa and metastatic lesions using affinity agents that link radionuclides to the PSMA [6], [7], [8]. Nevertheless, radionuclide imaging exposes patients to radiation and necessitates expensive infrastructure, hindering it from frequent surveillance of PCa patients. Another modality, optical imaging, has also brought opportunities for precise PSMA detection by employing contrast agents that binds to PSMA-positive (PSMA+) PCa [9], [10]. However, such agents only provide superficial imaging depth. Thus, optical imaging of PSMA is strictly confined to providing intraoperative guidance rather than being comprehensively useful in prostate cancer care. Hence, there is an urgent need for noninvasive, high-contrast imaging of PSMA expression in deep prostate tissues.

Photoacoustic (PA) imaging is an emerging molecular imaging modality that caters to the above requirements. It provides rich optical contrast at sub-millimeter spatial resolution and acoustic imaging depth over several centimeters [11], [12]. In PA signal formation, electromagnetic energy absorbed by the biological tissue converts to thermal energy, inducing a localized pressure increase due to the thermal-elastic expansion of the tissues. The pressure then propagates as acoustic waves and is recorded by acoustic detectors to form PA signals. Because most biological tissues possess unique optical absorption spectra, they can be unmixed by spectroscopic PA imaging presented in several novel clinical and scientific applications [13], [14], [15], [16], [17], [18], [19]. We recently validated the spectroscopic PA imaging of PSMA expression in vivo for the first time, using a PSMA-targeted agent [20], where spectral unmixing with ordinary least squares (OLS) discerned absorptive contrast of the agent specifically bounded to PSMA expressing PCa. However, there underline opportunities for significant improvements in the system level.

Here we focus on several challenges in our system-level optimization of the spectroscopic PA imaging frameworks. (1) The OLS method, the most common approach for spectral unmixing, suffers from many system noise sources that may annul the essential hypothesis of the flawless design matrix and zero-mean measurement noise. For example, normalizing the fluctuating pulse-to-pulse laser energy makes the system electronic noise dependent to the laser spectrum, therefore engenders spectral system error (SSE) and deteriorates spectral unmixing accuracy. (2) The redundant amount of spectral data ranging from 700 nm to 900 nm in 10 nm intervals in our previous study was still ineffective to show false contrast before the injection of the PSMA-targeted agent. Also, it suffers from low temporal resolution by using Nd:YAG pulsed laser equipped with optical parametric oscillator (OPO) commonly applied in the PA imaging research operating only at 10–20 Hz for sequential spectral tuning. Several wavelength selection algorithms were developed to reduce the computational complexity [21], [22], [23], but spectrotemporal instability of the tunable laser system were not considered, where the varying SSE at different wavelengths may affect spectral unmixing accuracy. (3) Frame averaging is a conventional method to improve the signal-to-noise ratio (SNR) in PA imaging, but another layer of analysis is still needed to understand how such SNR enhancement will improve the accuracy of spectral unmixing. Herein, we present a simulation-aided system-level optimization of the spectroscopic PA imaging of PSMA expression to achieve robust spectral unmixing accuracy and enhanced temporal resolution.

Methods

Conventional spectroscopic PA unmixing

The initial pressure generated by PA effect can be equated using the following equation [24].where is the PA efficiency based on the Grüneisen parameter, is the optical absorbance (m−1), and is the optical fluence (J/m2). In this work, the insignificant spatial variance of is omitted. Therefore, p0 is directly proportional to the optical energy () absorbed by a target.

By exploiting different optical absorbance characteristics of different tissue types, spectroscopic PA imaging enables quantification of biological tissues. If the PA effects of the tissues are assumed to be independent, then a linear mixed model (LMM) can be applied [24], [25]. By convention, we hereafter refer the tissues of interest as endmembers, and we assume the optical fluence has been compensated. Consequently, if there are endmembers and different wavelengths (), the inverse problem of LMM is given aswhere is the design matrix, , is the -th wavelength, is the measured PA spectrum, is the unknown concentrations of the endmembers, and is additive measurement noises. Since x is non-negative, we solve the problem by non-negative least squares (NNLS) using quadratic programming techniques [26], which is given byColumns of A must be linearly independent to achieve a unique solution, i.e., the optical absorbance of any endmember cannot be expressed as a linear combination of the others. And must be zero-mean to make the least-squares problem valid.

Spectral system error in spectroscopic PA imaging

The assumption of a linear relationship between the PA spectrum and the unknown concentration does not always hold in practice due to the presence of system-dependent noise sources [25]. We herein conclude them with two system noise models (NMs).

NM 1. Laser-dependent noise model. The emitted laser energy is heteroscedastic from wavelength to wavelength in the spectral domain (NM 1a) and varies from pulse to pulse in the temporal domain (NM 1b), which altogether makes fluctuate in the spectrotemporal domain. To illustrate, the coefficient of variance (CVs) of in the near-infrared band (700 – 900 nm) on our Nd:YAG laser equipped with a tunable optical parametric oscillator (OPO) (Phocus Inline, Opotek Inc., USA) ranges from 3.65% (709 nm) to 9.85% (751 nm).

NM 2. Ultrasound system noise model. A PA imaging system for PCa imaging may resemble clinical ultrasound imaging systems [27], [28], [29], [30]. The Johnson-Nyquist noise is dominant in the ultrasound imaging systems among various other electronic noises such as flicker noise and Poisson noise. Generally, Johnson-Nyquist noise follows Gaussian distribution in a finite positive bandwidth. It is independent of either the optical source (NM 1) or the anatomy [31].

Given that and are proportional, a conventional optimization of spectroscopic PA imaging tackles NM 1 by normalizing the acquired PA signal to the correspondent laser pulse energy. However, such straightforward compensation makes NM 2 no longer independent of NM 1 – It endows the NM 2 with spectral signature of the inverted laser spectrum and heteroscedasticity. In other words, spectral unmixing accuracy now becomes wavelength-dependent because the statistics of ultrasound system noise such as mean and standard deviation (STD) also becomes wavelength-dependent after the compensation. It annuls the assumption of least squares that the additive noise must be zero-mean and homoscedastic [33], [34] Here we entitle the combination of these noise models as ‘spectral system error (SSE)’.

System-level optimization framework of PA spectral unmixing

Spectral system error (SSE) segregation scheme

Fig. 1 shows the overall framework of the optimized spectral unmixing for effective spectroscopic PA imaging. Here we incorporate the individual NMs into the spectral unmixing problem. The direct measured spectroscopic PA signal can be compensated to denoise NM 1:where is a diagonal matrix with for , and is the pulse energy at each wavelength . Once weighted, becomes proportional to the optical absorbance as in Eq. (1), legitimizing the least-squares model [32].

Fig. 1

The proposed framework of optimal spectral unmixing for effective spectroscopic PA imaging of PSMA expression. Noise models (NM) for the laser-dependent noise (NM1) and ultrasound system noise (NM2) are constructed and denoised in spectral unmixing. Oxyhemoglobin (HbO2), deoxyhemoglobin (Hb), and a PSMA-targeted agent [33] are spectrally unmixed. Simulation-aided wavelength selection and frame averaging are also considered in the framework for further enhancement of temporal resolution and precision.

We also denote the distribution of the electronic noise in NM 2 as , where is the mean value (), is a vector with all ones, and is the covariance matrix with all its elements on the diagonal being . Equivalently, we can decompose NM 2 into a constant offset and a zero-mean Gaussian component such that . Once the spectrotemporal fluctuation in NM 1 is compensated, the distribution of the electronic noise becomes , whereand

converts the diagonal of a matrix into a vector. By taking the weighted offset as one of the endmembers and rearranging the least-squares equation, we getwhere is the new design matrix, , and . We name this modification as ‘SSE segregation’ because it segregates SSE from primary endmembers by spectral unmixing. To convert the overall problem into an NNLS model solved by quadratic programming technique, we have

Frame averaging scheme

Frame averaging scheme can be an additional strategy for further SSE suppression, but more theoretical description is needed to understand its impact to the spectral unmixing accuracy. Assume subsequent frames are acquired at wavelength for all , such that the measurement now becomes . To denoise NM1, we normalize the frames by , where the diagonal of constitutes the recorded energy of the th frame at , . Then in frame averaging, frames are averaged to form and . The revised covariance of NM2 becomes

Moreover, the design matrix is not always perfect. For perturbation analysis, we assume there is a small and additive Gaussian perturbation to and denote it by , where each entry of is independent and identically distributed and satisfies , represents Gaussian distribution. Then after averaging frames, we have

It is well known that frame averaging enhances image SNR by times in raw PA data [34], [35], but how the SNR improvement affects the spectral unmixing accuracy is nontrivial. To analyze the unmixing accuracy with respective to , we first analyze the sensitivity of the least squares. Suppose is a small perturbation in , the sensitivity of least squares is given by [32]:

Here is the 2-norm condition number of , and are the largest and smallest non-zero singular values of ,, and . For brevity, we use and omit the “2” in the subscript to represent the Euclidean norm throughout. From the definition, the sensitivity is proportional to RMSE by a factor of . An effective way to lower the bound of the sensitivity is to minimize the noise-related term and by frame averaging. See Supplemental Information for the proof that the expectation and the STD of are proportional to . Nonetheless, having closed-form expression of the distribution of and the distribution of with respect to is intractable. Therefore, we validate this proof with a numerical evaluation.

Iterative wavelength optimization scheme

A wavelength optimization allows us to use fewer wavelengths to enhance temporal resolution meanwhile preserve spectral unmixing accuracy. Mathematically, the problem is defined as: given the number of rows of , find the best permutation of rows which minimizes the least-squares error. The optimization essentially minimizes the collinearity among endmember absorption spectrums, at the same time keeping wavelengths that suffer less from SSE.

A brute-force search is combinatorial in exponential computation time. To avoid such high time complexity, we propose an iterative searching scheme in two stages. In the first stage, the best wavelengths that are minimally required for unmixing are included in set ; then, in the second stage, one best wavelength is searched and added to at a time until the user-defined number of wavelengths (denoted as ) is reached. A detailed workflow is as follows, and a pseudo-code is provided in the supplementary information.

Input: number of endmembers , the total number of wavelengths , user-defined number of wavelengths , design matrix , weight , covariance , number of iterations .

Output: the optimal wavelength subset .

First stage: Findthat containswavelengths.

Initialize a counter with bins where each bin is attributed to a permutation.

Initialize a random fraction of endmembers . Generate PA signal with noises, such that .

Solve and calculate the root mean squared error (RMSE) for all permutations that have out of wavelengths, where and contain out of rows of and , respectively.

Find the permutation with the smallest RMSE and count one to its bin, where

Iterate steps (2) – (4) times, the permutation with the highest counts is reported as which is the optimal subset found in this stage.

Second stage: Updateone wavelength at a time.

Initialize a counter with bins, where is the cardinality of a set, and is the complementary set of (i.e., ).

Generate a random PA signal .

For each element , solve and calculate RMSE, where and contain the rows of and indexed by , respectively.

Find the element that yields the smallest RMSE and count one to its bin.

Iterate steps (2) – (4) times, then select the element with the highest count and include it into .

Iterate (1) – (5) until , where is the desired number of wavelengths.

The proposed approach is significantly lower in computational complexity than exhaustively searching all choose permutations. The number of NNLS operations is reduced from to .

Evaluations of the spectroscopic PA imaging framework

Calibration of the system noise models

The SSE due to NMs 1 and 2 was calibrated on our research platform. The spectral and temporal laser fluctuations in NM 1 were calibrated on a tunable Nd:YAG OPO pulsed laser system, where the laser was set in fast-scanning mode with 20-Hz pumping frequency, and each pulse energy was recorded by an energy meter (PE50BF-DIF-V2, Ophir Photonics, USA) in real-time. Five thousand sequential emissions at 720 nm were acquired to characterize the temporal laser fluctuation. Probability distribution of the PA intensity of the sequential emissions was represented as a histogram, and its similarity to a Gaussian distribution was compared. On the other hand, 512 sequential emissions per wavelength were recorded at wavelengths from 700 nm to 900 nm in 10 nm intervals to characterize the spectral laser fluctuation. The change in the mean and the STD of the pulse energy for the wavelengths were calculated.

NM 2 was calibrated from an L7–4 transducer connected to the Verasonics research package, where 1344 subsequent frames of PA images (in the size of 256 by 128, axial by lateral) were acquired to contain only electronic noise signals without any ultrasound or laser transmission. The mean and STD of all 44,040,192 pixels were calculated, and the intensity distribution was represented as a histogram and fitted to a Gaussian function. The spectral signature of the electronic noise () and the covariance of were derived from the NM 2.

Optimization of wavelength subset and SSE segregation

To optimize the wavelength selection, we implemented the algorithm in Section 2.3.3. We set the number of iterations to ensure statistical significance. The random spectroscopic PA signals were generated in a single pixel in the k-Wave simulator [36]. Several hyperparameters were defined for PA signal generation: (1) The ground-truth design matrix was comprised of the absorbance signatures of four endmembers: Oxyhemoglobin (HbO2) and deoxyhemoglobin (Hb) obtained from [37], a polyamidoamine (PAMAM) dendrimer-based PSMA-targeted contrast agent [33], and the SSE derived in Section 2.4.1. Specifically, the contrast agent absorption was calibrated by spectrophotometry in the near-infrared (NIR) range (700 nm to 900 nm) using a SpectraMax i3x multi-mode system (Molecular Devices, LLC, USA). (2) Fractions of all endmembers were randomized such that , where represents uniform distribution. Afterward, they were normalized to sum to one. (3) The covariance  of was obtained by calibration in Section 2.4.1.

Once the wavelength subset was optimized, another 1000 random spectroscopic PA signals were generated for testing. The mean and STD of the RMSE of all test signals were compared between four methods: proposed two-stage searching method, a metric-based method (MB) in [23] using product of singular values, randomized selection method, and a determined method selecting wavelengths based on uniform (UNI) interval within the given range ( wavelengths in the interval of (nm), counting down from 900 nm).

To randomize the endmember composition, we set the upper bound of the fraction of the SSE as such that , and fractions of all the other endmembers as . Once all fractions of endmember were generated, was normalized so that . To test the effectiveness of SSE segregation, we compared spectral unmixing accuracy with and without SSE segregation. We further varied in different values to investigate the impact of different levels of SSE to the spectral unmixing accuracy, which corresponds to the highly spatially variant SNR in the imaging field-of-view due to absorber distribution. In the evaluation of different wavelength selection methods, was fixed to 1.

Numerical analysis of frame averaging scheme on spectral unmixing accuracy

We conducted a numerical experiment to validate that the expectation and STD of and are proportional to , as mentioned in Section 2.3.2. The mean and STD of and , and the unmixing RMSE were calculated from 10,000 randomized spectroscopic PA signals. Ten values of were tested, i.e., 1, 4, 9, 16, 25, 36, 49, 64, 81, and 100. No frame averaging was deliberately implemented. Instead, we equivalently set the variance of and the covariance of as and to generate test signals from the numerical simulation, we show that if not strictly, it is highly approximated that the expectation and the standard deviation of and are all proportional to , such that the sensitivity of least squares is also proportional to .

In vivo imaging of PSMA expression

In vivo experiments were set up using the same ultrasound research package and Nd:YAG OPO laser system as described in 2.4.1. Bifurcated outlets of a fiberoptic bundle with rectangular apertures (0.88 mm × 40.00 mm × 2) were attached on the sides of an L7–4 linear ultrasound imaging transducer and co-registered in-plane, which conformed to conventional linear TRUS imaging configuration. Xenografts with isogenic PSMA-positive (PSMA+) PC3-PIP and PSMA-negative (PSMA-) PC3-flu cells were prepared on the lower back of five 6-to-8-week-old male mice. Note that a small amount of PSMA-targeted contrast was expected from the PSMA- PC3-flu (control) [20] which is a human-engineered PCa cell line having negligible PSMA expression. Validation of contrast agent uptake by PSMA+ PC3-PIP and PSMA- PC3-flu through tumor sectioning followed by fluorescence imaging was reported in [33]. Both tumors were imaged axially at around 2 cm deep identical to the ultrasound transducer’s elevation focus. Pre-injection data and 24-h post-injection data were analyzed. Because columns of , i.e., the characteristic spectrums of endmembers are in very different magnitudes, to enhance the numerical sensitivity and stability, we normalized column 1 (the contrast agent), column 2 (Hb), and column 4 (electronic noise) of with respect to that column’s maximum value. As for column 3 (HbO2), we normalized it with respect to the maximum of column 2, so that Hb and HbO2 can be visualized in the same dynamic range.

In this study, the spectral unmixing in vivo was evaluated by two metrics. (1) The first metric is the SNR of the PSMA-targeted contrast over the background, which is given bywhere and denote the mean PA intensities in the region of interest (ROI) of PSMA+ PC3 PIP tumor and the background ROI, respectively. (2) The second metric is the ratio of PSMA-targeted contrast agent intensity of PSMA+ PC3 PIP vs PSMA- PC3 flu tumor, given bywhere denotes the mean PSMA-targeted contrast agent intensity in the ROI of PSMA- PC3 flu tumor. In the pre-injection phase, ground truth , , and are all 0. However, due to the spatially uniform SSE, , , and are expected to be a constant. Hence, and are expected to be 1 in pre-injection. In the post-injection phase, ground-truth is 0, but ground truth and are not quantifiable in biological tissues. is expected to be a small constant due to SSE, and and are expected to be high and low, respectively, such that and are expected to be large values in post-injection. Higher in post-injection especially indicates better spectral unmixing sensitivity.

Results

Calibration of the noise models

The system-dependent NMs calibrated on our imaging system are shown in Fig. 2. Fig. 2a is an exemplary character of the laser instability (NM 1), where the heteroscedasticity, spectral fluctuation, and temporal fluctuation are reflected by varied CVs, mean values, and STDs. Particularly, energy of laser pulses at 710, 720, 800, 810, 850, 870, and 880 nm is relatively more stable (CV < 5%). Fig. 2b further explains the pulse-to-pulse temporal fluctuation, where a histogram of the energy of 5000 sequential pulses at 720 nm is demonstrated. The fitting curve in blue suggests that Gaussian distribution is a good approximation to the probability distribution (R2 = 0.930). Fig. 2c shows the calibration result of the background noise in NM 2. A histogram of image intensity in 44,040,192 pixels (256 by 128 per frame, 1344 sequential frames) is delineated, where the mean ± STD of all pixels is 2825 ± 1528. The fitting curve in blue also suggests that Gaussian distribution is a good approximation (R2 = 0.982).

Fig. 2

Characteristics of system NMs in the spectroscopic PA imaging. (a) Mean, STD, and the coefficient of variance (CV) of the Nd:YAG OPO laser at different wavelengths. (b) Energy distribution of 5000 sequential laser pulses emitted by an Nd:YAG OPO laser. (c) Background noise intensity histogram in the ultrasound imaging system.

Simulation-aided optimization of wavelength subset and SSE segregation

Wavelength optimization considering system factors

In this section, a simulation-aided wavelength optimization incorporating system-dependent NMs is presented. In the first stage where only minimally required number of wavelengths were selected, the proposed method found as {700, 710, 800, 880} (nm) with 115 counts and {700, 780, 800} (nm) with 3125 counts out of 100,000 iterations for the SSE segregation scheme and the conventional scheme, respectively. In either scheme, the selected wavelength subsets had significantly more counts than the other permutations, as indicated in Fig. S2. In comparison, MB method reported {700, 710, 800, 900} (nm) and = {700, 800, 900} (nm), respectively. In the second stage, the proposed method in the SSE segregation scheme sequentially included the following wavelengths into : 760, 830, 870, 900, 790, 720, 770, 820, 780, 740, 850, 730, 810, 860, 750, 890, 840 (nm). The whole process of the proposed method took 24.9 and 5.6 h respectively for the SSE segregation and the conventional scheme in MATLAB running on a laptop with Intel Core i7-8550U CPU with 16.0 GB RAM. If all permutations were searched brute-force in the same setup, the whole process would take around 352 days on the same machine, considering that a single permutation took 145 ± 26 µs.

We performed the comprehensive performance evaluation of four-wavelength selection methods (the proposed algorithm, MB method, random wavelength selection, and UNI selection) at different number of wavelengths with and without the SSE segregation scheme. Fig. 3 shows the mean spectral unmixing accuracy in RMSE of 1000 test signals.

Fig. 3

Mean spectral unmixing accuracy as a function of the number of wavelengths. One thousand signals were tested using the proposed method, MB method, random selection, and UNI selection with and without the SSE segregation scheme.

As indicated in Fig. 3, a whole-wavelength spectral unmixing yielded 0.040 ± 0.029 and 0.096 ± 0.047 of RMSE with and without the SSE segregation, suggesting 58.7% enhancement due to the SSE segregation scheme. In reference to the whole-wavelength result, we specifically analyzed the scenarios under two user preferences: (1) having the highest temporal resolution with the most compact wavelength subset or (2) securing a specific RMSE regardless of the size of wavelength subset.

In the first scenario, by using the most compact wavelength subset with the SSE segregation (4 wavelengths), the proposed wavelength optimization method yielded RMSE at 0.070 ± 0.049, which indicated more precise unmixing than other methods: 0.073 ± 0.053, 0.150 ± 0.094, 0.150 ± 0.091 for MB method, randomized, and UNI selection, respectively. Whereas in conventional scheme without the SSE segregation (3 wavelengths), the proposed method yielded RMSE at 0.113 ± 0.055 and still outperformed other methods: 0.119 ± 0.059, 0.215 ± 0.126, 0.122 ± 0.061 for MB method, randomized, and UNI selection, respectively. Collectively, in terms of the wavelength optimization, the proposed method showed the highest spectral unmixing accuracy compared with other methods in statistical significance (p < 0.0001). In addition, the SSE segregation scheme provided further fractional improvements at 38.6%, 38.6%, 30.1%, and − 23.9% over the conventional spectral unmixing cases without the SSE segregation. When the number of wavelengths is greater than four, SSE segregation shows lower spectral unmixing accuracy than the conventional scheme in UNI selection, which indicates the significance of a proper wavelength selection especially when the desired number of wavelengths is minimum. Nevertheless, SSE segregation scheme improves the accuracy in all other cases, and it only requires one more wavelength which decreases the temporal resolution from 6.67 scanning sequences per second to 5 scanning sequences per second at 20-Hz spectral tuning frequency.

In the second scenario, we consider the SSE segregation and set a specific RMSE threshold as the criterion for wavelength selection. When the expected RMSE is no larger than 120% of full-wavelength selection’s (21 wavelengths), at least 9, 10, 17, and 16 wavelengths were required by the proposed method, MB method, random selection, and UNI selection. Thus, the proposed method would yield 11.1%, 88.9%, and 77.8% higher temporal resolution than the MB, randomized, and UNI selection methods at a comparable spectral unmixing accuracy.

Robustness of wavelength selection and SSE segregation

In this section, we evaluate the sensitivity of spectral unmixing to . The level of denoted by was tuned into different values to simulate its nonuniform spatial distribution: 0.01, 0.05, 0.1, 0.2, 0.3, 0.5, 1, 2, and 4. The purpose is to test the robustness of SSE segregation toward spatially variant endmember composition, where in absorber-abundant tissue regions, takes up a low portion so is small and in background regions takes up a high portion so is large. Note again that the fractions of all other endmembers were randomized (i.e., ), and was normalized afterwards to sum to one. Each was tested by 1000 randomly generated signals with and without the SSE segregation scheme in different wavelength selection methods, as demonstrated in Fig. 4a.

Fig. 4

Spectral unmixing accuracy as the level of the fraction of the SSE changes. (a) RMSE of each wavelength selection method as a function of the number of wavelengths with and without the SSE segregation scheme. (b) The relative RMSE change of MB method, random, and UNI selection methods over what obtained with the proposed method. Gray lines indicate results given by different values of . Red lines are the baselines at 1. (c) The relative RMSE change of MB method, random, and UNI selection methods between scenarios with and without the SSE segregation scheme.

Vertical comparison of the curves in Fig. 4a reveals the effectiveness of SSE segregation, where the first row (noise segregation scheme) showed significantly-less change of RMSE than the second row (conventional scheme) as varies. Ratios of RMSE between the conventional and the SSE segregation schemes under different are given in Fig. 4b. Overall, the SSE segregation scheme becomes more effective as increases. However, when is negligible (e.g. which corresponds to a fraction of ), the conventional scheme outperforms the SSE segregation scheme mainly due to the sensitivity of least-squares. SSE segregation scheme is more effective at most of the numbers of wavelengths when , , , and for the proposed method, MB method, random selection and UNI selection, respectively, which corresponds to fractions of the electronic noise at 6.3%, 9.1%, 14.3%, and 14.3%.

Horizontal comparison of the curves in Fig. 4a reveals the robustness of different wavelength selection methods to . Each data point represents average RMSE of 1000 test signals given the number of wavelengths and . Fig. 4c further quantifies the how spectral unmixing accuracy varies as changes, which indicates the robustness of different wavelength selection methods toward the spatial variance of . The proposed method showed smaller STDs than the other methods regardless of the number of wavelengths, which implies higher robustness to the spatial variance of .

Effect of frame averaging to spectral unmixing accuracy

In Section 2.3.2, we provided numerical analysis on with respect to . Given that it is hard to express the expectation and the STD of and in closed form, we herewith implemented numerical simulations to find how these values changes with respect to .

In simulation, we set , and . Ten values of (i.e., 1, 4, 9, 16, 25, 36, 49, 64, 81, and 100) were tested by 10,000 data. Fig. 5a demonstrates line plots with error bars of and , and 5b shows the line plot with error bars of spectral unmixing accuracy in RMSE. The linear fittings included an extra datapoint at (0, 0), for the expectations it yielded with = 1.000, with = 1.000, and with = 0.997. As for the STDs, the linear fitting yielded that with = 1.000, with = 1.000, and with = 0.996, respectively. The fitting lines indicated a proportional relation of the mean and STD of and , and a highly approximated proportional relation of the mean and the STD of the spectral unmixing accuracy with respect to . Therefore, by assuming and are zero-mean Gaussian vector/matrix, the frame averaging scheme is expected to improve spectral unmixing accuracy by times.

Fig. 5

, , and the spectral unmixing accuracy in RMSE as functions of . Solid lines indicate the expectation of the values, error bars indicate the STD of the values, and dotted lines are the linear fittings.

Since both the number of wavelengths and the frame averaging determine imaging speed, we are interested in which one of them enhances spectral unmixing accuracy more with given temporal resolution. Fig. 6 demonstrates a comparison tested on 1000 random signals in simulation, where wavelength selection adds more wavelengths as the given temporal resolution decreases, and frame averaging employs minimally required number of wavelengths but increases the number of frames for averaging. Fig. 6 indicates that when selected wavelengths were optimal (proposed algorithm), wavelength selection slightly outperformed frame averaging. Otherwise, wavelength selection significantly outperformed frame averaging (random and UNI). This emphasizes the importance of selecting proper wavelengths prior to frame averaging. In addition, the gap of RMSE between wavelength selection and frame averaging was consistent in conventional scheme and SSE segregation scheme in all wavelength selection methods.

Fig. 6

Comparison of effectiveness of wavelength selection (wave sel) and frame averaging (frame avr) on spectral unmixing as the given temporal resolution changes. Each subfigure represents a wavelength selection method. Test results in conventional scheme (CON) and SSE segregation scheme were indicated by black and blue curves, respectively.

In vivo experiment

We first evaluated the effectiveness of the SSE segregation in in vivo data collected before and after injection of the PSMA-targeted contrast agent. Whole wavelengths were applied for both the conventional scheme (three endmembers) and the SSE segregation scheme (four endmembers) to preclude the impact of the wavelength selection methods. Fig. 7a demonstrates representative spectral unmixing outcomes of the PSMA-targeted contrast agent in vivo for both before and 24-h after injection. Compared with the contrast agent maps obtained from the conventional scheme, those by the SSE segregation scheme exceptionally suppressed background noises. Meanwhile, they revealed notably fewer artifacts in the tumor ROIs before injection, and preserved remarkable image contrast in the PSMA+ PC3 PIP tumor 24-h after injection. These observations confirmed that a substantial amount of SSE could be segregated from the primary endmembers (i.e., contrast agent, Hb, and HbO2). Due to the sensitivity of least squares, the segregated SSE map presented a highly uniform spatial variance outside the body but less uniform inside the body where its fraction was negligible, which well agrees with our simulation results in Section 3.3.

Fig. 7

In vivo evaluation of the SSE segregation scheme in the whole-wavelength method utilizing a spectral range from 700–900 nm in 10 nm interval in NOD-SCID mice bearing subcutaneous PSMA+ PC3 PIP and PSMA- PC3 flu in the lower back right and left posterior flanks, respectively. (a) Representative in vivo ultrasound and spectral unmixed map of PSMA-targeted contrast agent in conventional and SSE segregation schemes, pre-injection and post-injection phases; Dashed white lines represent ROIs in PSMA+ PC3 PIP and PSMA- PC3 flu tumors. (b) and (c) are mean and standard error of and , respectively, in conventional (CON) and SSE segregation schemes, pre-injection and post-injection phases (n = 5). Red dashed lines are baselines at 1.

A quantitative analysis was performed to evaluate the effectiveness of the SSE segregation scheme in different ROIs, where and in conventional scheme and SSE segregation scheme were calculated using data acquired before and after injection of the PSMA-targeted contrast agent (n = 5). As shown in Fig. 7b, in pre-injection, the mean and standard error of in conventional scheme and SSE segregation scheme are 0.29 ± 0.06 and 1.56 ± 0.51, respectively. In conventional scheme, is lower than 1 because in some pixels system noise was more decomposed into other endmember maps (Hb and HbO2) than the contrast agent map, so that intensity in the body region is lower than the background region, as revealed in Fig. 7a. Contrarily, in SSE segregation scheme, is close but greater than 1 because system noise can be effectively segregated in the background region, but less effectively when strong optical absorber exists in the body region. In these body regions, the fraction of system noise is too low to be unmixed by spectral unmixing due to the sensitivity of least squares, as indicated in the noise map in Fig. 7a. In post-injection, the mean and standard error of in conventional scheme and SSE segregation scheme are 2.07 ± 0.30 and 9.12 ± 1.99, respectively. The mean in SSE segregation is 4.43 times of the one in conventional scheme, indicating that SSE segregation is significantly effective in improving SNR. As shown in Fig. 7c, in pre-injection, the mean and standard error of in conventional scheme and SSE segregation scheme are 0.77 ± 0.10 and 0.77 ± 0.11, both of which are slightly lower than 1. In post-injection, the mean and standard error of in conventional scheme and SSE segregation scheme are 6.79 ± 2.15 and 14.77 ± 7.76, respectively. The mean in SSE segregation is 2.18 times of the one in conventional scheme. The enhancement of is not as high as , possibly because Hb and HbO2 contributed false contrast or there was minor contrast agent uptake in PSMA- PC3 flu region. Nevertheless, it suggests that SSE segregation can greatly improve sensitivity of prostate cancer detection. Note that due to the lack of ground truth of Hb and HbO2 in vivo, we only conducted quantitative analysis on the unmixed contrast agent map. However, SSE segregation works on all unmixed endmembers. It can also effectively reduce background noise in the unmixed Hb and HbO2 map as demonstrated in Fig. S2.

Wavelength optimization scheme

We expanded our analysis to the wavelength optimization scheme using in vivo data before the injection of PSMA-targeted contrast agent (n = 5), given that ground truth of pre-injection data is known. In this section, proposed method, MB method, random selection, and UNI selection were analyzed in SSE segregation scheme.

Representative PSMA-targeted contrast agent maps and quantitative evaluation of and are demonstrated in Fig. 8. In Fig. 8a, the contrast agent maps were unmixed by four and eleven wavelengths, which respectively represent scenarios aiming at the highest temporal resolution and twice the temporal resolution of using whole wavelengths. In visual assessment, the proposed scheme demonstrated a substantial suppression of false contrast compared with random and UNI selection methods. Artifacts indicated by white arrows were found around the tumor regions even in whole-wavelength selection but not in the proposed method and MB method. Fig. 8b presents a quantitative evaluation of as a function of number of wavelengths, where the proposed method outperformed other methods in either accuracy or robustness. In particular, the mean and standard error of at four wavelengths were 1.20 ± 0.28, 2.19 ± 0.74, 4.80 ± 3.20, and 1.96 ± 0.62, for proposed method, MB method, random selection and UNI selection, respectively. Whole-wavelength selection yielded 1.33 ± 0.40. Fig. 8c reveals quantitative evaluation of as a function of number of wavelengths. The proposed method had a comparable performance with MB method, while outperforming random selection and UNI selection methods. Specifically, mean and standard error of at four wavelengths were 1.07 ± 0.27, 0.94 ± 0.19, 1.11 ± 0.14, 1.20 ± 0.24, for proposed method, MB method, random selection and UNI selection, respectively. And whole-wavelength selection yielded 0.89 ± 0.14. In all, the proposed method provided less false contrast of PSMA-targeted contrast agent before the injection of the contrast agent in all ROIs regardless of the number of wavelengths.

Fig. 8

Comparison of the wavelength selection methods in data acquired before injection of the PSMA-targeted contrast agent in NOD-SCID mice bearing subcutaneous PSMA+ PC3 PIP and PSMA- PC3 flu in the lower back right and left posterior flanks, respectively. (a) Representative spectral unmixed contrast agent maps using four wavelengths and eleven wavelengths. (b) and (c) are and as functions of the number of wavelengths. Red dashed lines are the baseline at 1.

We validated the impact of the frame averaging scheme on the spectral unmixing accuracy using in vivo data collected before injection of the PSMA-targeted contrast agent (n = 5). Note again that any intensity in the unmixed PSMA-targeted contrast map in the pre-injection phase is false and directly reflects absolute spectral unmixing error. A representative image is shown in Fig. 9a, where a PA image acquired at 780 nm and averaged for 64 times is demonstrated, and different values were tested for frame averaging: 1, 4, 9, 16, 25, 36, 49, and 64. In visual assessment, apparent trend of decreasing false intensity was found as increases.

Fig. 9

Evaluation of the relation between spectral unmixing accuracy and the number of frames for averaging () in vivo. (a) Representative contrast agent maps when 1, 4, 9, 16, 25, 36, 49, and 64 frames were averaged for spectral unmixing. (b) The mean and STD of false contrast in the ROIs of PSMA+ PC3 PIP, PSMA- PC3 flu, and background as functions of . The ROIs are indicated in (a) by dotted squares.

A quantitative analysis was further conducted. For 1, 4, 9, 16, 25, 36, 49, 64 frames of averaging, respectively, are 0.77 0.11, 0.71 0.27, 0.80 0.21, 0.84 0.36, 0.94 0.35, 1.07 0.29, 1.08 0.35, and 1.16 0.32, and are 0.98 0.11, 0.96 0.16, 0.91 0.20, 0.85 0.19, 0.87 0.31, 0.91 0.34, 0.89 0.34, and 0.92 0.35. Two relations were further verified on the images shown in Fig. 9: proportion of expected spectral unmixing accuracy vs. and proportion of STD of spectral unmixing vs. . Square-shaped ROIs in the same size (3 × 3 mm2 containing 200 pixels) were selected from the PSMA+ PC3 PIP, PSMA- PC3 flu tumors, and background regions (Fig. 9a) to keep an identical sample size for the following statistics. The mean and STD of false intensity within each ROI were plotted as a function of in Fig. 9b. The linear fitting of the mean false intensity suggested with in PSMA+ PC3 PIP; with in PSMA- PC3 flu; and with in the background, respectively. The linear fitting of the STD of false intensity suggested , in PSMA+ PC3 PIP; , in PSMA- PC3 flu; and , in the background, respectively. The results suggest strong linear relationship between spectral unmixing accuracy to with negligible bias in y-intercept. In addition, similar slopes at different ROIs indicated no substantial spatial variance of the spectral unmixing accuracy.

Discussion and conclusion

Here we presented a system-level optimization of the spectroscopic PA imaging to detect aggressive PCa using PSMA-targeted agent. The optimization framework consists of three schemes: SSE segregation, wavelength optimization, and frame averaging. We first constructed system-dependent NMs for laser and ultrasound systems, and incorporated the resultant ‘spectral system error’ into spectral unmixing. We then considered these NMs in the wavelength optimization, where an iterative searching method was proposed. We further analyzed the relation between the spectral unmixing accuracy and the number of frames for averaging.

To validate the framework, we first calibrated NMs on our imaging system. Gaussian distribution well described the temporal fluctuation of laser pulse energy in NM 1 and the electronic noise intensity distribution in NM 2. Even though the distribution of NM 2 was closer to a chi-square distribution, which is a sum of squares of two or more independent Gaussian variables, it was close enough to a single Gaussian such that the assumption of least squares was not annulled.

The proposed wavelength selection method showed superior spectral unmixing accuracy and robustness in our simulation-aided optimization. It is worth mentioning that MB method indicated comparable accuracy in in vivo validation, and it was less computationally complex. However, it showed lower robustness than the proposed method as it does not consider system factors. Moreover, the computation time is not a primary concern for wavelength selection if the algorithm's complexity is reasonable (e.g., timed in hours), because wavelengths selection is predetermined and will not be performed on a daily basis by clinicians [41], [42].

We are particularly encouraged by SSE segregation, which enhanced SNR and spectral unmixing accuracy by segregating background electronic noise. The SSE segregation reduced spectral unmixing errors by 58.7% in simulation (n = 1000) and enhanced the contrast between PSMA+ PC3 PIP and PSMA- PC3 flu tumors by 2.18 times, as well as the image SNR by 4.43 times in 24-h after injection in vivo data (n = 5). It suggests that the SSE segregation is efficacious in suppressing false contrast and improving specificity especially when detecting early lesion with low PSMA expression. In terms of SNR, note that one can possibly apply ultrasound-based delineation for co-registered PA images to get rid of background regions if the anatomy has a clear boundary in ultrasound images [38]. However, ultrasound-based delineation is not always applicable to clinical applications. For instance, standalone transrectal ultrasound provides low sensitivity (11–35%) and ambiguous morphologic feature to localize prostate cancer due to poor contrast resolution [39]. Moreover, SSE segregation does not conflict but further supplement ultrasound-based delineation because it is part of spectral unmixing before any segmentation. What is more important is SSE segregation segregates the system noise from other endmembers to avoid false contrast from being unmixed to the ROI. Despite of enhancements, profound analysis unveiled that SSE segregation is spatially variant. This is because in regions where the fraction of the unknown variable is too low, linear spectral unmixing is not capable to distinguish it from the additive noise in the measurement due to the sensitivity of least squares. Simulation suggested that SSE segregation instead lowered spectral unmixing accuracy when the electronic noise fraction was less than around 3%. In vivo validation also confirmed the spatial variance of SSE segregation, that intensity of segregated electronic noise was more uniform outside the body than inside. To tackle that challenge, one may consider an adaptive approach for spectral unmixing: the electronic noise would be included as an endmember in the initial spectral unmixing, but excluded afterwards if its unmixed fraction is lower than the threshold (e.g., 3%). Moreover, instead of implementing spectral unmixing pixelwise, spatial constraints such as total variation regularization can be added to minimize the spatial variance. Nevertheless, a regularization will increase the time complexity and potentially impair temporal resolution.

Finally, we quantified how frame averaging improves spectral unmixing accuracy in theoretical analysis and numerical validation. Assuming the system electronic noise is composed of a constant value plus a Gaussian noise , where is zero-mean and its entries are all independent, our study suggested that by having frames, frame averaging reduces spectral unmixing error by times. Note that we were assuming the perturbation in the design matrix and the measurement noise are all Gaussian noises, otherwise the quantitative relation between the spectral unmixing accuracy and the frame averaging remains an open question. Nevertheless, this relation is very useful for the trade-off between the imaging precision and the imaging speed depending on the application. Specifically, if the temporal resolution is fixed, this relation can be applied in simulation to decide whether to have more frames per wavelength or fewer frames but more wavelengths for optimal spectral unmixing accuracy, as demonstrated in Section 3.2.3. It turns out given a fixed temporal resolution, a proper wavelength selection is critical to be combined with frame averaging. When selected wavelengths are optimal, adding more wavelengths is slightly more effective in enhancing spectral unmixing accuracy than having more frames for averaging. Also note that frame averaging improves the effectiveness of SSE segregation, because it narrows the probability distribution of as shown in Fig. 2c, so has a less varying spectrum and it can be better unmixed as an endmember. Nevertheless, in the in vivo validation, increases from 0.77 to 1.16 as the number of frames for averaging increases from 1 to 64, while fluctuates around 0.9 and does not seem to have an obvious trend as the number of frames increases. This suggests that frame averaging is more effective for the background region, and again indicates the spatial variance of spectral unmixing due to the sensitivity of least squares. In background regions, SSE is the only composition, so after frame averaging it can be more easily unmixed. Whereas in the body region where Hb and HbO2 are the dominant contrast, SSE and the contrast agent may fall below the sensitivity threshold, thus they are less easily to be distinguished and frame averaging becomes less effective.

In all, our system-level optimization of the spectroscopic PA imaging of aggressive PCa successfully segregated the SSE while optimizing the temporal resolution. Note that we here focused on the system-dependent noise sources, and did not consider tissue-dependent spectral error in this study (e.g., spectral coloring, motion artifacts, etc. [24]). Including these factors in simulation-aided optimization may be time-consuming and arbitrary due to high anatomy-dependency. Nevertheless, for precision purposes, one can consider depth-resolved and model-based fluence compensation [40], [41], [42] or hardware solutions [43]. We will also further investigate whether a more realistic covariance matrix of the Gaussian process can approximate the fluence distribution. Once thousands or more ground-truth images are generated, they can be used for data-driven methods such as deep learning [44] to further enhance wavelength optimization. The proposed wavelength selection is a prototype to deep learning methods – the best combination of wavelengths was first ‘learned’ from 100,000 one-dimensional signals, then tested on real images. Also, other nonlinear spectral unmixing approaches studied in hyperspectral imaging [45] may provide additional insights.

Funding

This work was supported by the [CA134675; CA183031; CA184228; EB024495; R01HL139543]; . Yixuan Wu was supported by NIH Graduate Partnerships Program. Dr. Jeeun Kang was supported by the , United States Department of Defense [W81XWH-18-1-0188]. Dr. Emad Boctor was supported by NSF CAREER AWARD [1653322].

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.