Low frequency oscillations assessed by diffuse speckle contrast analysis for foot angiosome concept

An angiosome refers to a 3D tissue volume that is vascularized by a single artery and is a relatively new concept that is useful in vascular surgery; however, the direct relationship between arterial blood flow and micro-perfusion is still controversial. Here, we propose a diffuse speckle contrast analysis (DSCA), which is an emerging tissue perfusion monitoring modality, to investigate the correlations among low frequency oscillations (LFOs) measured from different areas on the feet of healthy subjects. We obtained reproducible results from the correlation analyses of LFOs, and their physiological implications were discussed. In order to confirm the changes in the frequency oscillations, we analyzed and compared the power spectral density changes due to heart rate variability in the electrocardiographic signal during reactive hyperemia and head-up tilt protocols.

www.nature.com/scientificreports/ creases according to cuff-occlusion and release of occlusion, corresponding to autoregulation of blood flow, from which we determined that the DSCA measurement was sensitive to hemodynamic changes in the human foot. Figure 1(c) shows the sectioned timeline BFI data for the baseline (supine), occlusion, and release periods after detrending, normalization, and band pass filter (BPF; pass band frequencies: 0.01-0.15 Hz). Since all BFIs on the four channels fluctuated in a similar manner (i.e., showed LFOs), it was hard to recognize the highly correlated pair intuitively. Figure 1(d) shows the six cross correlations from the four channels for the baseline, occlusion, and release periods, respectively. With the exception of the cross-correlation analysis for the occlusion period, the peaks from all six pairs were distributed around zero lag. In the occlusion period, since the amplitudes of the BFI and LFO decreased owing to the decrease of blood flow supply and autonomic nervous activity, the meaningful relationships disappeared, and we were unable to perform quantitative analyses of the cross correlations. Figure 1(e) illustrates how to obtain the maximum value and its lag time from the cross-correlation graph. As shown in Table 1, we analyzed the maximum values for each of the six cross correlations for the baseline and release periods. Then, the ranks and lag times for the maximum values were estimated; from the high ranks among the six relationship pairs, we noted that three pairs (Ch.2 ⊗ Ch.3, Ch.2 ⊗ Ch.4, and Ch.3 ⊗ Ch.4) were highly correlated regardless of the impact of the reactive hyperemia protocol. Compared to the previous hypothesis that Ch.1 ⊗ Ch.3 and Ch.2 ⊗ Ch.4 relationships will show relatively high cross correlations; the obtained result was not fully compatible with the foot angiosome concept.
From the cross-correlation analysis in this study, it is possible to estimate the similarities between the tissue perfusion data acquired from the optical channels by measuring how large their values are. Moreover, the time latencies (i.e., lag times) for hemodynamic responses can be estimated via the lead-lag time analysis. Briefly, from the locations of the maximum values with respect to zero lag, the leading or lagging of a specific channel pair can be determined. The negative sign in the lag time indicates that the second channel of the pair is faster, and vice versa. The meaning of zero in lag time analysis is that there is no lag latency between the two channels in the pair within our system resolution. Since the data sampling speed of the DSCA is 3 Hz, each data point of the lag number indicates a time duration of 0.33 (= 1/3) s, as shown in Fig. 1(e).
DSCA results of statistical analysis in all subjects. For the cross-correlation data from all the subjects (N = 10), we performed statistical analyses with the averages and standard deviations of the ranks for the maximum values from among the six cross correlations at two different periods and postures, as summarized in Table 2. From the rank analysis, we found that there were no significant differences compared to the previous representative results (Table 1), thereby showing high ranks for the average maximum values for three relationships (Ch.2 ⊗ Ch.3, Ch.2 ⊗ Ch.4, and Ch.3 ⊗ Ch.4). Furthermore, depending on baseline/release periods and supine/tilt postures, the ranks of the cross correlations were not significantly different from each other.
Once the statistical analysis for the lag time was performed with the average and standard deviation calculations for all subjects, the deviations were observed to be too high to obtain meaningful results because the lag time for each subject was different. Instead, as shown in Table 3, we counted the number of occurrences of each sign (+ : positive, − : negative, 0: zero, as shown in the top part of Fig. 1(e)) in the sixty outcomes for the baseline and release periods in all subjects (i.e., 60 = 2 (two different periods) × 10 (ten subjects) × 3 (repeat measurement 3 times). Then, we calculated the percentages for each sign depending on the supine and 90° head-up tilt postures. From this analysis, we derived two reproducible results related to hemodynamic responses among all the channel pairs, based on the angiosome concept. The first result was that Ch.1 and Ch.3 in the MPB angiosome Table 1. Maximum values, its rank, and lag times from six cross-correlation analysis according to baseline and release periods at supine posture in one subject's foot.

Power spectral density (PSD) Changes depending on reactive hyperemia. Representative results
in one subject. In order to analyze the PSD distribution changes for the VLF, LF, and HF ranges with respect to reactive hyperemia, we simultaneously measured the tissue perfusion and heart rate signals using DSCA and ECG, respectively. Figure 2 shows the representative results measured by DSCA (left side in Fig. 2) on the foot and by ECG (right side in Fig. 2) on the whole body from a single subject. In the DSCA measurements, we investigated the frequency oscillation changes of the VLFs (0.001-0.04 Hz), LFs (0.04-0.15 Hz), and HFs (0.15-0.4 Hz) rather than focusing on the change of the LFOs (0.01-0.15 Hz), in order to obtain the spectrum of VLF and to compare the PSD changes of LFs and HFs obtained by ECG. Firstly, we performed timeline selection and detrending/normalization without BPF to acquire entire oscillations during the periods of baseline (90° head-up tilt), occlusion, and release in raw BFI data from the channel 4. Then, we achieved the power spectral analysis with Welch's method as a non-parametric approach in the frequency range from 0 to 0.5 Hz, as shown in Fig. 2 , respectively. In order to examine the relative PSD distribution changes among three frequency bands for reactive hyperemia, based on the proportion of the absolute PSD, we performed pie chart analysis of the PSD distributions, as shown in Fig. 2(d). The relative PSD distribution of VLFs decreased by 26% in the occlusion period and recovered to the baseline values in the release period. In contrast to VLFs, the relative PSD distributions of the LFs and HFs increased by 10% (LFs) and 16% (HFs), respectively. For the ECG results, we first performed timeline selection and detrending/normalization with the raw ECG data according to reactive hyperemia, as shown in Fig. 2(e). Next, we analyzed the HRV signals by estimating the variations in the heart rates and R-R intervals from peak detection in the raw ECG, as displayed in Fig. 2(f). Figure 2(g) shows the results of the power spectral analysis using an autoregressive model as a parametric approach from the HRV in the frequency domain. In contrast to the PSD distribution of the DSCA, two distinct peaks were seen in the LFs and HFs. In the baseline period, the absolute PSDs of the LFs and HFs were 738.28 and 235.64 [ms 2 ], respectively. The absolute PSDs of the LFs and HFs were respectively 596. 13 2 ] in the release period. As shown in Fig. 2(h), we analyzed the relative PSD pie charts with LFs and HFs, except for the VLFs, because the absolute PSD distribution of the Table 3. Percentages for the number of signs of lag time at both baseline and release periods in all subjects for repeat measurement 3 times, according to supine and 90° head-up tilt postures. Bold numbers stand for the highest value among three different situations for the lag time expression. www.nature.com/scientificreports/ VLFs is relatively dominant among three frequency bands, and its origin has not been fully discovered in the ECG study. The changes in PSD distributions for the LFs showed 15% decrease between baseline and occlusion, and 12% increase between occlusion and recovery. For the PSD distribution of the HFs, this increased by 15% for the occlusion period.
Results of statistical analysis in all subjects. As shown in Table 4, we performed statistical analyses with the PSD distributions for the three different frequency ranges, and the data were acquired by both DSCA from channel 4 and ECG with respect to reactive hyperemia. Table 4 shows average and standard deviation of relative PSD distributions from DSCA and ECG measurements. In the DSCA results, because of high standard deviations of the absolute PSD distributions between inter-and intra-subjects, we analyzed relative PSD distributions as the proportion of absolute power within total frequency range. The relative PSD of VLFs among three frequency bands showed the highest changes from baseline to release. The changes in the PSD distributions for the VLFs and LFs were related to magnitude changes in the raw BFIs owing to the effects of autoregulation of blood flow. In the ECG results of Table 4, the absolute PSD distribution of the LFs changed following cuff-occlusion and release, while the absolute power of the HFs did not show any significant differences. Comparing the two www.nature.com/scientificreports/ results by DSCA and ECG for the reactive hyperemia protocol, both slight decreases in LF band were correlated between the relative (DSCA)-and absolute (ECG) PSD distributions, while the distributions in HF band were not correlated. Figure 3 shows the relative PSD distribution changes among three frequency bands by DSCA and two frequency bands by ECG according to reactive hyperemia. These results were similar to the previous representative results (Fig. 2(d) and (h)) obtained from one subject. The relative PSD for VLFs and LFs decreased and increased from baseline to release in both modalities. The relative PSD distribution changes have similar tendencies between the VLFs and LFs obtained by DSCA and the LFs obtained by ECG. From the PSD distribution change analysis from DSCA measurements, it is worth noting that the PSD distribution for VLFs is remarkably reduced and increased following reactive hyperemia. Therefore, we can state that the frequency oscillation measurements by the DSCA system are especially sensitive to PSD distribution changes of the VLFs in specific tissues, which cannot be estimated from the ECG methodology.
PSD changes depending on 90° head-up tilt protocol. Representative results in one subject. Several years ago, Malik et al. 13 reported that the relative PSD distribution of LFs measured by ECG exceeded that of HFs at the 90° head-up tilt, compared to the supine posture. In order to estimate the PSD distribution change of entire frequency oscillations (VLF, LF and HF) achieved by DSCA according to the head-up tilt protocol, we simultaneously performed the PSD analysis with the HRV signal, as shown in  Fig. 4(c), for the pie chart analysis of the relative PSD distributions among three frequency bands, the VLF distribution were changed in the tilt posture (VLF: 38% increase; LF: 18% decrease; HF: 20% decrease) compared to that of the supine state.
In the ECG results ( Fig. 4(d-f)), we observed similar tendencies for the PSD distribution changes of LFs and HFs as those in Malik et al. 13 . For the absolute PSD distributions, the LFs and HFs were respectively 430.83 and 365.63 [ms 2 ] at the supine posture and 422.79 and 156.45 [ms 2 ] at the tilt posture. In the PSD pie analysis from HRV ( Fig. 4(f)), the relative distribution of the LFs increased by 19% according to the posture change.
Results of statistical analysis in all subjects. Table 5 shows the statistical analyses of the relative (DSCA)-and absolute (ECG) PSD distributions for the three different frequency ranges, with respect to the supine and headup tilt postures. In DSCA, we analyzed the proportions of the absolute PSD within total frequency range, which is the same procedure in Table 4. The relative PSDs for both of VLFs and LFs at the tilt posture were higher than that of the supine posture, and the relative proportion of VLFs showed the highest change among the three frequency bands. However, the relative PSDs for the HFs were nearly unchanged. In the ECG results, the highest change in absolute power was noted for the PSD distribution of HFs. Figure 5 shows the pie charts for the relative PSD distributions among three frequency bands (DSCA) and two frequency bands (ECG). From the DSCA Table 4. Average and standard deviation of relative (DSCA)-and absolute (ECG) PSD distributions for VLF, LF, and HF, according to baseline, occlusion, and release periods in all subjects during 90° head up tilt posture. Note that DSCA data was obtained from channel 4. *P < 0.05 compared to baseline period.     www.nature.com/scientificreports/ result of Fig. 5(a), we noted that the relative PSD distribution of HFs decreased by 2%. From the ECG result of Fig. 5(b), the PSD distribution of the HFs decreased significantly. Although the change of the PSD distribution of the HFs in DSCA measurement was little, a decreasing tendency was found from the both results.

Discussions
In this study, we explored cross-correlation analysis with LFOs obtained using the DSCA system based on the foot angiosome concept 5 . From the cross-correlation analysis, it is observed that cross correlations do not strongly support the foot angiosome concept in terms of the micro-perfusion oscillations. As shown in Table 2, by analyzing the cross-correlation amplitudes among the four optical probes, we found that the only relationship between two probes in the LPB angiosome followed a previous hypothesis. For the angiosome concept to be valid, the Ch.1 ⊗ Ch.3 and Ch.2 ⊗ Ch.4 pairs should both have high maximum values among the six relationships because tissue perfusions in the same angiosome must share the same artery. Only the Ch.2 ⊗ Ch.4 pair located in the LPB angiosome space was compatible with our hypothesis. For slight inconsistencies in the average ranks for maximum values of the cross correlations with respect to the baseline/release periods and tilt/supine postures, one possible explanation is that reactive hyperemia and posture change may alter the origins of the LFOs, such as the myogenic or metabolic theory of autoregulation 26 . As shown in Table 3, in the lag time analysis of the cross correlations, we found two reproducible results. The first finding was that the LFOs in both channels 1 and 3 in the MPB angiosome space were faster than those in channels 2 and 4 in the LPB angiosome, which is compatible with arterial paths in the sole of the foot, as expressed by the red lines in Fig. 1(a). Since the PTA branches out from the heel of the foot to the MPB and LPB, and the MPB is closer to the CB of the PTA, the blood supplies to areas monitored by channels 1 and 3 involved in the MPB angiosome arrive earlier than those of the LPB angiosome (areas monitored by channels 2 and 4). The second finding was that in the relationships between both channels in the same angiosome (i.e., Ch.1 ⊗ Ch.3 and Ch.2 ⊗ Ch.4), high percentages for the lag time signs were ambiguous with respect to the supine/tilt postures. Although this is not compatible with physiological arterial blood flow, there may be two possible explanations: (1) In the same angiosome space, the blood flow may be fast enough that it is hard to distinguish the time latency between two channels within the sampling speed resolution. (2) The plausibility of the foot angiosome concept being incorrect may be attributable to the circulatory anatomies 27,28 . As shown in Fig. 1(a), the arteries in the regions of channels 1 and 2 meet each other between the LPB and MPB angiosome spaces at the top of the foot. Since the arteries that form the two different angiosomes are linked, the arterio-arterial anastomosis eventually renders their boundary blurry so that the accuracy of the lag latency analysis in the same angiosome may decrease. Depending on the supine and tilt postures, although the lag time analysis in the same angiosome seems unclear, the evident time difference in the hemodynamic responses between the different angiosome spaces seems to support the foot angiosome concept at the microcirculation level. Note that the validity of the angiosome concept is therefore still controversial 10-12 .
Based on Table 1, considering another correlation analysis, we quantified six Pearson's correlation coefficients which can estimate temporal similarity of the pairs. Table 6 shows the ranks between six cross correlations and six Pearson's correlation coefficients in a single subject, when the lag time is zero. Both ranks by cross-correlation and Pearson's correlation coefficients had similar tendencies at the baseline and release periods, although they were not exactly the same each other. This is because the coefficients of the cross-correlation analysis have been normalized. Compared to those of the baseline and release periods, the six Pearson's correlation coefficients at the occlusion period showed lower values close to zero, resulting in − 0.0695 (Ch.1 ⊗ Ch.  Table 7, from statistical analysis (N = 10) for the ranks of six Pearson's correlation coefficients, we observed that high ranks were still in three-relationships (Ch.2 ⊗ Ch.3, Ch.2 ⊗ Ch.4, and Ch.3 ⊗ Ch.4), which was also found in the high ranks from the six cross correlations. This implies that the analysis of Pearson's correlation coefficient also does not fully support to the foot angiosome concept. To date, regarding the angiosome concept, some studies have focused on assessment of peripheral tissue perfusion to identify the origin location of diabetic foot ulcers with bench-top optical modalities [29][30][31][32][33][34] , to evaluate the quantitative wound healing with high-resolution imaging modalities such as CT 35,36 , SPECT/CT 37,38 and MRI perfusion [39][40][41] , and to evaluate leg revascularization surgeries Table 6. Comparison of ranks between six cross-correlations and six Pearson's correlation (CORR) coefficients at zero lag and in supine posture according to baseline and release periods from one subject's foot.  11 . To the best of our knowledge, no studies have reported the assessment of cross correlation relationships between the inter-and intra-tissue spaces based on the angiosome concept. The main objective of this study was to explore the foot angiosome concept from the view of the microcirculation level by a novel approach compared to the previous studies. Therefore, we expect that our results, which partially support the angiosome concept, might help to consolidate the foot angiosome in tissue perfusion, with an understanding of the physiological background. From the PSD change analyses using DSCA, we confirmed that the DSCA system could effectively distinguish tissue perfusion and its fluctuations on the sole of the foot following reactive hyperemia and the 90° head-up tilt protocol, as effects of the autoregulation of blood flow and spontaneous systemic changes. Although we included the results for PSD distribution analyses with only channel 4 data, the results for all the other channels showed similar trends, and the sensitivity in channel 4 was the best among the four channels. As shown in Table 4 and Fig. 3, in the reactive hyperemia experiment, the remarkable reduction in PSD distribution of the VLFs among the three frequency bands at the cuff-occlusion period implies that the DSCA system probes the vasomotion tone inside the tissues in the VLF band and not the noise information. As shown in Table 4, we observed that the PSD distribution of HFs in the DSCA measurements following the reactive hyperemia were relatively increased because of dominant reduction of the relative PSD in VLF band. As shown in Table 5 and Fig. 5, however, the PSD distributions of VLFs and LFs measured by DSCA in the 90° head-up tilt protocol experiment, increased at the tilt posture and that of VLF band was mostly changed. In Table 5, the relative PSDs of HFs by DSCA did not change significantly from supine to head up tilt posture. As reported in Tong et al. 17 , we also experienced that all subjects had no clear peaks at the HF band which are closely related with the respiratory frequency. Moreover, it was limited to scope the HF band signal with the relatively low sampling rate of DSCA at the regional tissue.
In Table 4 and Fig. 3, for the analysis of PSD changes by ECG, in the reactive hyperemia protocol, we found that the LF distribution in the occlusion period somewhat decreased compared to the HF distribution for the baseline, thus showing that the ratio of LFs to HFs before and after occlusion were 2.57 and 1.56, respectively. By occluding the femoral artery in the leg, the reactive hyperemia protocol mimics a situation of peripheral arterial disease (PAD), which is a circulatory problem associated with narrowed arteries and reduced blood supply to the limbs. After cuff occlusion, the decrease in PSD distribution of the LFs was similar to a previous report that the LF distribution was lower in patients with cardiovascular morbidities with PAD than patients without PAD, as observed by Chen et al. 42 . As shown in Fig. 4(d-f), in the 90° head-up tilt protocol, the change in LF distribution by ECG was similar to that observed by Malik et al. 13 , where the absolute power of the LFs slightly decreased at tilt posture whilst the relative distribution of LFs considerably increased. From simultaneous measurements using DSCA and ECG, it was observed that there was a relationship between the LF distribution obtained from HRV analysis of ECG and VLF and LF ( ≈ LFO) distributions obtained using DSCA.
The reason why the results of the PSD analyses estimated using DSCA and ECG were not similar was attributable to the measurement scope of each system being different. The ECG measurements are used to estimate all systemic oscillations in the human body. In contrast, the DSCA estimates the relative blood flow and perfusion oscillations in specific tissues. Although the origin of the VLFs has not yet been fully realized from ECG studies, the basis by which we obtained the VLFs (0.001-0.04 Hz) and LFs (0.04-0.15 Hz) using DSCA might be affected by endothelium-related metabolic (0.008-0.02 Hz), neurogenic (0.02-0.05 Hz), and myogenic (0.05-0.15 Hz) regulations in skin blood flow oscillations 43 . In this work, the DSCA probed the LFOs at the sole of the foot. Therefore, we suggest that the ECG is useful for diagnosis of various diseases, such as cardiovascular diseases, diabetes, and so on, but the DSCA system has broad impact for diagnosing localized diseases (e.g., CLI, PAD, cancers, diabetic foot ulcer) especially at the tissue level, for which specific diseases cannot be diagnosed by ECG.
Thus far, to our knowledge, there are only a few reported studies on tissue perfusion fluctuation analysis with optical systems, which can measure hemodynamics, that focus on analyzing the relative changes or differences in blood oxygenation, volume, and flow between healthy subjects and patients. If the previous measurement parameters and LFOs (or VLFs; LFs) are together applied to patients, we expect to be able to scope these diseases more precisely and gain more resources in clinical situations. To this end, we are planning to implement a highspeed diffuse optical system that can estimate the tissue metabolic rate of oxygen consumption for scoping a broader power spectrum within the tissue.

Conclusions
We performed a pilot study for validation of the foot angiosome concept with a correlation analysis of the LFOs obtained using the DSCA system; our results in this study do not strongly support the angiosome concept. However, because the angiosome concept is gradually gaining more attention in terms of treatments for patients with CLI, our results can help improve insights on tissue perfusion assessments. As shown in this work, we verified Table 7. Average and standard deviation of the rank for Pearson's correlation coefficients at baseline and release periods in all subjects, according to supine and 90° head-up tilt postures.

Posture
Periods Ch.1 − Ch. www.nature.com/scientificreports/ that the DSCA can estimate VLF and LF ( ≈ LFO) changes in specific tissues during arterial cuff-occlusion and 90° head-up tilt protocols. In future work, we intend to compare the results between diabetic patients and healthy subjects to assess the key factors distinguishing these two groups, so that early diagnosis of diabetic complications may be possible. To date, some research groups have studied the early diagnosis of diabetic complications using optical systems with near-infrared sources [44][45][46] . However, these studies have focused on the time-averaged hemodynamic variables and disregarded LFO information. Since tissue perfusion and its oscillations are important biomarkers in human vascular health, it is expected that these factors will offer new insights into diagnosis and therapeutic monitoring of diabetic complications.

Materials and methods
Participants. The present experimental protocol was approved by the Institutional Review Board (IRB) of Daegu Gyeongbuk Institute of Science and Technology (approval no. DGIST-190829-h-072-01). Ten healthy subjects (mean age: 26.5 years old; 7 males and 3 females) participated in this study on tissue perfusion measurements on the foot. Prior to and during the experiments, we performed in accordance with relevant guidelines/ regulations for the approved research and obtained the informed consent for the experimental protocol from all participants, and they were asked to avoid smoking and drinking caffeinated and alcoholic beverages for minimizing the vasoactive effects and to limit movements as much as possible for minimizing motion artifacts in the DSCA and ECG signals.
Experimental situation and protocol. Figure 6(a) shows the experimental schematic for the hemodynamic response measurements from the sole of the foot using DSCA and via HRV measurements from the human body using ECG. The subjects' thighs were wrapped using a blood-pressure cuff to temporarily block the blood flow from the femoral artery to the foot. For simultaneous measurement of tissue perfusion and HRV, the four-channel optical probes of the DSCA system were attached to the right foot, and the three-channel ECG probes were attached to left and right wrists and right ankle, as shown in Fig. 6b and c. The reason behind the simultaneous measurements using DSCA and ECG was to evaluate the relevance between the frequency oscillations within specific tissues and systemic oscillations in the entire body. www.nature.com/scientificreports/ For the ECG record, three ECG sensor patches were linked to an ECG sensor module (AD8232, Heart Rate Monitor, EDUINO). The ECG signal obtained by the ECG module was sent to a computer using an Arduino Uno module. In the DSCA configuration, a laser diode (DL-785-100S, CrystaLaser, 785 nm, 100 mW) was linked to a 1 × 4 coupler with multi-mode optical fibers to extend the system to four channels. To reduce the harmful effects from large light power illumination, all lasers at the ends of the illumination optical fibers were adjusted to apply around 6.5 mW 47 . The scattered light was detected using a CCD camera (F-033B, Stringray) through four single-mode optical fibers. The optical probes (pairs of laser source and detector fibers) were attached on the foot through flexible polydimethylsiloxane (PDMS) for optical guidance. All of the inter-optode distances for the probes were set to 10 mm, so that the penetration depth was approximately 5 mm. More details on the probe configuration can be found in our previous paper 23 .
As shown in Fig. 6(b), all the measurement positions in each subject's foot were determined by the individual foot sizes. Channels 1 and 3 (Ch.1 and Ch.3) were located on areas vascularized by the MPB and channels 2 and 4 (Ch.2 and Ch.4) were located on areas vascularized by the LPB to investigate possible correlations between LFOs from different positions on the foot angiosome. We hypothesized that the optical probe pairs (Ch.1 ⊗ Ch.3 and Ch.2 ⊗ Ch.4) located at the same angiosome will have higher correlations than those from different ones (Ch.1 ⊗ Ch.4, Ch.1 ⊗ Ch.2, Ch.2 ⊗ Ch.3, and Ch.3 ⊗ Ch.4) since the microcirculation in any specific angiosome will share the same origin of LFOs.
In the experimental protocol, three periods were defined as baseline, cuff-occlusion, and release of occlusion, as shown in Fig. 6(d). In the baseline period, to analyze the PSD changes of VLFs, LFs, and HFs according to the head-up tilt protocol 13 , the measurements involve a supine posture for 5 min and 90° head-up posture for 5 min. Between the two postures, there was a 2 min measurement interval that was removed from analysis as a waiting period for the settling time to achieve stable vasomotion because dramatic posture changes can affect vascular tone. In the occlusion period, a cuff-occlusion pressure of about 200 mmHg was induced for 2 min to change the physiological blood flow to the leg and to verify that the DSCA instrument could detect hemodynamic changes according to reactive hyperemia. The release period is the post-occlusion of the pressure for 5 min. Therefore, the measurement time for each subject was about 19 min. For better statistical analysis, these measurement procedures were repeated six times with gaps of 2 h rest for each subject. In each measurement phase, the sequence of supine and head-up tilt postures was changed for a more accurate analysis of the PSD changes.
Data processing. In the ECG record, the signal was acquired via serial communication using an Arduino Uno with 9600 board rates, such that the data acquisition speed was 100 Hz. In the DSCA system, the image acquisition from the CCD camera was achieved via the LabVIEW program with an exposure time of 10 ms and sampling rate of 60 fps. The BFI within a volume of tissue was calculated using the speckle intensity on the region of interest (ROI). For detailed theoretical descriptions of the DSCA, we refer the readers to published literature 20 . In signal processing for speckle contrast calculation, there are two types of analysis for the temporal and spatial domains. For the temporal domain analysis, the speckle contrast can be obtained from the correlation of a group of pixels over a fixed position and at successively different times. In contrast, for the spatial domain analysis, the speckle contrast can be obtained from statistics over n × n pixels out of the entire image at a fixed time. In this work, we utilized the temporal domain analysis to simultaneously acquire multi-channel data with a 9 × 9 pixels window as the ROI. The temporal speckle contrast ( K t ) is defined as in Eq. (1).
where σ t and I are the temporal standard deviation and average of light intensities over 20 frames, respectively. Therefore, the data acquisition speed was 3 Hz when using 20 images at 60 fps. It is possible to increase the sampling rate above 3 Hz using a high-speed camera, but the camera used in this study was limited to 60 fps because of the maximum frame rate of the CCD camera. In addition, we did increase the sampling rate above 3 Hz because the LFOs did not exceed 0.15 Hz. Finally, we defined 1/K 2 t as the BFI with arbitrary unit, which varies linearly with physiological blood flow.

Data analysis.
Cross-correlation analysis from DSCA signals. After tissue perfusion measurements on the foot with the DSCA system, we analyzed the cross-correlation relationships among the four channels, as shown in Fig. 7. Within the time-series tissue perfusion raw data, we chose the timeline selections according to the baseline, occlusion, and release periods in the supine and head-up tilt postures. The reason for selections of the occlusion and release periods was to ascertain that the cross-correlation relationship still remained among the four probes after baseline period. Before applying the cross-correlation analysis, the selected matrices (4 channels × n-time-series: 4 × n matrix) were processed by detrending and normalization to remove linear trends and obtain better cross-correlation analysis. Next, the detrended and normalized data are transferred into a band pass filter with the pass band frequencies from 0.01 Hz to 0.15 Hz to deal with the LFOs range. Finally, we employed an unbiased cross correlation procedure using MATLAB software, for the cross-correlational analysis.
Heart rate variability (HRV) analysis from electrocardiographic signal. The HRV, which is obtained from the ECG record, represents one of the promising biomarkers to express oscillations in consecutive cardiac cycles 13 . Among several available measures of the HRV, we analyzed the variations of the instantaneous heart rate and R-R intervals. We obtained 256 consecutive R-R values in the baseline and release periods for the HRV analysis in the frequency domain. In the occlusion period, 128 consecutive R-R intervals were utilized because of limited measurement time (2 min).
(1) K t = σ t �I� Scientific Reports | (2020) 10:17153 | https://doi.org/10.1038/s41598-020-73604-0 www.nature.com/scientificreports/ Power spectral density (PSD) analysis from processed blood flow index (BFI) and HRV data. In both the processed BFI (detrended and normalized) and HRV data, the timeline according to the baseline, occlusion, and release periods in supine and head-up tilt postures were selected and processed by detrending and normalization. In order to acquire the entire oscillations in human body, we did not apply the filters such as low/high/ band pass filters in both the BFI and HRV data. The time-series signals for the processed BFI and HRV data were converted into frequency domain signals. Next, spectral analysis was performed using (1) Welch's method (triangular window; 900 samples) in a non-parametric approach 13 for the processed BFI and (2) the autoregressive model (order 20) in a parametric approach 13 for the HRV data in MATLAB software. Since the hemodynamic system in the sole tissue has not been modeled yet, the PSD analysis from DSCA was achieved by the nonparametric method. Then, both signals by DSCA and ECG were analyzed in the three frequency bands, such as VLFs (0.001-0.04 Hz), LFs (0.04-0.15 Hz), and HFs (0.15-0.4 Hz). In the case of HRV data, a pie chart analysis was done with the absolute power of the two frequency bands, to confirm the relative distributions between two frequency bands, LFs and HFs. In the processed BFI, we examined the relative PSD distributions in the three frequency bands from 0 to 1.5 Hz (sampling rate of DSCA: 3 Hz), due to the high variances of absolute PSD distributions from inter-and intra-subjects.
Statistical analysis. For statistical analysis, the averages and standard deviations of the ranks of the maximum values of the cross correlations and the absolute PSD distributions of the VLFs, LFs, and HFs were obtained with ten samples for the number of subjects (n = 10) and six-repeated measurement samples. A two-tailed paired t-test was performed with the significance level as a p-value < 0.05.