1.1.

Conventional Methods for the Assessment of Skin Microvascular Function

The examination of the microcirculation is a promising tool for the establishment of biomarkers of cardiovascular disease (CVD) risk, which may help to improve the diagnosis, prognosis, and early preclinical assessment of vascular pathologies.¹ However, the direct assessment of the microcirculation in vivo is limited by the requirement of invasive methods based on catheterization, e.g., coronary angiography and intracoronary Doppler.¹ In contrast, the examination of the microcirculation from peripheral organs, i.e., skin, represents a better modality to investigate microvascular function based on the use of noninvasive imaging laser technologies.^1–3

The conventional methods for the evaluation of skin microvascular function are laser Doppler and laser speckle contrast imaging (LSCI) in combination with reactive tests. Laser Doppler techniques rely on the detection of the Doppler shift of red or infrared light backscattered by moving erythrocytes for the estimation of the average blood flow in small tissue volumes.^2,4,5 Within these methods, laser Doppler flowmetry (LDF) is used for the continuous monitoring of microvascular perfusion through single-point laser probes in direct contact with the skin.⁵ An advantage of LDF is the high temporal data acquisition that allows studying the nonlinear fluctuations of blood flow.⁶ This can be achieved by applying the continuous wavelet transform (CWT) spectral analysis on LDF time series,^7,8 which reveals indirectly the activity of specific biological factors in the microcirculation, e.g., endothelial cells (ECs), vascular smooth muscle cells (VSMCs), and the local sympathetic innervation. However, relevant limitations of laser Doppler and LSCI are the impossibilities to visualize blood vessels and the expression of data in arbitrary units making it possible only relative assumptions based on the combination of perfusion measurements with reactive tests.^4–6 An example of a widely used reactive test for the examination of human cutaneous microcirculation is the postocclusive reactive hyperemia (PORH), which evaluates the increase in blood perfusion from the forearm or fingertips after a temporary period of ischemia stimulated by artificial occlusion of the brachial artery.^5,6

1.2.

Techniques for the Examination of Cutaneous Microvascular Structure

Skin microvascular structure can be visualized in vivo by several methods based on optical microscopy, three-dimensional (3-D) photoacoustic imaging or optical coherence tomography (OCT). A common relevant limitation of the optical microscopy techniques is the generation of two-dimensional (2-D) images with uncertain depth resolution.^9,10 3-D photoacoustic imaging is based on the generation of pressure acoustic waves in all the directions employing a pulsed laser light absorbed by chromophores in the skin, which act as contrast agents to visualize microvessels.^11,12 The main benefit of photoacoustic imaging is the full light penetration depth.¹³ However, the method suffers from poor spatial resolution, even though this aspect has been improved by implementing the optical resolution photoacoustic microscopy (OR-PAM).⁹ Moreover, the method requires a coupling medium in direct contact with the tissue, which is a drawback for several applications, i.e., image-guided brain surgery,^14,15 wounds/burns assessment,¹⁶ and ophthalmology.^14,15 This might also be a limitation to assess the activity of skin microvessels, which is sensitive to the local cutaneous temperature¹⁷ and may be affected by the temperature of the coupling medium. However, the recent implementation of contact-free PAM^18–20 is promising to overcome this shortcoming. Additionally, some limits of photoacoustic imaging have been solved with the further implementation of the OCT method to allow microvessels imaging. OCT is a noninvasive contactless technology employing a broadband near infrared light source to provide a low coherence gate for the selection of backscattered light from different tissue depths, producing in-vivo 2-D and 3-D structural images.²¹ The main advantage of OCT is the generation of cross-sectional images comparable to conventional optical biopsies with micrometer-level structural resolution.^22,23 However, the technique is limited by low penetration depth ($\sim 1$ to 2 mm).²⁴

1.3.

Concurrent Evaluation of Microvascular Structure and Function

Although the visualization of microvascular structure is powerful for the detection of specific organizational patterns or damage to the microvasculature that may be associated with CVD, laser Doppler and LSCI are still the preferential choices for the examination of skin microvascular function. The reason is that functional abnormalities in the microcirculation precede structural changes and clinical manifestations, thus representing a hallmark of the early stages of the pathology suitable for the preclinical prediction of CVD risk.^1,2 However, laser Doppler and LSCI are limited by both ambiguous spatial and depth resolution. Therefore, the ideal tool for a complete assessment of the microcirculation would be a technology providing simultaneous functional and structural information. Extensions of OCT utilize the Doppler shift (DOCT)²⁵ and the speckle variance (svOCT)²⁶ principles for providing a contrast mechanism to visualize blood flow, or processing techniques relying on the separation of static areas and flow dynamic regions of the tissue, i.e., Doppler optical microangiography (DOMAG).²⁷ The combination of OMAG with the PORH reactive test has been described previously as an optimal method for providing quantitative functional data of the human peripheral skin microcirculation.²⁸

Another technique for the in vivo imaging of the microcirculation is correlation mapping OCT (cmOCT). This method employs correlation coefficient statistics to determine blood flow from the intensity of the reflected OCT signal,^9,10 taking advantage of the time-varying speckle of moving scatters in skin regions containing active blood vessels and the constant reflectance of the stationary scatters in the bulk tissue lacking of blood vessels.^9,10 Therefore, the correlation map obtained from the cross-correlation analysis is able to visualize the location of blood flow.^9,11 The advantages of cmOCT compared with DOCT, svOCT, and DOMAG, respectively, are the apparent angle independence of the technique,¹⁰ nonrequirement of previous knowledge of the structure to determine flow regions, and very fast data processing to obtain the flow maps.⁹ However, cmOCT is limited by high decorrelation caused by noisy background regions in the images, which can be misinterpreted as a blood flow signal. This necessitates the use of a binary mask to suppress these regions⁹ leading to reduced quality images compared with those generated with OMAG or svOCT techniques.

1.4.

Aims of the Study

The objectives of this research were providing evidence on the potential application of cmOCT for assessing microvascular function and demonstrating for the first time that optical coherence tomography based angiography (OCTA) techniques are suitable to monitor vasomotion and the activity of specific biological factors of the microcirculation. This was achieved by generating sequential cmOCT blood flow maps from the human forearm at different cutaneous depths during a 25-min PORH reactive test and applying the CWT spectral analysis on continuous blood flow signals reconstructed from cmOCT maps.

The results displayed that cmOCT is able to detect functional flow dynamics and vasomotion oscillatory patterns reflecting the activity of the microvascular endothelium and local sympathetic innervation. Our results demonstrate that the study of microvascular function by OCTA techniques is not restricted only to the general observation of flow dynamics but can be extended for gaining information on the behavior of specific biological components in the microcirculation network.

2.1.

Study Population

The experiments were performed on eight women and four men with an average age of $33\pm 10$-years-old, who signed an informed consent form according to the rules of the Declaration of Helsinki. The measurements were collected in the morning time in a laboratory with a constant temperature ($20°\mathrm{C}\pm 1°\mathrm{C}$), allowing at least 15 to 20 min of acclimatization of the subjects. The volunteers were asked to observe a fasting period of 4 to 6 h prior to the examination to avoid physiological differences related to nutritional habits. The measurements were performed from the left volar forearm maintained at the heart level while the volunteer was laying in supine position in a clinical bed. A custom-built soft support was used to fix the arm, ensuring a static and comfortable position during the test to avoid bias in the experimental data due to movement artifacts (Fig. 1). The skin site for the measurements was chosen while avoiding hairy and injured regions. A pressure cuff was placed in the upper part of the left arm for performing the PORH test during the acquisition period of the images. The occlusion of the brachial artery to stimulate PORH response was induced by inflating the pressure cuff at 200 mmHg.

Fig. 1

Experimental setup. A custom-built mount was used to immobilize the arm as much as possible during testing. A contact support was applied in the skin surface with gentle pressure to stabilize the scan location between datasets avoiding mismatches. The OCT probe was placed above the contact support to collect data from the chosen region. The image shows also additional equipment (an LDF optical probe and a camera) that is not relevant for the research presented in this work.

2.2.

OCT Device for Images Acquisition

A commercial OCT system (TEL2200C1, Thorlabs Inc.) was used for imaging 3-D volumes of the cutaneous tissue. This spectral domain (SD) device consists of a superluminescent diode (SLD) with a center wavelength of 1300 nm, which supports an imaging depth up to 3.5 mm (in nonscattering samples), an axial scan rate from 5.5 to 76 kHz, and an axial resolution of $5.5\mu \mathrm{m}$ in air. The sample arm contains a LSM003 (Thorlabs Inc.) scanning lens providing a transverse resolution of $13\mu \mathrm{m}$ in air. Figure 1 shows the experimental setup employed to collect OCT data. For each volunteer tested in the study, an initial 3-D volume of $1.50\times 1.50\times 1.82{\mathrm{mm}}^3$ was acquired at a speed of 48 kHz with a scanning density of $1024(\text{length})\times 1024(\text{width})\times 512(\text{depth})\text{pixels}$, to visualize the microvascular structure of the volar forearm. The time employed for the acquisition of the structural volume was $\sim 51\mathrm{s}$, which is too long for the assessment of blood flow functional dynamics. For this reason, the examination of microvascular dynamics was performed by reducing the sampling density to $256\times 256\times 512\text{pixels}$ that can be acquired at a temporal rate of $\sim 6\mathrm{s}$, allowing the fast scanning of sequential 3-D volumes suitable for the assessment of microvascular function. However, an important requirement for visualizing blood flow by cmOCT is a dense sampling below the lateral resolution of the adjacent B-frames of the OCT volume, which ensures a strong correlation between the stationary scatters of consecutives B-scans. Therefore, to maintain a spatial separation between the B-frames suitable for the cmOCT analysis the size of the sampled area was reduced to $0.75\times 0.75\times 1.82{\mathrm{mm}}^3$, which means a spatial separation of $2.9\mu \mathrm{m}$. Although the possibility of performing each scan in $\sim 6\mathrm{s}$, the volumes were acquired with a sampling frequency of one volume every 16.6 s because the software provided with the OCT device (ThorImage OCT 4.3, Thorlabs Inc.) required $\sim 10.6\mathrm{s}$ for saving each volume. In total, 90 sequential low-resolution volumes were acquired from the forearm of each volunteer in a time window of $\sim 25\min$ while performing the PORH reactive test. An initial 10-min acquisition of 36 volumes was performed to determine the baseline flow, then 18 volumes were recorded for a 5-min ischemia period stimulated by blocking blood flow through the brachial artery, and finally 36 volumes were acquired for further 10 min to monitor flow dynamic changes during the hyperemic response.

2.3.

cmOCT Data Processing

The cmOCT analysis to generate blood flow maps at different depths from each 3-D OCT volume was performed by employing the same in-house designed Java algorithm and principles described by Enfield et al.⁹ and Jonathan et al.¹⁰ The software is able to perform the correlation analysis between consecutive adjacent 2-D B-scans extracted from the 3-D OCT volume, by estimating the cross correlation of a square grid from the first B-scan 1 ($I_1$) to the same grid from the adjacent B-scan 2 ($I_2$), according to the following equation:

A $7\times 7$ kernel size of the grid was employed to generate the correlation maps, which represents an optimal compromise of producing good sensitivity to blood flow without loss of spatial resolution due to kernel size. As mentioned in the introduction, the cmOCT maps are affected by a background noise with low correlation coefficients between consecutive B-scans that make the identification of the flow regions difficult. The noise was removed by processing the correlation maps with a structural mask produced from the original OCT structural images by the application of a kernel blur and binary threshold. The maps were presented as maximum intensity projection (MIP) images, which allow at the same time a better visualization of microvessels morphology and location, and the detection of the maximum flow intensity. To examine blood vessels from different skin microvascular beds (capillaries and arterioles/venules), the cmOCT maps were extracted from different depths (80 to $180\mu \mathrm{m}$, 180 to $280\mu \mathrm{m}$, and 300 to $400\mu \mathrm{m}$) of the OCT volumes, i.e., in the enface plane. This was obtained by determining the position of skin surface from the structural images and using the superficial location as a reference to establish the various depths of the underlying tissue.

2.4.

Reconstruction of Blood Flow Time Series

To investigate vasomotion dynamics by wavelet analysis, a time series of microvascular blood flow was reconstructed at different tissue depths from the cmOCT maps using the same methodology previously described by our group.²⁹ First, the average blood flow was estimated from each of the 90 sequential MIP maps generated at different tissue depths during the reactive hyperemia test. This was achieved by calculating the mean value across the region of interest (ROI) of each map. Then, the mean values of flow extracted from the sequential MIP maps were employed as discrete data points to reconstruct a continuous curve, covering the 25-min time window used in the experiment. The reconstruction was performed by the piecewise cubic spline interpolation²⁹ method using the cftool instrument in MATLAB R2015a (The Mathworks Inc.). Specifically, we have created a vector containing the 90 sequential average values of blood flow extracted from the cmOCT maps and a time vector defining an interval of 16.6 s between each pair of discrete data points corresponding to the sampling frequency adopted during the temporal acquisition of the 3-D OCT volumes. More details on the principles of spline interpolation can be found in our previous publication.²⁹

2.5.

CWT Spectral Analysis

The CWT is a widespread method for the investigation of blood flow fluctuations, revealing the contribution of specific biological components in the microcirculation to the nonlinear dynamic changes of microvascular flow. In this work, the application of CWT on blood flow signals reconstructed from the cmOCT maps allowed the identification of oscillations in the frequency intervals 21 to $52\times {10}^{-3}\mathrm{Hz}$, 9.5 to $21\times {10}^{-3}\mathrm{Hz}$, and 5 to $9.5\times {10}^{-3}\mathrm{Hz}$. According to the literature,^8,30 these oscillations reflect the activity of local sympathetic nerves, the ECs activity nitric oxide (NO)-dependent, and the ECs activity NO-independent.

The CWT analysis was performed according to the principles described by Stefanovska et al.^7,8,30 The method is based on the use of the Morlet mother wavelet function $\psi$ scaled by a factor $s$ and a time $t$ to provide an adequate frequency window, which is shifted along the time domain of the analyzed signal for obtaining an optimal time–frequency resolution.^7,8,30 The CWT spectrum of a signal $g(u)$ is calculated as follows:^7,8,30

In this work, a central frequency $f$ of the Morlet wavelet equal to 1 was used,³⁰ noting that the scale parameters $s$ and $f$ were related according to the following equation,³⁰ $s=1/f$. Specifically, a logarithmic scale array of 346 scales ranging from 0.5 to 200 with Morlet factor 1.03 was used, which corresponds to a frequency interval between 5 and $2000\times {10}^{-3}\mathrm{Hz}$. The minimum size of the time window used in the study was $\sim 200\mathrm{s}$ to allow the investigation into the slowest oscillations of interest (frequencies up to $5\times {10}^{-3}\mathrm{Hz}$). A timescale of minutes is sufficient to detect oscillatory phenomena with a period of 200 s. However, while for periodic signals a single period is sufficient to observe the phenomenon under investigation, this is not true for nonlinear quasiperiodic time series, i.e., blood flow signals.⁸ In this case the period of a specific phenomenon fluctuates constantly, thus the trace should be long enough to contain at least several periods.⁸ For this reason, here a 25-min time window was used for data collection to sample a sufficient number of periods for observing oscillations up to $5\times {10}^{-3}\mathrm{Hz}$.

3.1.

Visualization of Microvascular Structure

To select the cutaneous areas for the study of functional dynamics, an initial 3-D OCT volume was acquired from a $1.50\times 1.50\times 1.82{\mathrm{mm}}^3$ region in the volar forearm of each volunteer to visualize microvascular structure. The acquisition was performed in $\sim 51\mathrm{s}$ employing $1024\times 1024$ A-scans, and the MIP flow maps were obtained at different tissue depths by cmOCT data processing using a $7\times 7$ kernel grid. Figure 2 shows examples of structural enface MIPs at three different depths obtained from a volunteer, describing the typical skin microvascular structure reported by Enfield et al.,⁹ even if in this case the size of the analyzed region was smaller.

Fig. 2

Structural visualization of cutaneous microvessels from the volar forearm. (a) 2-D structural B-frame from a $1.50\times 1.50\times 1.82{\mathrm{mm}}^3$ 3-D OCT volume. The region marked in green indicates the epidermal surface (0 to $40\mu \mathrm{m}$), the yellow area represents roughly the location of the dermal–epidermal junction (80 to $180\mu \mathrm{m}$), and the red region outlines a part of the dermal plexus (300 to $400\mu \mathrm{m}$). (b) Enface MIP projection of the epidermal surface at 0 to $40\mu \mathrm{m}$. (c) Enface MIP projection of the dermal–epidermal junction around 80 to $180\mu \mathrm{m}$, where the capillary loops are located. The small yellow square indicates the $750\times 750\mu \mathrm{m}$ area chosen for the study of capillaries functional dynamics. (d) Enface MIP projection of the dermal plexus at a depth of 300 to $400\mu \mathrm{m}$ characterized by a horizontal network of microvessels. The small red square outlines the $750\text{-}\times 750\text{-}\mu \mathrm{m}$ region chosen for the study of arterioles/venules functional dynamics.

Figure 2(a) shows a structural 2-D B-scan of the OCT volume, showing a longitudinal view of the tissue regions from which the MIPs were generated. The depths were selected using the skin surface, marked with a green line, as a reference for the estimation of the specific microvascular beds. Figure 2(b) shows the MIP map from the epidermal surface at 0 to $40\mu \mathrm{m}$ that, as expected, does not show any flow region due to the absence of blood vessels. Figure 2(c) shows the MIP generated from the dermal–epidermal junction, which is located around a depth of 80 to $180\mu \mathrm{m}$. The image shows small blood vessels with vertical orientation corresponding to the finger-like capillary loops arising from the ascending microvessels of the dermal layer and appearing as colored spots in the blue background. The function of the capillary loops is the nourishment of the epidermal layer that is lacking of blood vessels. Figure 2(d) shows an MIP projection of the skin layer around 300 to $400\mu \mathrm{m}$ corresponding to the dermal plexus. This microvascular bed is characterized by a network of arterioles and venules with horizontal orientation from which, respectively, the nutrients are delivered to the capillaries and the waste products directed to the venous system. To study the microvascular dynamics from the dermal–epidermal junction and dermal plexus regions, an enface plane of $0.75\times 0.75{\mathrm{mm}}^2$ was selected from the structural images to perform sequential temporal acquisitions of 3-D OCT volumes during PORH test. Figures 2(c) and 2(d) show examples of $750\text{-}\times 750\text{-}\mu \mathrm{m}$ regions selected for the study of capillary loops and arterioles/venules functional dynamics, respectively, marked by yellow and red squares.

3.2.

Assessment of Flow Dynamics from cmOCT Maps

Ninety sequential volumes ($256\times 256\times 512\text{pixels}$, $0.75\times 0.75\times 1.82{\mathrm{mm}}^3$) of the small yellow/red area selected in Figs. 2(c) and 2(d) were acquired from each volunteer during a 25-min PORH functional test at a sampling rate of 16.6 s. The first 36 volumes were acquired during a 10-min resting period to determine the baseline flow, 18 volumes were collected during a 5-min occlusion of blood flow through the brachial artery, and 36 volumes were acquired for 10 min after the removal of the occlusion for monitoring blood flow changes during PORH response. The volumes were processed by cmOCT to obtain sequential flow maps of the dermal–epidermal junction (80 to $180\mu \mathrm{m}$) and dermal plexus (300 to $400\mu \mathrm{m}$) regions for tracking the temporal flow changes of capillaries [Fig. 3(a)] and arterioles/venules [Fig. 3(b)], respectively.

Fig. 3

Sequential cmOCT flow maps generated at different tissue depths from the OCT volumes acquired during PORH test. (a) Temporal enface MIPs from the $750\text{-}\times 750\text{-}\mu \mathrm{m}$ yellow area selected in Fig. 2(c) for the study of skin capillaries microvascular dynamics at a depth of 80 to $180\mu \mathrm{m}$. (b) Temporal enface MIPs from the $750\text{-}\times 750\text{-}\mu \mathrm{m}$ red area selected in Fig. 2(d) for the study of cutaneous arterioles/venules microvascular dynamics at a depth of 300 to $400\mu \mathrm{m}$. The maps were obtained from $0.75\times 0.75\times 1.82{\mathrm{mm}}^3$ volumes sampled every 16.6 s with a scanning density of $256\times 256$ A-scans. The baseline 1 to 4 MIPs describe the flow at rest detected prior to the occlusion stage of the PORH functional test, the ischemia 1 to 2 and 3 to 4 MIPs refer, respectively, to the flow detected immediately after the beginning of the occlusion process and before the removal of the occlusion, and the PORH 1 to 4 MIPs outline the flow detected during PORH response. The images contain residual background noise, which is clearly visible in the ischemia MIPs at both 80- to $180\text{-}\mu \mathrm{m}$ and 300- to $400\text{-}\mu \mathrm{m}$ cutaneous depths.

For simplicity, Figs. 3(a) and 3(b) show only the cmOCT maps extracted from 12 volumes of the total 90 temporal samples, which represent the key steps of the reactive test. The maps were characterized by residual background noise represented by small snow-like dots [see the ischemia MIPs in Figs. 3(a) and 3(b)], which may be confused and misinterpreted as a capillary blood flow at the 80- to $180\text{-}\mu \mathrm{m}$ depth. Nonetheless, the relevant increase of flow observed in specific points of the MIP maps at 80 to $180\mu \mathrm{m}$ during PORH response [Fig. 3(a), PORH 1 to 4] was helpful to identify the specific regions of the images containing the capillaries. Additionally, because all the MIPs were characterized by the same degree of background noise at all the stages of the reactive test, this made it possible monitoring reliably the relative dynamic changes of blood flow over time. Indeed, although the images were characterized by poor resolution and residual background noise compared to the structural images, they clearly outlined the temporal flow changes associated with the different stages of the functional task. In detail, the capillaries did not show relevant changes during occlusion [Fig. 3(a), ischemia 1 to 4] compared with baseline [Fig. 3(a), baseline 1 to 4], whereas the network of blood vessels in the dermal plexus displayed a clear reduction and disappearance of blood flow during occlusion [Fig. 3(b), ischemia 1 to 4] compared with baseline [Fig. 3(b), baseline 1 to 4]. These observations probably reflect the delivery of the residual flow from arterioles to capillaries during the 5-min ischemia period, ensuring the nourishment of the epidermal tissue under hypoxic stress conditions. In contrast, a massive increase of flow intensity typical of the hyperemic response was detected during PORH for both capillaries and arterioles/venules [Figs. 3(a) and 3(b), PORH 1 to 4], associated with the vessels reactivity for the restoration of a normal microcirculation after removal of the occlusion. These results demonstrate that cmOCT is suitable for tracking the general microvascular functional dynamics.

3.3.

Reconstruction of Time Series from cmOCT Maps

In this study, the cmOCT flow maps were also employed to detect the key contribution of biological components in the microcirculation to the rhythmic variations in vessels diameter (vasomotion) in response to PORH stimulus. The goal was achieved by reconstructing continuous time series of blood flow at different depths from the intensities of the sequential temporal cmOCT maps (see Sec. 2.4 in Methods for more details), and processing the reconstructed signals with the CWT spectral analysis that is a widespread method for the study of microvascular oscillations. Figure 4(a) shows an example of continuous time series obtained by using the average flow extracted from 90 sequential cmOCT MIP maps as discrete data points for the reconstruction of a curve covering the 25-min period of the PORH test. The method employed for the reconstruction process was the piecewise cubic spline interpolation that is advantageous because no oscillation artifacts are introduced between each pair of discrete data points. This means that the obtained time series reflects the real dynamics of the process under investigation. The reconstructed curve in Fig. 4(a) clearly shows the typical trend of the reactive hyperemia test, characterized by a reduction of blood flow during ischemia and a massive increase of flow during PORH response reaching a peak before the gradual restoration of the normal microcirculation. This is a further confirmation that cmOCT is able to detect microvascular functional dynamics.

Fig. 4

Reconstruction of continuous time series from cmOCT maps. (a) Example of curve reconstruction by piecewise cubic spline interpolation of the discrete average flow values estimated from the ROI of the 90 sequential cmOCT MIP maps generated at a tissue depth of 300- to $400\text{-}\mu \mathrm{m}$. A time interval of 16.6 s between each pair of data points was applied to cover the 25-min period of the PORH test. The flow is expressed in arbitrary units (AU). (b) Average absolute blood flow from 12 healthy subjects during reactive hyperemia test reconstructed at three different tissue depths. The flow is expressed in AU. (c) Average relative blood flow from 12 healthy volunteers during PORH test reconstructed at three different depths. The flow is expressed as a function of baseline.

For each volunteer enrolled in the study, the microvascular flow was reconstructed from cmOCT maps generated at three different depths: 80 to $180\mu \mathrm{m}$ (dermal–epidermal junction), 180 to $280\mu \mathrm{m}$ (upper dermal plexus), and 300 to $400\mu \mathrm{m}$ (median dermal plexus). Figures 4(b) and 4(c) show, respectively, the absolute and relative average reconstructed blood flows at the selected depths for the 12 tested individuals. As shown in Fig. 4(b), the absolute flow reflects the increasing density of blood vessels in the deepest layers of the cutaneous tissue. Indeed, the absolute flow intensity was greater in the deepest microvascular beds at all the stages of PORH test. The relative average blood flow in Fig. 4(c), expressed as a function of the baseline, provides a more reliable comparison of flow dynamics from different tissue depths in response to the reactive test. The graph shows the same findings discussed for Figs. 3(a) and 3(b), consistent with the maintenance of a normal flow in the capillaries and a gradual reduction in the flow intensity going deeper in the tissue during ischemia. This may reflect the delivery of the residual flow from the deeper arterioles network to the upper capillaries for feeding the epidermis during ischemia. In contrast, PORH response was characterized by higher flow intensity and peak response in the deepest vascular beds, reflecting the path of flow directed from the deepest network of arterioles to the superficial capillaries and a major vasodilator response in the deepest layers associated with greater size and density of microvessels.

3.4.

CWT Spectral Analysis of the Flow Signal Reconstructed from cmOCT Maps

The study of vasomotion at different tissue depths was performed by applying the CWT spectral analysis on the reconstructed time series. The method allowed the identification of oscillations in the frequency ranges 21 to $52\times {10}^{-3}\mathrm{Hz}$ (neurogenic), 9.5 to $21\times {10}^{-3}\mathrm{Hz}$ (endothelial NO-dependent), and 5 to $9.5\times {10}^{-3}\mathrm{Hz}$ (endothelial NO-independent).^8,30 Figure 5 shows an example of CWT analysis of the blood flow signal reconstructed from a tissue depth of 180 to $280\mu \mathrm{m}$. Figure 5(a) shows the CWT scalogram graph that describes the distribution of the wavelet power of the signal in the time–frequency domain, visualized as a gradient colored map. The graph clearly shows an increase of power (increase of brightness) during PORH response in all the oscillatory frequency ranges under investigation, reflecting the activation of the endothelial and sympathetic control mechanisms associated with vasodilation. The time-averaged spectrum in Fig. 5(b) discriminates the maximum wavelet power peaks at various frequency intervals corresponding to the specific blood flow oscillators. This allows monitoring the activity of each oscillator and making comparisons between different subjects, by extracting values of the maximum amplitude $A_i$ and energy $E_i$ (area under the curve) from each wavelet peak during the various steps of the reactive test. For instance, in Fig. 5(b) an increase in the wavelet maximum amplitude and power of all the oscillators was clearly observed during PORH (black line) compared with ischemia (dark red line) and baseline (dark blue line).

Fig. 5

CWT analysis of a continuous blood flow signal reconstructed from cmOCT maps generated at a depth of 180- to $280\text{-}\mu \mathrm{m}$. (a) CWT scalogram showing the distribution of the wavelet power in the time–frequency domain using a colored map ranging from dark (low energy) to bright (high energy) colors. The colored areas in the scalogram outline the “cone of influence,” which is a time-frequency area where distortions of the CWT due to the limited duration of the time series are not relevant.³¹ In contrast, the transparent areas at the bottom-right and bottom-left of the graph represent regions external to the cone in proximity of the time edges of the signal, where the wavelet transform is characterized by boundary effects making the estimations from these regions ambiguous.³¹ (b) Time-averaged wavelet spectra computed for all the stages of the reactive hyperemia test. The plots distinguish the CWT maximal amplitude peak and energy (area under the curve of the peak) associated with each specific blood flow oscillator. The neurogenic (21 to $52\times {10}^{-3}\mathrm{Hz}$), ECs NO-dependent (9.5 to $21\times {10}^{-3}\mathrm{Hz}$), and ECs NO-independent (5 to $9.5\times {10}^{-3}\mathrm{Hz}$) oscillators were identified.

Figure 6 shows the median time-averaged CWT spectra obtained from the 12 healthy volunteers, by wavelet analysis of the cutaneous flow signals reconstructed at 80 to $180\mu \mathrm{m}$ (yellow line), 180 to $280\mu \mathrm{m}$ (blue line), and 300 to $400\mu \mathrm{m}$ (red line) tissue depths, respectively. The CWT spectral power was normalized as a function of the minimum baseline power, to compare the relative spectral energy between different cutaneous depths during PORH test. All the wavelet oscillators showed greater power with increasing depth at all the stages of the reactive hyperemia task [Figs. 6(a)–6(c)]: baseline, ischemia, and PORH. These observations may be explained by a major number of microvessels involved in vasomotion in the deepest cutaneous layers due to greater vessels density.

Fig. 6

Median relative time-averaged CWT spectra of reconstructed flow signals from 12 healthy subjects during PORH test. (a) Baseline median CWT spectra of blood flow time series reconstructed from cmOCT maps generated at 80 to 180, 180 to 280, and 300 to $400\mu \mathrm{m}$. (b) Median CWT spectra during ischemia. (c) Median CWT spectra during PORH response. All the spectra were normalized by the respective minimum baseline spectral power to compare the relative increase of wavelet energy during the hyperemic response between different cutaneous depths.

The wavelet energy $E_i$ and amplitude $A_i$ of all the oscillators, extracted from the time-averaged CWT spectra according to the principles described by Shiogai et al.,⁸ were significanty increased during PORH response compared with baseline for all the analyzed tissue dephts (Table 1), which is in agreement with the results reported for the spectral analysis of LDF signals recorded during the hyperemic challenge.³²

Table 1

Comparison of the amplitude and power of the wavelet oscillators between baseline and PORH time points. The variables were nonparametric, thus data are expressed as median (interquartile range). The amplitude Ai represents the maximum wavelet power in the interval of interest, whereas the energy Ei is the overall wavelet power in the frequency range under investigation. The comparison was performed at all the analyzed depths: 80 to 180, 180 to 280, and 300 to 400  μm. The p-values (p) were estimated by paired Wilcox test performed using the software R-Studio.

Although OCTA methods are flow-based imaging tools not suitable to monitor the architecture of blood vessels, the increase of blood perfusion observed in this study during PORH response should reflect vasodilation dynamics. Indeed, literature reports have demonstrated that the increase of flow observed during cutaneous PORH is mediated by the activation of vasodilator mechanisms regulated through the major activity of the endothelium and the sympathetic nerves.⁵ This is also supported by the results of the wavelet spectral analysis of cmOCT blood flow signals discussed in Table 1. The wavelet analysis is a method suitable for the indirect detection of vasomotion dynamics by its application on blood flow signals that do not contain architectural information, i.e., LDF signal.^7,8,30 Therefore, the increase in the wavelet spectral power observed for the endothelial and neurogenic oscillators should reflect indirectly the rythmic vasorelaxation process mediated by ECs, VSMCs, and sympathetic nerves for the fast reperfusion of the tissue after an ischemic period. Damages to the microvascular endothelium may impact negatively the hyperemic response with a reduction of the endothelial-mediated vasodilation, which can be reflected in the wavelet spectrum as a reduced increase of the endothelial wavelet power.

3.4.1.

Mechanisms at the basis of the vasomotor perturbations of skin blood flow signal

The endothelial and neurogenic oscillators and the myogenic oscillation, which was not detected in this study, are classified as low-frequency dynamic patterns of flow motion associated with local vasomotor phenomena occurring on the timescale of blood distribution to the tissue surrounding microvessels.^8,33 These slow vasomotor rhythms cause changes in the vascular diameter as part of the local regulation of microvascular flow, and they are reflected in the blood flow signal as naturally occurring perturbations representing a correlation of volumetric flow to a third-order blood pressure wave.³³ The vasoconstriction or vasodilation induced by these regulatory rhythms has an effect on red blood cells (RBCs) velocity with a consequent decrease or increase of velocity leading to lower or higher blood perfusion to tissue,³⁴ respectively. These oscillatory phenomena have been detected also on capillary blood flow, as confirmed by our results (Fig. 6), probably reflecting the degree of vasomotion in the terminal arterioles connected to the capillary loops.³⁴

The mechanism at the basis of the vasomotor perturbations of blood flow signal is the mechanical contraction or relaxation of VSMCs, due to depolarization or hyperpolarization induced by cellular rhythmic fluctuations of ${\mathrm{Ca}}^{2+}$ concentrations throughout voltage-dependent channels in the plasma membrane and sarcoplasmic reticulum.^35,36 This process may be modulated through different regulatory mechanisms (Fig. 7) that are characterized by a different oscillatory rhythm reflected by the low-frequency vasomotor waves detected in the blood flow signal.

• 1. Myogenic mechanism (52 to $145\times {10}^{-3}\mathrm{Hz}$) based on spontaneous contraction/relaxation of VSMCs in response to changes of the intra-arteriolar blood pressure.^8,34,35

• 2. Neurogenic mechanism (21 to $52\times {10}^{-3}\mathrm{Hz}$) based on a feedback between VSMCs and the local sympathetic or parasympathetic innervation.³⁴ The nerves are activated by changes in the intravascular pressure and/or tissue perfusion,³⁴ which cause the release of substances targeting VSMCs to regulate the vascular diameter.⁸

• 3. Endothelial mechanisms (5 to $21\times {10}^{-3}\mathrm{Hz}$) based on a feedback between VSMCs and ECs.^8,34 The activity of ECs may be stimulated by the concentration of metabolic mediators in the blood or by the shear stress,^8,37 which cause the release of vasoactive substances modulating the relaxation/contraction of VSMCs.

Fig. 7

Mechanisms causing the vasomotor perturbations of blood flow signal. The slow vasomotor waves of skin microvascular flow are caused by changes in the vascular diameter mediated through three main mechanisms occurring with different rhythms: (a) The spontaneous myogenic contraction/relaxation of VSMCs in the microvessel wall induced by changes in the intra-arteriolar blood pressure.^8,34,35 (b) The neurogenic feedback between the local nerves and VSMCs stimulated by changes in intra-arteriolar blood pressure or blood perfusion,³⁴ which cause the release of substances from the nerves modulating the contraction/relaxation of VSMCs.⁸ (c) The myoendothelial feedback between the ECs and the VSMCs.^8,34 The activity of ECs may be stimulated by changes in the concentration of metabolic substances in the blood or by shear stress,^8,37 which cause the release of vasoactive agents to modulate VSMCs tone, i.e., NO, PGs, EDHF (mediated by release of metabolites) and endothelin-1 or electric signals, e.g., EDHF (mediated by gap junctions signaling). The endothelial oscillations of blood flow signal may reflect the rate of release of vasoactive agents targeting VSMCs, or the rhythmic response of a cellular oscillator in the VSMCs (e.g., ${\mathrm{Ca}}^{2+}$) to these substances.³⁷ NO, nitric oxide; PGs, prostaglandins; EDHF, endothelial-derived hyperpolarizing factor.