1. Introduction

Microvascular health is maintained by both an adequate vessel function and an optimal branching architecture, that change in cancer [1], inflammation, diabetes, heart infarct, stroke, organ transplantation, tissue regeneration and burns. The vascular bifurcation density has recently been used as a biomarker in photoacoustic (PA) images of primary breast cancer in patients [1]. Indeed, this parameter was higher in comparison to the contralateral breast. Conventional medical imaging techniques, like magnetic resonance imaging, lack the spatial resolution for the detection of this biomarker and other early functional microvascular changes that arise in the early angiogenic sprouting [1]. Photoacoustic imaging give access to them in the early onset of tumor angiogenesis at penetration depths until 5 cm in humans [2] when the excitation wavelength is adapted to the optical transparency window of a tissue. Therefore, this technique is of interest for tumor diagnosis on or near the surface and for the monitoring of tumor recurrence in the resection cavity wall, which is i.e. the case in 90% of glioma-bearing patients [3].

Hemoglobin (Hb) in RBC is an endogenous photoacoustic dye and can be excited in vivo at different wavelengths for the determination of the hemoglobin concentration and the oxygen saturation of hemoglobin [4]. Galanzha et al [5] already demonstrated an in vivo PA flow cytometry approach for the detection of RBC aggregates, but their studies were not focused on in vivo analysis on functional (micro)-vascular networks as presented here.

OR-PAM has the potential to study these processes at a microscopic resolution [6,7] at limited depth (<1 cm) and can be made handheld [8,9]. The aim of the present study is to demonstrate that 2D bifurcation maps of the vasculature in tissues can be easily generated after a simple and fast processing of OR-PAM images. For this purpose, our images were obtained with a PA microscope that uses intensity modulated laser diodes for the simultaneous excitation of Hb and other endo- or exogenous PA dyes at different wavelengths [10]. In particular, this multi-wavelength approach may be used to get rid of the signal from interfering molecules such as melanin in the skin [11].

2. Methods

2.1 The PA microscope

For a detailed description, see [10]. The intensity of a CW laser diode (λ = 415 nm) is sinus modulated at f = 5 MHz and expanded to fill the entrance pupil of the microscope objective (Zeiss LD Epiplan 20X/0.4). At the output of the objective, the beam waist is about 1 µm and the mean optical power is typically set to 2.3 mW. The light absorbed by the target is converted into acoustic waves by a thermoelastic effect. An image is obtained by moving the beams on the target with a two-axis galvanometric mirror scanner, with a pixel dwell time of 150µs. In this study, the typical imagesize was 128 × 128 pixels with an acquisition time of 2,5s. These parameters were chosen to make a compromise between acquisition time, signal-to-noise ratio and tissue preservation. We experimentally verified that no vascular damage appeared after long exposures: ∼ 20 min [10].

On the other side of the target, the PA signal is collected through an ultrasound transmission gel by a focused piezoelectric transducer in a confocal arrangement. The electric signal of the transducer is finally amplified and demodulated by a lock-in amplifier running at f to recover the amplitude of the PA signal.

2.2 Animal model

PA experiments were performed on vasculature networks in ears of white CD-1 IGS mice (Charles River, Écully, France). All efforts were made to minimize the number of mice used in the experiences. Mice were housed in ventilated cages with food and water ad libitum in a 12 h light/dark cycle at 22 ± 1°C. Before the PA microscopy, mice were anesthetized by inhalation of a gas mixture of 1% isoflurane in air/30% O₂. Their body temperature was maintained at 36–37° C using an electric heating pad with a feedback system. Ears were cleaned with a hair-removing cream (VEET, France) and rinsed with 70% ethanol, before positioning the ear on the lid of a petri dish (Ø=10 mm) in the focal plane of the objective (ZEISS LD Epiplan 20x, NA = 0,4) with ultrasound gel.

In accordance with the policy of Clinatec (permit number: B38-185 10 003) and the French legislation, experiments were done in compliance with the European Parliament and the Council Directive of September 22, 2010 (2010/63/EU) on the protection of animals used for scientific purposes. The research involving animals was authorized by the Direction Départementale des Services Vétérinaires de l’Isère – Ministère de l’Agriculture et de la Pêche, and the Ministère de l’Enseignement Supérieur et de la Recherche, France permit number: 2015051914157522_v1 (PI: Boudewijn van der Sanden, PhD, permit number 38 09 40 for animal experiences).

2.3 Image processing and analysis

2.3.1 ViTuBi analysis

For this analysis, all PA images were treated and analyzed using Image-J software: version ImageJ 1.51u, Wayne Rasband, NIH, USA, http://imagej.nih.gov/ij, Java 1.8.0_66 (64 bit). In a vessel, the pixel intensity distribution is linked to the probability that the beam intercepts a RBC, which is supposed to be frozen during the pixel dwell time [10]. Thus, the sum of a sufficient number of acquisitions represents the local RBC density or hematocrit in the vascular network during a time interval. In normal tissue, the fluctuations of the RBC density and their linear velocities are most important in the micro-vasculature. Here, in the normal microvasculature of the mouse ear, short timeseries (< 2 min) were long enough to generate a homogeneous PA signal after the sum of all images for calculations of vascular bifurcation maps. However, timeseries may become longer to generate a homogeneous PA signal in e.g. tumor microvessels that are often immature, tortuous and are irregular in diameter, which may decrease the linear velocity of RBCs.

In a next step, multiple copies of the sum of the PA-images in a time block were gathered together in a 3D stack before applying a 3D tubeness filter. The minimum number of copies (here 8) must be large enough compared to the σ parameter of the filter (see next paragraph). However, increasing this number beyond a value of 5 times σ does not yield better results. After this stacking operation, bifurcations become tubes while vessel segments between bifurcations become sheets. The first are selected and highlighted by a tubeness filter, whereas the latter are excluded by the filter.

The tubeness filter (imageJ > Plugins > Analyses > Tubeness) produces a score for how “tube-like” each point in the image is (https://www.longair.net/edinburgh/imagej/tubeness/). The plugin uses the eigenvalues of the Hessian matrix. If the larger two eigenvalues (λ₂ & λ₃) are both negative, the score is $\sqrt {{\lambda _2}{\lambda _3}} $, otherwise it is 0. Here, the Hessian matrix (H) is a result of a convolution between the signal intensities I(x,y) in a 3D stack with the second order partial derivatives of Gaussians (which is an approximate of the second partial derivative of I(x,y) [12]). The Full Width at Half Maximum (FWHM) of the Gaussian is defined by σ and should be adapted to the diameter of the vessels in an image. The σ parameter should be equal to the smallest vessel size and no vessel should be more than 5 times wider than the smallest in a single image, which is true in most cases for small regions of interest. Here, a σ value of 2 pixels and 8 copies highlighted all bifurcations, but need be adapted for each new dataset. Further, the tubeness filter may select vessel segments with strong curvatures, which is a disadvantage of the filter, see discussion.

In summary, the ViTuBi analysis for the highlighting of the bifurcations was as follows: 1) background correction with a rolling ball algorithm (radius 20 pixels) of all PA-images in time series, 2) creation of a sum image of at least 50 PA images in time, 3) ∼ 8 copies of the sum are gathered together in a stack to properly generate tubes at bifurcations, 4) using a 3D tubeness filter with σ = 2 pixels for the highlighting of bifurcation and the elimination of vessel segments. A color merge (composite image) of the total vascular network (green) and the vascular bifurcations images (red) is often generated to visually check real bifurcations.

2.3.2 Classical analysis

The ViTuBi analysis using PA Hb images was compared to classical morphological analysis of this parameter using AngioTool after step 1) and 2) in the previous paragraph [13]. Briefly, this software needs a segmentation of the (fluorescent) vessel staining and estimates of the vessel diameters, before skeletonizing (central line prediction) of the vascular network and detection of bifurcations. Segmentation is not without risks as it excludes areas of low vessel staining intensities and vessels that are closely located to each other or that are interlaced are often seen as a single structure.

3. Results

3.1 ViTuBi analysis of bifurcations in a vascular network

The micro-vascular networks are constructed in a way that they optimize oxygen and metabolite delivery to the surrounding tissue. The number of bifurcations or bifurcations/unit area determine the sprouting activity in tissues and indirectly the minimal O₂ diffusion distances between capillaries, which differ in tumors [1] or other pathologies like type 2 diabetes mellitus [14]. Measurements of the RBC distribution (or Hb PA signals) in a vascular network by PA microscopy permit to localize all bifurcations instead of using morphological analysis of vascular networks as proposed in e.g. AngioTool [13]. AngioTool software calculates, among other morphological vessel parameters, the so-called ‘‘bifuraction index’’ (bifuractions/unit area), providing a measure of the sprouting activity. In ViTuBi analysis, copies of individual Hb PA sum images are piled up in a stack and generate tubes at vascular bifurcations and sheets for vessel segments between bifurcations. In Fig. 1, the sum of all 250 Hb PA images in a time series resulted in a homogeneous PA Hb signal in the whole vascular network. When we apply a tubeness filter with σ = 2 pixels on the pile of at least 8 copies of the sum PA Hb images, the result is a bifurcation map, see Fig. 1(C). The composite image of 1B (green) & C (red), nicely highlights the bifurcations (orange dots) in the vascular network (green) without further image processing, see Fig. 1(D).

 

Fig. 1. Steps of the ViTuBi analysis. Virtual tubeness analysis of bifurcations in a vascular network. A) A representative Hb PA image (128 × 128 pixels, 2.46 s per image) in the time series. The scale bar = 50 µm. B) The sum of 250 Hb PA images in time. C) The bifurcation map is generated after applying a tubeness filter with sigma = 2 pixels on the stack of 8 copies of A). Bifurcations become tubes in the virtual 3^rd dimension, which is accentuated by the filter that simultaneously exclude the vessel segments between the bifurcations. D) is a composite image of B (green) and C (orange), which highlights the bifurcations in a vascular network (orange dots) without complex imaging processing.

Download Full Size | PDF

3.2 Comparison with a classical analysis of vascular bifurcations

In Fig. 2, our fast ViTuBi analysis of bifurcations was tested on 6 different vascular network configurations in the mouse ear (see Fig. 2(A), (C), (E) and (G)) and compared to Angiotool analyses in Fig. 2(B), (D), (F) and (H): network A shows close located and interlaced vessels, C is a representative network with a homogeneous vessel size distribution (see detail in Fig. 1), E has a heterogeneous vessel size distribution with vessels in parallel and G is a representative capillary network.

 

Fig. 2. Comparison of ViTuBi with classical analysis. Proof of principle that fast ViTuBi analysis of the RBC density (hematocrit) by OR-PAM can select bifurcations (orange dots) for different vascular networks (green) in the mouse ear with simple imaging processing in A), C), E) and G). For comparison, the number of bifurcations was classically analyzed using “Angiotool” in B), D), F) and H). I) The number of bifurcations (black squares) in 128 × 128 pixels images found with Angiotool (y-axis) was compared to the results of virtual tubeness analysis (x-axis) for 6 different mice (only 4 vascular networks are shown in A) - G)). The linear regression line (dotted black line, slope = 0.76 and r²=0.78) of the results (black squares) was forced to go through zero and the dashed red lines indicate the hypothetical regression line if the number of bifurcations would have been equal for both methods. Linear regression analyses were generated by GraphPad PRISM version 3.02, 2000 (Graphpad Software, Inc).

Download Full Size | PDF

The number of bifurcations obtained by ViTuBi analysis was compared to the results obtained by Angiotool in Fig. 2I. In general, Angiotool underestimates the number of bifurcations, especially in zones of low PA signal intensities and in zones where vessels are closely located to each other, see white dashed ROI in B), D), F) and H).

The ViTuBi analyses tend to slightly overestimate the number of bifurcations, because vessel sections that are perpendicular to the imaging plane or that have strong curvatures, irregular widths or intensity variations may produce a significant tubeness score: see the white dashed rectangular ROIs in 2A), C), E), G) and Fig. 3, which is a detail of Fig. 2(E). Perpendicular vessels were not observed in our experiments, because the vascular network of the mouse ear is very little developed in the third dimension. The tubeness scores of the vessels irregularities which are not bifuractions are often lower than in real bifurcations, but they are hard to separate by intensity thresholding, see Fig. 3 and discussion. In Fig. 2(A)): 5% of the orange dots are not bifurcations. The F1-score of the ViTuBi analysis is 0.97 with the precision (p) = 0.95 and recall (r) = 1. In Fig. 2(C)) 7% of the orange dots are not bifurcations, the F1-score is 0.96 with p = 0.93, r = 1. In Fig. 2(E) 6% of the orange dots are not bifurcations, the F1-score is 0.97 with p = 0,94, r = 1 and in Fig. 2(G) only 2% of the orange dots are not bifurcations, the F1-score is 0.99 with p = 0.98, r = 1. The mean F1 score of all 6 networks is 0.97, which indicate that the accuracy of the ViTuBi analysis is high.

 

Fig. 3. Minor disadvantages of ViTuBi analysis with a tubeness filter. 3A) Detail of Fig. 2(E) showing curvatures in the white dashed rectangular ROIs that are selected by the tubeness filter (σ = 2 pixels) in the ViTuBi analysis. 3B) These curvatures: C1–C3 have comparable intensities as the bifurcations: B1-B3 after filtering and therefore, they cannot be separated be simply intensity thresholding, see 3C. However, they can be distinguished by their elongated shape. 3C) In the surface plot of 3B) e.g., the intensities of C1 and B1-B3 are comparable.

Download Full Size | PDF

4. Discussion

In the previous sections, OR-PAM is validated as a powerful tool for analyzing parameters of micro-vascular networks. The implementation of the equipment and the image processing and analysis are straightforward and adapted for clinical applications.

The current configuration of the PA microscope with the piezo-electric transducer in the transmission mode should be converted in an epi-collection mode for a hand held configuration [8]. Both the objective and the transducer will use a selective dichroic mirror for the transmission of light and reflection of ultrasound waves. This facilitates the alignment of both the excitation with light and detection of the generated ultrasound waves in the focal plane. The use of pulsed laser diodes should be compared to the current frequency encoded laser diode excitation [9], as well as the sensitivity of CMUT arrays for PA detection in stead of piezoelectric transducers with a limited sensitive area. Laser diodes are more compact and less expensive than the pulsed lasers that are classically used in PA imaging for a better penetration in light scattering tissues. The absorption coefficient of HbO₂ is maximum at 415 nm, but a better compromise between efficient light absorption and penetration depth in tissues will be a shift to longer wavelengths ≥ 567 nm.

In this paper, we focus on ViTuBi analysis of the bifurcation density, which is an important hallmark for the sprouting activity in tumor angiogenesis [1] and other pathologies such as: in age-related macular degeneration and diabetic rethinopathy [14–16]. In previous clinical studies on breast cancer patients [1], the density of vascular bifurcations was used as a biomarker to distinguish between normal breast tissue and tumor angiogenesis with a higher bifurcation density. Here, manual counting of the vascular bifurcation density in 2D photoacoustic images of patients could be easily replaced by our ViTuBi analysis of vascular bifurcations after a short Hb PA experiment in time.

In a comparison between ViTuBi - and classical analysis with Angiotool [13], we observed that the first method may slightly overestimates the number of bifurcations and the second method often underestimates this parameter. Nevertheless, the precision of the ViTuBi analysis is high: the overall F1 score of all 6 networks in Fig. 2 is 0.97. The analysis does not generate false negative results (see r-values), but may include false positives, because a 3D tubeness filter may select irregularities such as strong curvatures (C1-C3) in vessels, see white dashed ROIs in Fig. 3. These curvatures have comparable intensities as the bifurcations, e.g.: compare C1 to intensities of B1-B3 in Fig. 3. Thus, a simple intensity thresholding is not possible, but curvatures can easily be distinguished from bifurcations, because they have elongated shapes in comparison to bifurcations that are large dots. Circularity image analysis of the latter will exclude the low percentage of elongated shapes.

In classical analysis, a user interaction is required for different steps: thresholding of fluorescence images, and estimation of the vessel diameters for the correct estimation of skeletal vessel segments and bifurcations. Region with low fluorescence intensities may be missed and several vessels that are located close to each other or that are interlaced are detected as a single vessel. Therefore, classical analysis tends to underestimate the number of bifurcations. Both methods are not perfect, but ViTuBi analysis are faster and only require a user interaction for the setting of the σ value of the tubeness filter.

In conclusion, fast and robust ViTuBi analysis of bifurcations in a vascular network is possible with OR-PAM using the endogenous Hb PA signals in short time series. When at least 8 copies of 2D PA Hb sum images of the series are piled up in a stack, vascular bifurcations become tubes and can be highlighted by the 3D tubeness filter. This method requires fewer interactions with the user in comparison to morphological analysis as proposed by e.g. Angiotool. This method often underestimates the number of bifurcations, especially in tiny (tumor) capillary networks. Indeed, tumor angiogenesis includes in the early stage more tiny vessels with low PA signal intensities than in a mature vessel network; therefore, the ViTuBi analysis may be more robust than classical approaches. This will be tested in a next study on small tumors (< 5 mm diameter) with an early stage of angiogenesis. These are representative for early recurrence or metastatic outgrow.