Wavelength optimization in the multispectral photoacoustic tomography of the lymphatic drainage in mice.

Multispectral photoacoustic tomography provides mapping of the tissue chromophore distributions using sets of tunable laser wavelengths. With the overall goal of studying the dynamics of cerebrospinal fluid in mice in vivo, our work aims to minimize the number of wavelengths to reduce scanning time, improve the temporal resolution, reduce the energy deposition and avoid the tracer photobleaching while maintaining high image quality. To select small sets of wavelengths we directly searched for the combinations of wavelengths providing the best and worst image quality in comparison with a reference image obtained using 131 closely spaced wavelengths between 680 and 940 nm in terms of the peak signal-to-noise ratio (PSNR). We have shown that using the PSNR optimization method, additional improvements could be achieved over the wavelength set selected using the method of the minimization of the extinction matrix condition number.

Introduction

Multispectral photoacoustic tomography (MPAT) is a rapidly evolving imaging modality, which offers improved resolution of tissue chromophores due to the combination of high resolution of ultrasound imaging with the high molecular specificity resulting from differences in the light absorption by different molecules at different near-infrared wavelengths [1,2]. Due to its ability to facilitate functional, molecular, anatomical and oxygen metabolism expression, MPAT is a widely used tool for many applications among which are the ophthalmology, cardiology, neurology, vascular biology, molecular imaging, multiscale functional imaging, whole-body imaging of small animals, real-time dynamic imaging using tracers [[1], [2], [3], [4]], and in cancer studies [5,6].

Lymphatic vascular system, which comprises a network of vessels extending to almost every part of the body, drains the fluid from the interstitial tissue to regional lymph nodes and then to blood system to maintain overall fluid balance. Lymphatic vessels transport this clear fluid against a pressure gradient between tissue and venous system with an active pumping system. The lymphatic system is also involved in the metastatic spread of cancers. The discovery of the lymphatic vessels in the human eye [7] opened the door to further investigations of ocular lymphatic drainage [8]. The mouse appears to be an adequate study model due to the similarity to human aqueous dynamics and pharmacology [8]. Using near-infrared contrast agents, MPAT can be used in vivo to map the active lymphatic drainage in mice [8]. In this study, to trace the lymphatic drainage from the eye subconjunctival region into the neck lymph nodes of the mouse, we used the near-infrared dye QC-1. The ultimate goal of the study was to ensure the accurate separation from the endogenous chromophores (such as oxy- and deoxyhemoglobin) and mapping of the NIR tracer using a minimum number of light wavelengths. Minimization of the number of wavelengths was important for this goal because it increased the temporal resolution, reduced the motion artifacts and the tissue energy deposition and decreased the photobleaching effect i.e. the irreversible fading of the dye contrast during measurements [9].

In terms of the wavelengths optimization in MPAT, Yuan et al [10] proposed using the approach earlier developed for the diffuse optical tomography of tissue [11], which was based on the minimization of the condition number of the matrix of the extinction coefficients of tissue chromophores. Luke & Emelianov [12] further examined this method along with two other methods for reduced wavelengths sets selection: they compared the evenly separated wavelengths and the wavelength set selected based on the maximum singular value () of the extinction matrix and the set selected based on minimum condition number () algorithm [12]. Among the three methods, the set of three wavelengths selected based on the maximum singular value () as well as the set based on the minimum condition number showed better results than the set of evenly separated wavelengths [12].

While the idea of using the minimum condition number approach was based on the general theory of the optimal least square fitting [11], in this study we directly found the reduced wavelength sets producing the best quality of chromophore images (oxy and deoxy-hemoglobin and the tracer) in comparison with the images obtained using 131 wavelengths between 680 nm and 940 nm using the maximum peak signal-to-noise ratio criterion [13].

Methods

The animal model and protocol

The study was performed after the approval by the St. Michael's Hospital’s Animal Care Committee. Images were acquired using a commercially available MPAT system (MSOT inVision 128, iThera Medical GmbH, Germany) equipped with the tunable 532 nm Nd:YAG solid-state laser that pumped an optical parametric oscillator. One microliter of QC-1 (1 mM, LI-COR Biosciences, USA) conjugated with Bovine Serum Albumin (Sigma-Aldrich, USA) was injected into subconjunctiva of the right eye of CD1 albino mice (Charles River Laboratories, Canada) (n = 4). The animal preparation was performed as previously described [8]. The mouse (∼30 g body weight) was wrapped in the waterproof plastic wrapping and placed into the imaging chamber filled with the warmed (34 °C) water. The image was acquired 2 h after the tracer injection into the subconjunctival space of the right eye. At that time the molar concentration of QC-1 compared to the blood chromophores was on the order of ∼1/106.

MPAT imaging

The ultrasonic signals from the array of 128 ultrasound transducers placed around the tissue as shown in Fig. 1, were digitized with double precision accuracy, at10 shots for averaging, and then used to determine the initial acoustic distribution of the source, and to map the light absorption properties at 131 wavelengths between 680 nm and 940 nm with the step of 2 nm. The set of wavelengths was chosen between 680 nm and 940 nm to avoid the effects of the water absorption. The acquisition time per each wavelength was 1 s per slice (position), the total acquisition time for a set of 131 wavelengths was 131 s per slice. The back projection image reconstructions for each wavelength at 332 × 332 pixels resolution were performed using the native ViewMSOT™ software (iThera Medical GmbH, Germany). The further multispectral image deconvolution using the linear regression algorithm produced images of the distributions of three chromophores: oxyhemoglobin (HbO2), deoxyhemoglobin (Hb), and the tracer QC-1.

Fig. 1

Photographs of experimental setup. (a) Photograph of ultrasound transducer ring (red arrows) and laser fibers (yellow asterisks). (b) Photograph of animal preparation. An anesthetized albino mouse is placed in a specialized mouse holder (yellow arrow). (c) View from the top of the mouse (red arrow) in the holder placed in the MSOT scanning chamber.

Chromophore spectra

Fig. 2 shows near-infrared spectra of oxy- and deoxyhemoglobins and of the lymphatic tracer QC-1. The hemoglobin specific absorption (molar extinction) spectra were taken from [14]. A near-infrared quencher dye QC-1 was used to trace the lymphatic drainage. The specific absorption of the tracer QC-1 was taken from [15]. The maximum absorption of QC-1 was at 737 nm which was very close to the local minimum of Hb and low absorption region of HbO2. Thus HbO2, Hb, and QC-1 dominated light absorption at different regions without severe overlapping of their absorption maxima.

Fig. 2

Absorption spectra of oxyhemoglobin (HbO2), deoxyhemoglobin (Hb) and of QC-1 tracer according to [14,15]. QC-1 had maximum absorption at 737 nm, where Hb had minimum absorption, while 756 nm corresponded to the local maximum of Hb absorption. 800 nm corresponded to the isosbestic point where Hb and HbO2 had equal specific absorption.

Spectral deconvolution and image rendering

Spectral deconvolution was performed using ViewMSOT™ software and custom MATLAB® (MathWorks, Inc., USA) code, both based on the linear regression [16]. In particular, in our MATLAB® code we used the standard regress function. The chromophores included in our spectral model were oxy-, deoxy-hemoglobin and QC-1, which is sufficient for the albino mouse. Using a custom MATLAB® code, the images were thresholded and the Hb and HbO2 chromophore crosstalk filters were applied to remove oxy-, deoxy-hemoglobin cross-talk signals from the tracer images.

Fig. 3 explains the importance of the accurate mapping of both the tracer and hemoglobin. Fig. 3(a) shows a typical chromophore map of the coronal sections of the mouse neck two hours after the tracer (QC-1) injection into the subconjunctiva of right eye obtained using 131 wavelengths between 680 nm and 940 nm. The bright green color shows the tracer in the right neck lymph node (shown by the arrow). Red and blue colors show HbO2 and Hb in the neck arteries and veins, respectively. In actuality, due to the low concentration of the tracer and minor spectral cross-talk with the hemoglobin the tracer may spuriously appear not only in the lymphatic vessels but also in the blood vessels. Fig. 3(b) shows only the tracer distribution. Since the goal of our study was mapping the lymphatic vessels, the pixels corresponding to the blood vessels in Fig. 3(b) were identified using the HbO2 and Hb images, and their values were set to zero so that only the lymph node appears bright in Fig. 3(b).

Fig. 3

(a) ViewMSOT™ image of the coronal sections of the mouse neck two hours after the tracer (QC-1) injection into the subconjunctiva of right eye obtained using 131 wavelengths between 680 nm and 940 nm. R(right),L(left),S(superior) and I(inferior). The bright green color shows the tracer in the right neck lymph node (shown by arrow). Red and blue colors show HbO2 and Hb in the neck arteries and veins, respectively. (b) A separate tracer map reconstructed using MATLAB® and the same set of 131 wavelengths as in (a). In (b) the signals corresponding to the locations of the blood vessels were removed.

Wavelength optimization and image quality assessment

The peak signal-to-noise ratio (PSNR) and the structure similarity metric (SSIM) were introduced in [13] as the error sensitivity metrics and the structural similarity, respectively, of an image against a reference image. For two identical images PSNR tends to infinity and SSIM tends to 1 [13]. To find the initial three - wavelength PSNR-optimized sets we calculated PSNR for the images obtained using all possible combinations of three wavelengths out of 131, and the reference image was the one obtained using all 131 wavelengths. We also found a three-wavelength set using the minimization of the condition number (CNM) of the chromophore extinction (molar absorption coefficient) matrix [[10], [11], [12]] using our custom MATLAB® code. This set was not the same as in [12] though, since in our case one chromophore was QC-1 while in [12] the dye was MMPSense 750 FAST, Perkin Elmer. We also used PSNR and SSIM to measure the quality of images reconstructed using additional sets of four, five, and six wavelengths obtained by combining various wavelengths obtained using the PSNR and CNM methods. Since the final image data format can affect both PSNR and CNM values, we analyzed these measures for the unsigned 8-bit integer (the coarsest) and double precision (most precise) images. The peak value for the 8-bit images was 255, and for the double precision images it was taken as a maximum pixel value of the image.

Results

Table 1 showed all wavelength sets included in our study. MCN minimization for the three-chromophore extinction matrix yielded a set of 3 wavelengths (Set 1). PSNR optimization was first performed for each of five slices, and then averaging was applied to determine three sets of wavelengths corresponding to best PSNRs for QC-1, HbO2 and Hb (Sets 2, 3, 4). One can see that each chromophore had its own optimal wavelength set different from the others.

Table 1

Wavelength sets compared in the study.

#	Method	Wavelengths, nm
1	Minimum condition number (MCN)	680, 724, 930
2	QC-1 PSNR optimization	691 ± 8, 731 ± 8, 913 ± 3
3	Hb PSNR optimization	680 ± 1, 828 ± 7, 886 ± 3
4	HbO2 PSNR optimization	728 ± 9, 842 ± 13, 905 ± 6
5	Combination of PSNR- optimized and MCN	684,728, 852, 916
6	Combination of PSNR- optimized and MCN	686, 730, 836, 896, 922
7	Combination of PSNR- optimized and MCN	680, 700, 730, 840, 886, 930
8	Recommended by iThera	700,730,760,800,850
9	Close to minimum average PSNR	724, 782, 882
10	Minimum average PSNR	726, 784, 884
11	Spectral features of chromophores	736, 756, 800

In addition to the optimized three-wavelength sets 1–4 we also created four-, five-, and six-wavelength sets (sets 5–7), by partitioning all different wavelengths found using PSNR and MCN optimizations into four, five and six clusters, so that the wavelengths in each set corresponded to the cluster centroid locations [17]. In addition, we included the five-wavelength set recommended by iThera Medical GmbH (set 8). Finally, the three-wavelength sets 9 and 10 were selected based on the minimum average PSNR for Hb, HbO2 and QC-1 maps. Set 10 corresponded to the absolute minimum of the modulus of the three-wavelength extinction matrix determinant. In addition, we also considered set 11 of wavelengths corresponding to the spectral features indicated in Fig. 2. Note that the wavelengths in “bad” sets 9, 10, and 11 (according to the values in Table 2) were well separated similar to those in “good” sets 1–4.

Table 2

PSNR and SSIM values averaged for five slices and two (8-bit images) or three (double precision images) chromophores.

Set #	8-bit images (Hb and HbO2, no QC-1)		Double images (Hb,HbO2 and QC-1)
	Average PSNR (dB)	Average SSIM	Average PSNR (dB)	Average SSIM
1	54.4 ± 3.8	0.994 ± 0.003	29.7 ± 5.4	0.97 ± 0.02
2	56.3 ± 3.2	0.996 ± 0.001	30.8 ± 4.1	0.97 ± 0.02
3	56.2 ± 2.7	0.996 ± 0.002	29.4 ± 3.4	0.97 ± 0.02
4	54.1 ± 3.4	0.989 ± 0.004	27.6 ± 5.2	0.96 ± 0.04
5	58.7 ± 3.3	0.998 ± 0.002	33.3 ± 4.0	0.98 ± 0.02
6	57.7 ± 3.8	0.997 ± 0.001	33.3 ± 4.8	0.98 ± 0.02
7	57.1 ± 3.4	0.996 ± 0.002	31.9 ± 4.2	0.97 ± 0.02
8	55.8 ± 3.0	0.997 ± 0.003	29.7 ± 3.5	0.97 ± 0.02
9	43.2 ± 14.6	0.976 ± 0.026	15.3 ± 15.2	0.94 ± 0.05
10	31.4 ± 11.2	0.958 ± 0.010	1.0 ± 11.5	0.89 ± 0.04
11	50.7 ± 2.93	0.989 ± 0.006	18.9 ± 13.1	0.95 ± 0.02

Table 2 showed PSNR and SSIM values for each of the elleven sets averaged for six slices and chromophores. In addition, Table 3 showed PSNR and SSIM values of each chromophore for each of the elleven sets averaged for six slices.

Table 3

PSNR and SSIM values of Hb, HbO2, and QC-1 images averaged for five slices.

Set #	Double images Hb		Double images HbO2		Double imagesQC-1
	AveragePSNR (dB)	AverageSSIM	AveragePSNR (dB)	AverageSSIM	AveragePSNR (dB)	AverageSSIM
1	28.8 ± 1.8	0.95 ± 0.02	24.3 ± 1.3	0.97 ± 0.02	36.0 ± 3.4	0.98 ± 0.01
2	28.9 ± 1.3	0.95 ± 0.01	27.9 ± 1.8	0.98 ± 0.01	36.7 ± 3.0	0.98 ± 0.01
3	30.8 ± 1.9	0.96 ± 0.02	25.9 ± 2.0	0.97 ± 0.01	31.3 ± 1.3	0.96 ± 0.01
4	21.2 ± 1.4	0.91 ± 0.01	31.2 ± 2.0	0.99 ± 0.01	30.5 ± 3.5	0.96 ± 0.01
5	31.0 ± 1.5	0.96 ± 0.01	32.0 ± 0.6	0.99 ± 0.01	36.7 ± 3.8	0.98 ± 0.01
6	31.9 ± 0.9	0.96 ± 0.01	30.1 ± 1.6	0.99 ± 0.01	37.9 ± 1.4	0.98 ± 0.01
7	30.6 ± 2.0	0.96 ± 0.01	29.5 ± 0.8	0.99 ± 0.01	38.2 ± 4.2	0.99 ± 0.01
8	29.8 ± 1.8	0.95 ± 0.01	28.9 ± 1.1	0.99 ± 0.07	36.8 ± 3.3	0.98 ± 0.01
9	29.1 ± 1.5	0.95 ± 0.01	26.8 ± 2.1	0.98 ± 0.01	33.1 ± 3.1	0.97 ± 0.01
10	−11.7 ± 5.0	0.90 ± 0.03	14.1 ± 2.8	0.92 ± 0.02	0.70 ± 3.7	0.84 ± 0.05
11	11.9 ± 2.7	0.93 ± 0.02	10.2 ± 3.4	0.97 ± 0.02	34.6 ± 1.4	0.95 ± 0.02

The 8-bit QC-1 images for all sets except set 10 had infinite PSNR and SSIM equal to one. For this reason, the PSNR and SSIM values corresponding to 8-bit images were averaged only for Hb and HbO2, while the values corresponding to the double precision images were averaged for all three chromophores. Fig. 4 shows PSNR values from Table 2 sorted from largest to smallest for both unsigned integer and double precision images. One can see that in terms of relative PSNR results for both 8-bit and double precision images were similar within standard deviations.

Fig. 4

Averaged PSNR values for all sets: (a) unsigned 8-bit Hb and HbO2 images; (b) double precision Hb, HbO2 and QC-1 images. Error bars are the standard deviations.

PSNRs for all optimized three-wavelength sets were close within their standard deviations. This showed that MCN method provided a relatively good set of three wavelengths (set 1), although not the best of all (set 2 had a slightly higher PSNR). The highest PSNRS were found for sets 5 (four wavelengths) and 6 (five wavelengths), obtained by partitioning wavelengths initially found by the three-wavelength optimization methods. The best sets 5 and 6 provided PSNR values approximately 10% higher than set 1. The 6-wavelength set 7 showed lower PSNR than sets 5 and 6. The five-wavelength set 8 recommended by the scanner vendor produces images of the same quality as the three-wavelength optimized sets. The spectral feature-based three-wavelength set 11 had poor PSNR and SSIM values compared to sets 1-9. The worst set 10 had PSNR close to 30 time smaller than set 1 for double precision images and about half of set 1 PSNR for the 8-bit images. The main distinctions between PSNR values for double precision and 8-bit images were larger values for 8-bit images and larger differences between different sets for the double precision images.

To visually assess the results in terms of the image quality, the 8-bit chromophore maps obtained using various wavelength sets with ViewMSOT software were shown in Fig. 5 along with the reference image, which was obtained using 131 wavelengths.

Fig. 5

8-bit chromophore maps obtained using different wavelength sets with ViewMSOT software. Red, blue, and green show oxyhemoglobin, deoxyhemoglobin, and QC-1, respectively. The arrow shows the tracer concentrated in the right neck lymph nodes of the mouse; L,R,S,I indicate left, right, superior, and inferior sides.

One can see that the “good” wavelength sets 1, 5, and 8 produced images structurally close to the reference one with main differences in the hemoglobin distribution rather than in the tracer locations. The image from the “best” set 5 shows slightly better similarity with the reference image in terms of the deoxyhemoglobin distribution than images from other sets. The image for the “worst” set 10 shows no similarity with the reference, and the image for the “bad” set 9 shows some poor similarity in terms of the oxy-hemoglobin distribution and no similarity in terms of other chromophore distributions.

Discussion

Unlike previous similar studies relying on the CNM method to minimize the number of wavelengths for MPAT [[10], [11], [12]], in this study we found the best wavelength sets for imaging specific chromophores by a direct search through all possible three-wavelength combinations to find the combinations providing the highest and the lowest PSNR values. Similar to the previous work [12] we used the images obtained using a large set of 131 closely spaced wavelengths as a reference image. We did not consider the aspect of wavelength optimization from the standpoint of the optimal optical fluence evaluation within the blood or lymphatic vessels [18], which will be addressed in our future work.

Apart from PSNR values, we also calculated SSIM values, which provided another measure of the comparison of two images. All values were calculated for five slices covering the mouse neck providing the samples for the statistical analysis. Since images could be exported in formats with different precision levels between 8-bit unsigned and double precision, we calculated PSNR and SSIM values for those two extreme cases of most coarse and most precise image data.

Although set 1 selected using CNM did not provide the absolute best image quality among all optimized three-wavelengths sets, it provided PSNR values within the standard deviation from the PSNR-optimized sets. We also examined larger wavelength sets obtained by clustering [17] all wavelengths found using the optimization methods. Adding one wavelength (sets 5 and 6) resulted in about 10% improvement in terms of PSNR, which was significant given that PSNR was a logarithmic measure. No further improvement was achieved using six- wavelength sets. A “standard” set of five wavelengths (set 8) recommended by iThera Medical GmbH provided images of the quality close to those obtained using the optimized three-wavelength sets 1–4.

Three-wavelength sets close to a minimum of the modulus of the extinction matrix determinant provided lowest PSNRs and SSIMs and also the corresponding images were visually identified as the worst. It was notable that the wavelengths in “bad” combinations were not very close to each other and could be accidentally selected based on the idea to cover the range of wavelengths corresponding to the spectral features of the chromophores.

There were two sources for the distinctions between PSNR (and SSIM) values for double precision and 8-bit images: larger values for 8-bit images and larger differences between different sets for the double precision images. The first reason was the difference in the range of peak values: the standard value of 255 for the 8-bit images and the maximal pixel value of the reference image for the double-precision images. The second reason was that since 8-bit images were coarser, the overall errors were smaller. In particular, for all sets except “bad” ones (9, 10, and 11) and almost for all slices PSNR values for QC-1 images were infinite and SSIM values were equal to one, indicating that the 8-bit QC-1 images obtained using small sets were identical to the reference images (obtained using 131 wavelengths). However, this was not the case for the Hb and HbO2 images. In the case of the double precision images all PSNR values were finite and all SSIM values were smaller than one, indicating that the double precision images obtained using small sets were not identical to the reference images.

While for all sets QC-1 images showed higher PSNR and SSIM values than Hb and HbO2 images, one should note that the accurate imaging of all three chromophores was equally important for the lymphatic drainage studies because the tracer could penetrate not only the lymphatic vessels but also the blood vessels, and in addition some spectral cross-talk among all chromophores was possible.

Although the optimization of SSIM values could also be used for the wavelength selection, we found that PSNR values were much more sensitive to the differences between images than SSIM values.

Using our custom MATLAB code the PSNR optimization took about 16 h on a PC equipped with the Intel® Core™ i7 processor and 16 GB RAM. This could be done much faster using a code implemented in a faster language than Matlab, such as C, or employing a parallel computation scheme on a computer cluster. Also instead of running through all possible combinations of wavelengths, one can reduce the number of combinations by applying restrictions such as a minimum separation between wavelengths.

Conclusion

Among the three-wavelengths sets, the set selected using CNM provided the average PSNR within the standard deviation from the PSNR-optimized sets. Therefore, we are confirming that CNM method [[10], [11], [12]] can be used if one needs to limit the number of wavelengths to the number of major chromophores. Adding one wavelength resulted in about 10% improvement in terms of PSNR, which was significant given that PSNR was a logarithmic measure. No further improvement could be achieved using six-wavelength sets. Some “arbitrary” choices of wavelengths or choices based on the spectral features of chromphores (set 11) resulted in significant image distortions. The approach that we describe here may help to improve the planning of the experiments by the selection of the minimum number of wavelengths to obtain good quality images with greater temporal resolution, which was particularly helpful for the dynamic imaging of chromophores distributions required for the lymphatic drainage studies.

Conflict of interest

The authors declare that there are no conflict of interest.

Funding

None.