Quantitative Characterization of Age-Related Changes in Peripheral Vessels of a Human Palm Using Raster-Scan Optoacoustic Angiography

: The analysis of age-related changes in skin vessels based on optoacoustic angiographic images during the in vivo skin monitoring of healthy volunteers at different ages is reported. As a result of a quantitative analysis of the three-dimensional OA images, the age-associated differences in the following image parameters were revealed: image intensity, ratio of blood content at different characteristics depths, total vessel length, and number of branches. The reported approach can be effectively employed for automatic assessment and monitoring of age-related vascular changes in the skin and underlying tissues.


Introduction
Age-related changes in peripheral vessels are manifested by a number of structural and functional disorders that reduce the blood supply to tissues, thereby weakening the processes of their gas exchange and hindering the intake of nutrients and the outflow of metabolic products. To date, a number of primary criteria for age-related changes in skin vessels have been identified, which include a decrease in vascular density, violation of branching geometry, a decrease in the number of collaterals, and suppression of vascular reactivity [1]. Similar changes in the peripheral vascular bed can accompany the development of angiopathies common in the population, being a predictor of cardiovascular death [2,3]. Thus, the geometric features of the vessels of different organs in diabetes mellitus are described, and the prognostic value for the development of the disease is shown [4,5]. In atherosclerotic lesions of the arteries, developed collaterals are able to compensate and improve the condition, and a decrease in capillary density is considered a negative prognostic sign [6,7]. Chronic venous diseases with the formation of venous insufficiency are associated with a violation of the architectonics of peripheral vessels, giving positive dynamics in response to adequate therapy [8]. It has been proven that various stages of Raynaud's syndrome, affecting up to 10% of the population [9], is characterized by vasodilation or vasospasm and decreased circulation in the periphery [10,11].
A combination of natural age-related and other pathological changes can mask the true picture of the disease, leading to a late diagnosis. At the same time, vascular dis-eases significantly increase with the aging of the population [12][13][14][15], which also makes it necessary to differentiate the natural age-related changes in angioarchitectonics.
Currently, age-related changes in the microvasculature of surface tissues are studied using both ex vivo and in vivo methods. The method of immunohistochemical analysis with markers of vascular endothelium showed a gradual decrease in the vascularization of human skin with age [16][17][18]. Decreased capillary loop density and increased vessel length were demonstrated using capillaroscopy and Laser Doppler flowmetry [19,20]. A promising modern method for solving the problem of studying age-related changes in the vascular system of the skin is angiographic optical coherence tomography. This method allowed for the noninvasive identification of age-related decreases in the vessel density and diameter of skin vessels [21], but the visualization depth was up to 350 µm.
Raster-scan optoacoustic (OA) angiography is method being actively developed that is often used to visualize blood vessels, providing optical [22] or acoustic [23] resolution of tens of microns at the imaging depths of several millimeters, which makes it possible to obtain in vivo morphological and functional information on the microvessels' state. However, due to the large amount of 3D data acquired by modern high throughput raster scanning OA systems, the resulting OA vasculature images can be difficult to analyze visually. Therefore, the development of algorithms that extract quantitative information about vessels [24] from OA data is important for analyzing the tissue state.
In this paper, age-related changes in the vascular bed of the skin of healthy volunteers are assessed in vivo for the first time using the OA angiographic technique. As a criterion for analyzing age-related changes in the skin, various parameters of the vascular network of the skin were calculated using a recently developed [25] and successfully applied in [26] algorithm for processing the three-dimensional OA images based on graph representation of the vascular net. The proposed algorithm was used to calculate the parameters characterizing the average density of the vessels, the length of the vessels, and the branching of the vessels on the basis of a graph constructed from a skeletonized image.

In Vivo Experiments
In vivo studies were conducted with healthy volunteers 19 to 74 years of age with Caucasian skin (Fitzpatrick type II-III) without confirmed angiopathies, neuropathies, and injuries of the studied areas. The volunteers received no medications during the 3 days prior to the measurement. All studies were carried out under the same conditions; before the start of the investigation, the volunteer was kept the horizontal position for 10 min with a horizontally-located limb. A noninvasive assessment of age-related vascular features of the skin and underlying tissues from the thenar area of both palms was performed. Initially, 8 measurements were taken; however, the images with artifacts (affected by a person's breath, movements, etc.) were excluded by manual selection at the stage of dataset preparation. As a result, for each volunteer, from 2 to 8 measurements were included providing a total of 128 3D images of thenar vasculature. To analyze age-related changes, the volunteers (both female and male) were divided into three age groups: group 1: <35 years (n = 10), group 2: 35-49 years (n = 10), and group 3: ≥50 years (n = 10). All in vivo experiments were approved by the Ethical Committee of Privolzhsky Research Medical University (Protocol #17 from 11 October 2019).

Raster-Scan Optoacoustic Angiography System
Experimental data were obtained using a custom-made bright-field OA probe (Figure 1) consisting of a polyvinylidene difluoride [27] ultrasonic (US) detector (focal distance F = 6 mm, numerical aperture NA = 0.6, and bandwidth 1-40 MHz) and a multimode optical fiber (diameter 400 µm, NA = 0.22). The scanning OA probe was integrated into an OA microscopy system [28] consisting of a 532 nm laser (pulse repetition rate of 2 kHz, a pulse duration of 1 ns, and a pulse energy of 0.2 mJ) and XY scanning stages (scanning area Photonics 2022, 9, 482 3 of 10 of 5 × 5 mm 2 , scanning step 25 µm). The acquisition time for one three-dimensional image was~5 min.
Photonics 2022, 9, x FOR PEER REVIEW 3 of 10 OA microscopy system [28] consisting of a 532 nm laser (pulse repetition rate of 2 kHz, a pulse duration of 1 ns, and a pulse energy of 0.2 mJ) and XY scanning stages (scanning area of 5 × 5 mm 2 , scanning step 25 µm). The acquisition time for one three-dimensional image was ~5 min.

Quantitative Algorithm for Assessing Vascular Changes
The first stage of the OA image analysis included preprocessing. Two-dimensional delay-and-sum image reconstruction in the frequency domain [29] was applied to the raw OA B-scans in two planes XZ and YZ. The reconstructed 3D images were further Frangifiltered [30] using the "vesselfilter" function from the k-Wave software package [31]. In the second stage, the quantitative algorithm described in [25] was employed. A mask with a manually selected threshold (the same for all images) was constructed from the Frangifiltered OA images, converted to binary form, and then skeletonized. A three-dimensional binary voxel skeleton was transformed into a graph [32], characterized by vertices and edges, consisting of "vascular trees" (connected graph structures). The algorithm [25] enabled extraction of the following parameters: the total number of nodes per unit volume, the total number of branches per unit volume, the vessel amount per unit volume (further referred to as "vessel density"), the total length of all vessels per unit volume, rib average length (calculated as the mean distance between two branching points in a blood vessel), and the average length of one vessel. In addition, for a quantitative comparison of the OA images of different age groups, the total intensity over the entire image (further referred to as "image intensity") and the total intensity over the vascular mask (further referred to as "vessel intensity") were calculated. In addition, for each image, the blood content indepth profile was calculated from the surface to the depth of 2 mm as a ratio of the amount of voxels corresponding to blood vessels at a particular depth to the total number of voxels at this depth.

Quantitative Algorithm for Assessing Vascular Changes
The obtained quantitative parameters were presented as the average values of each parameter calculated from the OA images and their standard errors of the mean (SEM). The statistically significant differences between age groups were assessed using one-way ANOVA with a posteriori analysis of variance followed by the Duncan's test. Differences

Quantitative Algorithm for Assessing Vascular Changes
The first stage of the OA image analysis included preprocessing. Two-dimensional delay-and-sum image reconstruction in the frequency domain [29] was applied to the raw OA B-scans in two planes XZ and YZ. The reconstructed 3D images were further Frangifiltered [30] using the "vesselfilter" function from the k-Wave software package [31]. In the second stage, the quantitative algorithm described in [25] was employed. A mask with a manually selected threshold (the same for all images) was constructed from the Frangifiltered OA images, converted to binary form, and then skeletonized. A three-dimensional binary voxel skeleton was transformed into a graph [32], characterized by vertices and edges, consisting of "vascular trees" (connected graph structures). The algorithm [25] enabled extraction of the following parameters: the total number of nodes per unit volume, the total number of branches per unit volume, the vessel amount per unit volume (further referred to as "vessel density"), the total length of all vessels per unit volume, rib average length (calculated as the mean distance between two branching points in a blood vessel), and the average length of one vessel. In addition, for a quantitative comparison of the OA images of different age groups, the total intensity over the entire image (further referred to as "image intensity") and the total intensity over the vascular mask (further referred to as "vessel intensity") were calculated. In addition, for each image, the blood content in-depth profile was calculated from the surface to the depth of 2 mm as a ratio of the amount of voxels corresponding to blood vessels at a particular depth to the total number of voxels at this depth.

Quantitative Algorithm for Assessing Vascular Changes
The obtained quantitative parameters were presented as the average values of each parameter calculated from the OA images and their standard errors of the mean (SEM). The statistically significant differences between age groups were assessed using one-way ANOVA with a posteriori analysis of variance followed by the Duncan's test. Differences were considered statistically significant at p < 0.05. Statistical data processing was carried out in Excel and Statistica 12.0 programs.  were considered statistically significant at р < 0.05. Statistical data processing was carried out in Excel and Statistica 12.0 programs.  In accordance with the modern anatomy paradigm, normal skin microcirculation is organized in two plexuses parallel to the skin's surface. The upper plexus, consisting of small vessels, is located near the papillary dermis. The lower plexus is located at the dermal-hypodermal interface and consists of larger arteries and veins [33].  Figure 2i was added for visualization of this separation) for a volunteer from group 3, while for the younger volunteers this separation was not manifested. (Figure 2g,h). In Figure 3a the average blood content in-depth profiles for In accordance with the modern anatomy paradigm, normal skin microcirculation is organized in two plexuses parallel to the skin's surface. The upper plexus, consisting of small vessels, is located near the papillary dermis. The lower plexus is located at the dermal-hypodermal interface and consists of larger arteries and veins [33].  Figure 2i was added for visualization of this separation) for a volunteer from group 3, while for the younger volunteers this separation was not manifested. (Figure 2g,h). In Figure 3a the average blood content in-depth profiles for all three groups are shown together with corresponding SEM ranges. In the superficial layers, the blood content demonstrated a high correlation with age group: group 1 demonstrated low blood content at typical depths up to 300 µm as compared to the lower plexus, while group 2 demonstrated a comparable value of blood content for both the superficial layers and lower plexus. Group 3 demonstrated two distinguished blood layers separated at a depth of about 300 microns with a pronounced dip of blood content in this area (Figure 3a, red curve). The ratio of blood content at the two characteristic depths z 1 = 120 µm and z 2 = 570 µm, corresponding to the local maxima of red curve in Figure 3a, demonstrated statistically significant differences between the three age groups and increased with age (Figure 3b). Photonics 2022, 9, x FOR PEER REVIEW 5 of 10 all three groups are shown together with corresponding SEM ranges. In the superficial layers, the blood content demonstrated a high correlation with age group: group 1 demonstrated low blood content at typical depths up to 300 µm as compared to the lower plexus, while group 2 demonstrated a comparable value of blood content for both the superficial layers and lower plexus. Group 3 demonstrated two distinguished blood layers separated at a depth of about 300 microns with a pronounced dip of blood content in this area (Figure 3a, red curve). The ratio of blood content at the two characteristic depths z1 = 120 µm and z2 = 570 µm, corresponding to the local maxima of red curve in Figure 3a, demonstrated statistically significant differences between the three age groups and increased with age (Figure 3b). The reconstructed 3D OA images after the Frangi-filter were used to construct the 3D graphs describing the vascular net. Typical en face OA angiograms presented in the projection of the maximum amplitude (Figure 4a-c) and OA sinograms built in the projection of the maximum amplitude onto the surface perpendicular to the scan axis (Figure 4d-f) with superimposed constructed 3D graphs (white lines) are shown in Figure 4 for the three age groups. The reconstructed 3D OA images after the Frangi-filter were used to construct the 3D graphs describing the vascular net. Typical en face OA angiograms presented in the projection of the maximum amplitude (Figure 4a-c) and OA sinograms built in the projection of the maximum amplitude onto the surface perpendicular to the scan axis (Figure 4d-f) with superimposed constructed 3D graphs (white lines) are shown in Figure 4 for the three age groups.

Results and Discussion
The calculation of the parameters of the vascular network structure was carried out for each 3D data set. Various parameters obtained from the 3D OA angiographic images were averaged over the three age groups; the results are summarized in Figure 5. One can see that the image intensity and vessel intensity increase with age (Figure 5a), while the vessel density decreases with age (Figure 5b). At the same time, the number of nodes in the reconstructed vessel skeleton decreases with age, which can be explained by a decrease in the number of small convoluted vessels.
The corresponding p-values (probability values) for the analyzed parameters differences between groups 1, 2, and 3 that were statistically significant (p < 0.05) are shown in Table 1.
Statistically significant age-related differences between groups 1 and 3 were revealed in the total intensities summated over the entire image and over the vascular mask only, as well as in the number of branches in blood vessels and total vessel length. Photonics 2022, 9, x FOR PEER REVIEW 6 of 10 The calculation of the parameters of the vascular network structure was carried out for each 3D data set. Various parameters obtained from the 3D OA angiographic images were averaged over the three age groups; the results are summarized in Figure 5. One can see that the image intensity and vessel intensity increase with age (Figure 5a), while the vessel density decreases with age (Figure 5b). At the same time, the number of nodes in the reconstructed vessel skeleton decreases with age, which can be explained by a decrease in the number of small convoluted vessels.
The corresponding p-values (probability values) for the analyzed parameters differences between groups 1, 2, and 3 that were statistically significant (p < 0.05) are shown in Table 1. Table 1. Analysis of parameters variance between different groups by the Duncan's test. Differences considered statistically significant (р < 0.05) are shown in bold. Group 1: <35 years, group 2: 35-49 years, group 3: ≥50 years. p12 is a p-value between group 1 and 2, p13 is a p-value between group 1 and 3, p23 is a p-value between group 2 and 3.   A decrease with increasing age was revealed for the following parameters: the total number of branches, the total number of nodes, the vessel density, the total length of all vessels, and the average length of one vessel. The average rib length varied insignificantly. The total intensity over the entire OA image in the group of volunteers over 50 years old exceeded that in the group of volunteers < 35 years old by 1.3 times (Figure 5a). The vessel intensity in group 2 was higher than in group 1 by 1.3 times on average, and the values in group 3 were higher than those in group 2 by 1.3 times on average, resulting in a 1.7 times difference between group 1 and group 3. The total vessel length and the number of branching nodes in the group of volunteers over 50 years old were lower than those in the group of volunteers < 35 years old by 1.2 times.

Parameter
As a result, five parameters showed statistically significant differences between groups 1 and 3: the total intensity over the entire image ( Figure 5a) and over the vascular mask (Figure 5a), the total vessel length (Figure 5b), the number of branches (Figure 5d), and the ratio of blood content (Figure 3e). In addition, the total intensity over the vascular mask demonstrated a statistically significant difference between groups 2 and 3, while the ratio of blood content demonstrated a statistically significant difference between groups 1 and 2. Overall, the parameters connected with intensities increased with age, while the indicators characterizing the number of blood vessels decreased with age.
In paper [21], the authors revealed, based on OCT angiography images, that vessel density decreased with age; however, the imaging depth was up to 350 µm limited to the upper plexus. In this study, the vessel density was calculated in all imaging volume up to 2 mm and was demonstrated to decrease with age ( Figure 5c); however, no statistically significant difference was revealed. However, the ratio of blood content between the lower and upper (above 400 µm) plexus demonstrated a statistically significant difference, which was in line with the OCT observations [21]. Color-encoded images (Figure 2) show that the density of the superficial blood vessels decreased with age, which was also consistent with [21]. Deceleration of angiogenesis due to impaired VEGF (vascular endothelial growth factor) signal transmission is considered among the reasons for the decrease in vascular density [1,34]. It is important to note that the gradual decrease in average vessel length and branching with age, revealed in our work, accompanies an increase in the OA signal level. A loss of superficial capillary loops, an increase in the transparency of the skin, and the availability of deep vessels for observation accompanying the aging process [35], inhibition of blood flow, and stasis [36] may be the reasons for the revealed increase in the OA signal level with age. Statistically significant age-related differences between groups 1 and 3 were revealed in the total intensities summated over the entire image and over the vascular mask only, as well as in the number of branches in blood vessels and total vessel length.
A decrease with increasing age was revealed for the following parameters: the total number of branches, the total number of nodes, the vessel density, the total length of all vessels, and the average length of one vessel. The average rib length varied insignificantly. In this paper, the experiments with Fitzpatrick skin types above III were not performed. Since a recent publication [37] demonstrates certain limitations of OA imaging of deep vessels due to high melanin content in Fitzpatrick skin types IV and V, we expect that a high concentration of melanin in the superficial skin layer could serve as a limitation of imaging deeper vessels and, as a result, to performance of the vessel analysis.
Thus, we can successfully use optoacoustic diagnostics to assess blood vessels in healthy subjects, to monitor premature skin aging, and to identify various stages of angiopathies, in particular chronic venous and arterial insufficiency, regardless of age. In addition, the developed analysis can be used for multispectral OA monitoring of blood vessels, which will allow for verification of arterial and venous pathologies.

Conclusions
Age-related changes in the skin of healthy volunteers were analyzed based on OA angiographic images. As a result of the quantitative analysis of three-dimensional OA images based on the skeletonization of the vasculature net, the statistically significant age-associated variations between group 1 (age < 35 y.o.) and group 3 (age > 50 y.o) were revealed for the following parameters: the intensity of the OA images, total vessel length, and number of branches in the vasculature skeleton. In addition, the ratio of blood content at different characteristics depths corresponding to the upper and lower blood plexuses demonstrated a statistically significant difference between group 1 (age < 35 y.o) and groups 2 and 3 (age > 35 y.o.). The proposed algorithm for vasculature analysis can be effectively used for automatic assessment and monitoring of age-related vascular changes in the skin and underlying tissues based on OA diagnostics. In the future, when assessing the state of the vascular architectonics of peripheral tissues, it will be possible to differentiate age characteristics and angio-associated diseases, more reliably assess the severity of the pathological changes, and predict the progression of the disease.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/photonics9070482/s1, Figure S1: OA angiograms built in the projection of the maximum amplitude onto the skin surface of palms of healthy volunteers from groups 1 (first raw), group 2 (second raw) and group 3 (third raw). All bars are 1 mm.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.