In vivo imaging of blood vessels and lymphatic vessels is important for monitoring physiological parameters and early diagnosing of fatal diseases.^1–3 The blood and lymphatic vessels are two sets of circulatory networks that have distinct but interdependent functions. Lymphatic vessels are found in most tissues and play an important role in many physiological processes.⁴ The blood circulation supplies oxygen and nutrients, removes waste, and is highly associated with the lymphatic circulation. Recent research has shown that blood flow, local oxygen content, and lymphatic vessels are important factors in tumor growth and metastasis.^5–14 Simultaneous imaging of blood and lymphatic vessels is of great value in disease diagnosis and treatment.³ Most existing bioimaging techniques are limited to either low spatial/temporal resolution or unable to simultaneously image multi-contrasts in the blood and lymphatic vessels. A high-resolution multi-contrast imaging technique that can image both the blood and the lymphatic vessels is highly needed.

Recently, optical-resolution photoacoustic microscopy (OR-PAM) has been developed for in vivo imaging of hemoglobin concentration (CHb), oxygen saturation (sO2), blood flow speed (v), metabolic rate of oxygen, and lymphatic vessels.^15–43 Although promising, limited by the laser source, simultaneous imaging of blood and lymphatic vessels has not been implemented in OR-PAM. In a single scanning, OR-PAM can only quantify one or two contrasts. Repetitive scanning must be done to map multiple functional or molecular contrasts, which lowers the throughput and temporal correlation among different quantities.

Here, we present five-wavelength OR-PAM that can image the blood and lymphatic vessels in a single scanning. Multiple contrasts, including hemoglobin concentration, blood flow, oxygen saturation, and exogenous dyes, can be quantified. Based on the stimulated Raman scattering (SRS) effect,⁴⁴ we develop a five-wavelength (532, 545, 558, 570, and 620/640nm) nanosecond pulsed laser. At each scanning spot, five A-lines at these wavelengths are sequentially acquired. The 532-nm wavelength is used to image hemoglobin concentration; the 532- and 545-nm wavelengths are used to measure the blood flow speed using a dual-pulse method; the 532-, 545-, 558-, and 570-nm wavelengths are used to determine the oxygen saturation (sO2); and the 620/640-nm wavelength is used to image the dye-labeled lymphatic vessels. Using the five-wavelength OR-PAM system, we demonstrate high-resolution multi-contrast in vivo imaging of blood and lymphatic vessels in tumor growth, lymphatic dye clearance, and brain function monitoring.

The five-wavelength pulsed laser needs to have sufficient pulse energy for each wavelength and nanosecond wavelength switching time. Figure 1 shows a schematic of the laser and the imaging probe. We use a pulsed laser (532-nm wavelength, VPFL-G-20, Spectra-Physics) to pump the other wavelengths via the SRS effect in optical fibers. The Raman frequency shift of the silica fiber is ∼13.2THz, corresponding to a 13-nm wavelength shift near 532 nm. The five wavelengths are 532, 545, 558, 570, and 620/640nm. The five pulses are temporally separated with different optical delays in the fiber. The pulse energy on the sample surface for each wavelength is above 90 nJ.

In the SRS shifter, two polarizing beam splitters (PBS1 and PBS2, PBS251, Thorlabs Inc.) and one 10/90 beam splitter (BS, BSN10, Thorlabs Inc.) split the pump laser beam into four paths, i.e., a 532-nm direct path, a 545-nm Raman path, a 620/640-nm Raman path, and a 558/570-nm Raman path. The pump pulse width is 7 ns. The total pump energy is 23.5μJ per pulse. The incident pulse energies for the split paths are 0.25μJ for the 532-nm path, 2.25μJ for the 545-nm path, 8μJ for the 620/640-nm path, and 13μJ for the 558/570-nm path. Half-wave plates (HWP1−6, WPH10E-532, Thorlabs Inc.) are used to control the polarization so that the SRS efficiency can be adjusted. Neutral density filters (NDF1−5, NDC-50C-2, Thorlabs Inc.) are used to adjust the laser pulse energies.

The 532-nm path is coupled into the OR-PAM probe with minimal delay. The 545-nm path uses a 25-m polarization-maintaining single-mode fiber (PM-S405-XP, Nufern) to shift the wavelength and delay the pulse by 121 ns. A long-pass filter (LPF, RET537lp, Chroma) is used to select the 545-nm wavelength. In the 620/640-nm Raman path, a 50-m fiber (PM-S405-XP, Nufern) is used to generate the 620/640-nm wavelength and delay it by 243 ns. In the 558/570-nm Raman path, the pump beam is coupled into a 100-m graded-index multimode fiber (MMF) [GIMMSC (50/125)HT, Fibercore] to generate the 558- and 570-nm wavelengths and delay them by 486 ns. Because the multimode fiber has a larger core diameter (50μm), it has a higher SRS threshold and thus can output higher pulse energy than a single-mode fiber at the same length.⁴⁵ Furthermore, the MMF can generate the SRS wavelengths with at least 75% output energy in the fundamental mode.⁴⁶ A long-pass dichroic mirror (DM1, T565lpxr-UF1, Chroma) is used to separate the 558- and 570-nm wavelengths. Another 25-m fiber (PM-S405-XP, Nufern) is used to delay the 570-nm pulse by 121 ns. A short-pass filter (SPF, 580SP, Omega) is placed at the output end of the 25-m fiber to select the 570 nm wavelength.

The five wavelengths are combined using four beam combiners. First, a long-pass dichroic mirror (DM2, T565lpxr-UF1, Chroma) is used to combine the 558- and 570-nm paths. Second, the 558/570-nm beam is combined with the 620/640-nm beam via a short-pass dichroic mirror (DM3, #69-204, Edmund Optics). Third, a long-pass dichroic mirror (DM4, T550lpxr, Chroma) combines the 545-nm path and the 558/570/620/640-nm path. Last, the four SRS wavelengths are combined with the direct 532-nm path via a 10/90 BS and are coupled into the OR-PAM probe via a 2-m single-mode fiber (P1-460B-FC-2, Thorlabs Inc.).

The coupling efficiencies for all single-mode fibers are above 50%. Because the 558-nm beam has high-order modes, when coupling it to the 2-m single-mode fiber, the coupling efficiency is ∼30%. In the 532-, 545-, 558-, and 570-nm light paths, the pulse energies on the sample surface are above 90 nJ. The pulse energy for the 620/640-nm wavelength can reach 300 nJ on the sample surface. Detailed pulse energies and time delays for all wavelengths are summarized in Table 1. The shortest wavelength switching time is 121 ns, which corresponds to ∼200-μm imaging depth. This is acceptable for superficial imaging. For future deeper in vivo application, we may consider separating mixed signals from different wavelengths. For example, a deconvolution method has been developed to separate two overlapping signals generated from ultrafast dual-wavelength excitation.⁴⁷

In the OR-PAM probe, the laser beam from the 2-m fiber is focused by a pair of achromatic doublets (AC064-013-A, Thorlabs Inc.). The focused optical beam is reflected on an optical/acoustic beam combiner, is transmitted through a planoconcave lens (45-697, Edmund Optics), then illuminates the sample. Induced ultrasonic waves are collimated by the planoconcave lens, transmitted through the optical/acoustic beam combiner, and detected by a 50-MHz broadband piezoelectric transducer (V214-BC-RM, Olympus). To optimize the detection sensitivity, the focused optical beam is coaxially and confocally aligned with the focused ultrasonic detection beam. The detailed information about the OR-PAM probe can be referred to in previous publications.^46,48–50

Fiber coupling and SRS efficiency may be affected by airflow and temperature changes, causing random errors and long-term drift in the pulse energy. To improve the pulse-energy stability, we isolate the airflow using a cover and maintain the temperature fluctuation within ±0.1°C.

As shown in Table 1, we measured the pulse energy fluctuation and drift at each wavelength for one hour. The fluctuation is calculated from the standard deviation of the pulse energy. The drift is determined from the averaged pulse energy difference between the first and the last minute in the hour. With airflow isolation and temperature control, the pulse energy fluctuations at all wavelengths are reduced by 4% to 20%, and the drifts for all wavelengths are reduced by 11% to 23%. Because the pulse energy ratios among 532, 545, 558, and 570 nm are used in the calculation of blood flow and sO2, we list the relative pulse-energy fluctuations, i.e., the standard deviations of the pulse energy ratios of 545/532, 558/532, and 570/532, in Table 1. The relative pulse energy fluctuations of 545, 558, and 570 nm are consistently smaller than their absolute pulse fluctuations. Because the MMF has a large mode area and is more sensitive to perturbation, the pulse energy fluctuations and drifts are larger in the multimode fiber (558 and 570 nm) than in single-mode fibers (545 and 620/640nm). With the airflow cover and temperature control, the pulse energy fluctuations of the SRS wavelengths (relative for 545, 558, and 570 nm and absolute for 620/640nm) are 26% to 41% higher than the pump wavelength, and the drifts are 20% to 75% higher than the pump one.

We measured the imaging depths for each wavelength by obliquely inserting a black human hair into fresh chicken breast tissue. The imaging depths are 0.75, 0.76, 0.77, 0.78, and 0.81 mm for 532, 545, 558, 570, and 620/640nm, respectively, as shown in the Figs. 2(a)–2(e). The lateral resolutions for the five wavelengths are measured by scanning the stainless-steel blade and measuring the amplitude profiles across the sharp edges. By fitting the PA signals to an edge spread function (ESF), we can obtain the line spread functions (LSFs) and the full widths at half maximum (FWHMs) as the lateral resolutions, as shown in Figs. 2(f)–2(j). The lateral resolutions are 3.12, 3.25, 3.37, 3.51, and 3.94μm for 532, 545, 558, 570, and 620/640nm, which agree with the simulation results shown as Fig. S1 in the Supplementary Material. We also measured the axial resolutions by scanning a 10-μm-diameter tungsten filament. The FWHMs of a Hilbert-transformed A-line are shown in Fig. 2(k). The axial resolutions for all wavelengths are 38μm.

Figure 3(a) shows the laser spectrum measured with an optical spectrometer (USB 2000+, Ocean Optics). The linewidths of the SRS wavelengths are within ±1.5nm.^46,48 We use the molar extinction coefficients at the peak wavelengths for sO2 calculation. Figure 3(b) shows the absorption spectra of oxyhemoglobin, deoxyhemoglobin, and Evans blue (EB).⁵¹ The hemoglobin molecules and the EB dye have distinct absorption coefficients at 620 to 640 nm. With EB dye uptake, the lymphatic vessels can be imaged with the 620/640-nm wavelength. The other wavelengths are used to image the CHb, sO2, and flow speed in the blood vessels.

Under thermal and stress confinements, the PA amplitude in the blood vessel can be approximated as $$\begin{gathered} P=kΓηFμa, \\ μa=rCHb[sO2εoxy+(1−sO2)εde], \end{gathered}$$(1)where k is a constant related to the PA detection sensitivity, Γ is the local Grüneisen coefficient, η is the efficiency of the absorbed optical energy that converts into heat, F is the local optical fluence, μa is the absorption coefficient, r is a constant,¹CHb is the total hemoglobin concentration, and εoxy and εde are the molar extinction coefficients of HbO2 and HbR at the excitation wavelength, respectively.

At the 532-nm wavelength, εoxy and εde are almost the same, thus μa is not sensitive to the sO2 change.⁵² We can determine the relative change of CHb from the PA amplitude at 532 nm or the absolute CHb after calibration.

The 545- and 532-nm wavelengths are used to measure the blood flow speed using a dual-pulse photoacoustic flowmetric method that is published in Ref. 46. When we excite two PA signals with a sub-microseconds delay, the residual heat from the first pulsed excitation elevates the second PA amplitude via the Grüneisen relaxation effect.⁵³ Blood flow may accelerate the heat dissipation and thus alter the second PA amplitude as follows: P2=r1F2F1P1+r1F2F1P12Ae−(τα+bv)δt,(2)where r1=μa2μa1, P1,2 and F1,2 are the PA amplitudes and the optical fluences at λ1=532nm and λ2=545nm, δt is the time delay between the 532- and 545-nm laser pulses, and A, τα, and b can be determined via calibration.⁴⁶ The detailed calibration method is shown as part 2 in the Supplementary Material. Here, A, τα, and b were calibrated as 0.09±0.04 (SD), 9.68±0.85 (SD), and 0.25±0.06 (SD). Via measuring the two PA signals at the 532- and 545-nm wavelengths (two isosbestic points), we can determine the blood flow speed using only two A-lines. After measuring the vessel diameter and the average flow speed, we can further determine the blood flow rate.^54–56 The detailed calculation method is shown as Eq. S(1) in the Supplementary Material.

Using 532-, 545-, 558-, and 570-nm wavelengths, we can measure oxygen saturation. Considering the Grüneisen relaxation effect among successive pulse excitations, a PA amplitude is affected by previous excitations. Because the 620/640-nm PA signal (the third pulse) from the blood vessel is only ∼2.5% of those by other pulses, we neglect the heating by the third pulse. The PA amplitudes for the first, second, fourth, and fifth pulses can be written as $$\begin{gathered} P1=kΓ0ηF1μa1,P2=kΓ0ηF2μa2[1+AP1e−(τα+bv)δt12], \\ P4=kΓ0ηF4μa4[1+AP1e−(τα+bv)δt14+AP1F2μa2F1μa1e−(τα+bv)δt24], \\ P5=kΓ0ηF5μa5[1+AP1e−(τα+bv)δt15+AP1F2μa2F1μa1e−(τα+bv)δt25+AP1F4μa4F1μa1e−(τα+bv)δt45], \end{gathered}$$(3)where δtij is the time interval between the i’th and j’th pulses. With the calibrated A, τα, and b, and determined flow speed v, we can solve the above equations and determine the sO2 value from the four measurements.

Four-week-old female ICR mice were used in experiments. All procedures involving animals were approved by the animal ethical committee of the City University of Hong Kong. The mice were anesthetized by inhaling isoflurane gas. EB dye solution (0.5%, mass fraction, 20μL) was injected into the ear tip to label the lymphatic vessels. The mouse ear was placed on a flat holder and was imaged in 5 min after EB injection. The laser repetition rate is 4 kHz for each wavelength. The pulse energy is 160 nJ for the 620/640-nm wavelength and 90 nJ for the other four wavelengths. The scanning step size is 2.5μm. The imaging area is 2.5mm×2.5mm. It takes 250 s to acquire a volumetric image. Maximum-amplitude-projected images are shown in Fig. 4. Figure 4(a) is the sO2 image. The sO2 images without and with compensation for the Grüneisen relaxation effect are shown as Fig. S2 in the Supplementary Material. The average sO2 values calculated with and without compensation for the Grüneisen relaxation effect are significantly different (P<0.001). With compensation, the sO2 values in the arteries and veins are closer to their normal physiological values.^19,20,23,46 The compensation for the Grüneisen relaxation improves the averaged sO2 values by 4% to 17% in the arteries and by 21% to 24% in the veins. The sO2 in the veins is corrected more than that in the arteries. The possible reason is that the vein has slower flow speed and thus has a more severe Grüneisen relaxation effect than the arteries. Due to the same reason, the sO2 is corrected more in the small arteries than in the big ones. Because the flow speed and the Grüneisen relaxation effect in the veins are relatively uniform, the sO2 correction shows no obvious difference in the big and small venous vessels. Figure 4(b) shows the blood flow speed image. Figure 4(d) depicts the average sO2 and flow speeds in different vessel segments. From the root to the tip of the mouse ear, the averaged blood flow speed in the trunk arteries decreases from 5.8±0.24 (SD) to 5.0±0.21 (SD) mm/s. In the small arterial branches, the averaged blood flow speed decreases to 3.4±0.36 (SD) mm/s. The averaged blood flow speed in the veins maintains at ∼1.9±0.18mm/s. After measuring the blood vessel diameter, we can further compute the flow rate. In Fig. 4(b), the profiles across five arteries and six veins are labeled with white dashed lines as A1−5, and V1−6, respectively. As shown in Fig. 4(e), the vessel diameters are 86, 74, 58, 47, and 43μm for A1−5 and 189, 112, 97, 84, 48, and 123μm for V1−6. The corresponding flow rates are 1.03, 0.65, 0.41, 0.25, and 0.15μL/min for A1−5 and 1.68, 0.63, 0.42, 0.28, 0.14, and 0.68μL/min for V1−6. The veins V2, V3, and V6 are converged to V1, and the total flow rate in V2, V3, and V6 matches the flow rate in V1. The total flow rate in V4 and V5 is approximately equal to the one in V3. Similarly, the arteries A2 and A3 are the branches of A1, and thus the total flow rate in A2 and A3 matches the one in A1. The total flow rate in A4 and A5 matches the flow rate in A3. Four pairs of arteries and veins (A2−5 and V2−5) are compared. The flow rates in the arteries match the ones in the veins. The flow rates in A1 and V1 are not very close. We guess the possible reason might be that the two vessels are different in shape, and the measured diameter cannot be used to accurately estimate the cross-sectional area. Figures 4(c) and 4(d) show the blood and lymphatic vessels and the relative lymphatic concentration in vessel segments. With distinct absorption spectra, the blood and the dye-labeled lymphatic vessels can be identified with high contrast. The highest SNRs (signal amplitude over noise standard deviation) of the blood and lymphatic vessels are 20.1 and 22.3 dB.

To demonstrate the potential applications, we first conducted in vivo early cancer detection experiments using five-wavelength OR-PAM. Before tumor implantation, the multifunctional vascular images of hemoglobin concentration, compensated oxygen saturation, and blood flow speed were acquired in the mouse ear and are shown as Fig. S3 in the Supplementary Material. Then, 4T1 breast cancer cells were injected to the mouse ear. After the tumor developed for 5 days, the same region was imaged to demonstrate the changes of the aforementioned multi-functional parameters. As a hallmark of cancer,^57,58 obvious angiogenesis is observed in Fig. 5(a). Figures 5(b) and 5(c) show the sO2 and blood flow speed in the tumor region. Compared with the baseline, the averaged sO2 in the tumor vessels increases by 45% to 51%, and the blood flow speed increases by 32% to 37%, indicating hyperoxia and elevated blood perfusion in the early-stage tumor due to fast growing.^55,57,58Figure 5(d) shows pseudo-color-encoded vascular depth, manifesting height changes in the tumor region. Figures 5(e) and 5(f) show the vessel dilation and vessel tortuosity in the tumor region. Compared with the non-tumor region, the average diameter and tortuosity of the venules are increased by 48% to 52% and 10% to 12% in the tumor region, which has been validated in other reports.^57,58 Detailed quantification methods and comparisons are shown as Fig. S4 in the Supplementary Material.

Enabled by the new wavelength, five-wavelength OR-PAM can image the dye-removal processes in lymphatic vessels. Figures 5(g)–5(i) show the simultaneous imaging of the blood and the dye-labeled lymphatic vessels in the mouse ear at 0, 10, and 20 min after injecting EB dye solution (0.5%, mass fraction, 20μL). The normalized concentrations of hemoglobin and EB are encoded with red and blue colors, respectively. After 20 min, the average EB concentration decreases by 87% in the lymphatic vessels, indicating the waste clearance ability of the local lymphatic vessels. Taking advantage of simultaneous multi-contrast imaging, five-wavelength OR-PAM can realize the simultaneous in vivo dye-labeled lymphatic and multi-functional vascular imaging, which is of great importance to imaging the dynamic blood perfusion and lymphatic circulation.

We also demonstrated five-wavelength OR-PAM in multi-contrast brain imaging. A 2mm×2mm window was opened in the mouse skull after the mouse was anesthetized. Using 60-nJ pulse energy, we acquired functional brain images and quantified normalized hemoglobin concentration, sO2, and blood flow speed, as shown in Figs. 5(j)–5(l). In the arteries and veins, the average sO2 values are 0.95±0.06 (SD) and 0.68±0.08 (SD); the average blood flow speeds are 7.2±0.62 (SD) mm/s and 3.9±0.45 (SD) mm/s, respectively. Because the new five-wavelength system has red light, other exogenous photoacoustic contrast agents can be used in the future to label neurons or tumors in the brain.

We developed five-wavelength OR-PAM that can simultaneously image the EB-labeled lymphatic vessel, as well as hemodynamic and functional information in the blood vessel in a single scanning. A five-wavelength pulsed laser was developed based on the SRS effect. Fiber-based optical delay switches the wavelengths within sub-microseconds, which reduces the misalignment in multi-wavelength imaging. We isolate airflow and control the temperature of the SRS laser to further stabilize the pulse energy. The 532-, 545-, 558-, and 570-nm wavelengths are suitable for blood vessel imaging, and the 620/640-nm wavelength offers high contrast for lymphatic vessel imaging. With the blood flow speed measured by the dual-pulse flowmetric method and corresponding vessel diameter, the blood flow rates in several arteries and veins were estimated, and the results manifest the flow conservation. In vivo experiments demonstrated the simultaneous multi-contrast imaging of blood and lymphatic vessels. We showed the potential applications of this technique in early cancer detection, monitoring lymphatic circulation, and imaging the brain. Other molecular or nano-probes that have strong absorption at 620/640nm can also be used as the contrast agent in the future. Co-registered hemodynamic, functional, and molecular information can be acquired in a single scanning. This technical advance enables simultaneous multi-contrast imaging in the living body, offering a new biomedical imaging tool.