2.1.

Subject Profiles

The study protocol was approved by the institutional review board of Gangnam Severance Hospital, Seoul, Korea. Written informed consent was obtained from each subject. A total of 19 subjects (9 diabetic patients, 10 normal controls) underwent near-infrared ICG fluorescence imaging. The age of the subjects ranged from 50 to 65 years (mean: $58.78\pm 5.52\text{years}$) for diabetic patients and 51 to 60 years (mean: $55.90\pm 3.14\text{years}$) for controls. For diabetic patients, exclusion criteria included the presence of acute coronary syndrome (unstable angina, acute myocardial infarction), history of heart failure, and pregnancy/lactation. Exclusion criteria for controls included Raynaud’s syndrome and pregnancy/lactation. Demographic data of the subjects are summarized in Table 1.

Table 1

Demographics of subjects enrolled in the study.

2.2.

Dynamic Fluorescence Imaging

DyFI with ICG was used to measure various vascular parameters, as previously reported.²¹ The ICG fluorescence imaging system for clinical applications was manufactured by Vieworks Corporation (Anyang, Gyeonggi-do, Korea)²² and consisted of a charge-coupled device digital camera (RXD-500, Vieworks Co.) with an 830-nm bandpass filter (RG830; Edmund Optics, Barrington, New Jersey) and 760-nm light-emitting diode (LED) arrays (SMC760; Marubeni America Co., Sunnyvale, California). The sum of the current of the three LED arrays was kept under 1.3 A. Time-series images ($768\times 512\text{pixels}$) of both dorsal feet were obtained while subjects were in a supine position. Each image was taken at 5-s intervals for 600 s, immediately following an intravenous bolus injection of ICG ($0.16\mathrm{mg}/\mathrm{kg}$; Dongindang Pharm. Co., Gyeonggi-do, Korea). Characteristic features of the ICG fluorescence images are presented in Fig. 1, such as the temporal sequence [Fig. 1(a)]. Figure 2(b) shows that each pixel has different ICG pharmacokinetic dynamics. A three-dimensional (3-D) plot of ICG fluorescence dynamics is presented in Fig. 1(c). Data with movement artifacts were excluded from the analyses.

Fig. 1

Characteristics of ICG fluorescence dynamics: (a) temporal sequence of ICG fluorescence images after the intravenous injection of ICG at 0 s; five representative time-points are shown ($768\times 512\text{pixels}$, 120 frames, total 600 s). (b) Time-series of ICG fluorescence intensity at different regions of interest (ROIs) on the foot. Three ROIs were selected: the left hallux (red, ROI1); the vein (green, ROI2); and the peripheral vessel (blue, ROI3). Each pixel has a different ICG pharmacokinetic pattern according to type of vessel. (c) Spatiotemporal profile of ICG dynamics showing the high dimensionality of the raw data. Normalized fluorescence intensity is shown in a 3-D space ($x$-axis, pixels; $y$-axis, frames; $z$-axis, fluorescence intensity). Each pixel has different ICG fluorescence dynamics.

Fig. 2

Flowchart of PCA for ICG fluorescence dynamics. Fluorescence intensities of both feet were extracted using 120 sequential ICG fluorescence images from each subject. For the mathematical analysis, normalized fluorescence intensity was converted to a matrix. The rows of the input matrix correspond to spatially successive pixels, while columns correspond to temporally sequential frames. The covariance matrix was calculated from the input matrix. The components were then extracted as its eigenvector. The eigenvector with the highest eigenvalue was taken as the first PC. The feature vector consisted of the first three PCs. By calculating the inner product of the input matrix with the PCs, high-dimensional raw data could be converted into low-dimensional data without significant loss of information. $i_{\mathrm{pk}}$, fluorescence intensity of raw signal at $k$’th frame of $p$’th pixel; $S$, vectorized fluorescence intensity; $C$, covariance matrix; F.I., fluorescence intensity; $\lambda$, eigenvalue; $U$, eigenvector; $e$, function of feature extraction; $x$, projection data.

2.3.

Principal Component Analysis

A schematic diagram of the ICG dynamics analysis using PCA is shown in Fig. 2. The first step involves loading ICG sequential images, which consist of ICG fluorescence intensities. The region of interest was both feet except the background and shadow. $i_{\mathrm{pk}}$ is the normalized fluorescence intensity at the $k$’th frame of the $p$’th pixel. $I_p$ is the vectorized fluorescence dynamics of the $p$’th pixel in 120 sequential frames [Eq. (1)]. To apply the PCA algorithm, the image pixels were interpreted as a matrix of integers. We made an input matrix $S$ that incorporates the whole spatiotemporal profile, composed of all fluorescence dynamics of all $N$ pixels. Rows of $S$ correspond to the dynamics of every pixel, while columns of $S$ correspond to frames [Eq. (2)].

The covariance matrix $C$ is calculated using Eq. (3), while the eigenvector ($U$) and eigenvalue ($\lambda$) are calculated from $C$ [Eq. (4)].

PCs were extracted according to their eigenvalue. The eigenvector with the highest eigenvalue was considered the first PC (PC1). The PC2 vector was the second highest eigenvalue and orthogonal to PC1. By calculating the inner product between the input matrix $S$ and eigenvector ($U$) sorted by decreasing order, we could project the entire fluorescence dynamics onto PC space ($x$) [Eq. (5)].

We wrote a C++ (Visual Studio 2010, Microsoft) program that can extract ICG fluorescence dynamics from each pixel. PCA was performed using the princomp function included in MATLAB^® software (MATLAB^® 2014b, Mathworks).

2.4.

Statistical Analysis

Statistical differences were analyzed by two-tailed Student’s $t$ test, or Spearman’s correlation where indicated. Data are expressed as the $\mathrm{mean}\pm \text{standard}$ deviation, and a $p$ value of $<0.05$ is considered statistically significant.

3.1.

Distinctive Patterns of Principal Components in Diabetic Patients and Normal Controls

PCA was applied to extract the important elements from raw ICG fluorescence dynamics. From the input matrix, which contains the sequential ICG fluorescence intensities of all pixels on the foot, PCs were selected based on their eigenvalue. The eigenvector with the highest eigenvalue was defined as the first PC. The second PC was defined as the second-highest eigenvalue with a direction uncorrelated to the first PC. Each eigenvector has 120 dimensions because the raw data consisted of 120-frame sequential images. Time to the PC1 maximum coefficient value was faster in normal controls than in diabetic patients, corresponding to general ICG pharmacokinetic patterns [Fig. 3(a) and 3(b)]. The PC2 curves extracted from normal controls had a sharp peak in the early time phase and decreased exponentially after the peak point [Fig. 3(c)]. This feature of PC2 curves in normal controls is similar to that observed with arterial input function (AIF), as an early bolus arrival, a steep rise, and a narrow peak are considered arterial-like features.²³ In diabetic patients, PC2 curves were smoother and more dispersed than in normal controls [Fig. 3(d)].

Fig. 3

Representative PCs time-courses. Each curve was extracted from a different individual subject. The PC2 vector is the second-highest eigenvalue and is orthogonal to PC1. (a) and (b) Average time-course of the first three PCs. PC1 is similar to ICG pharmacokinetic curve. PC2 of controls represents arterial input function. (c) PC2 time-courses extracted from four representative normal controls. A sharp peak is seen in the early time phase (before 100 s); the peak decreases exponentially with image acquisition time. (d) PC2 time-courses extracted from four representative diabetic patients. A delayed peak is observed, with lower coefficients than normal controls.

The meaning of PCs can be revealed by investigating the relationship between pixel values projected onto PCs and other known dynamics features. The dynamics features compared with PCs included blood flow index (BFI), mean transit time (MTT), and time to max ($T_{\max }$). It has been reported that these parameters are related to vascular conditions. BFI and MTT showed significant differences between normal and middle cerebral artery occlusion in mice.²⁴ Modified $T_{\max }$, which means the time from onset of ICG fluorescence to time for maximum intensity (MI), was reported as a diagnosis parameter of RP.¹⁴ We calculated Spearman’s correlation coefficients between pairs of features. The maximum fluorescence intensity and area under curve (AUC) showed positive correlation with the PC1 coefficient ($r=0.74$ and $r=0.73$, respectively, $p<0.0001$, Spearman’s correlation) [Fig. 4(a)]. MTT, calculated as the center of gravity of the dynamics,²⁵ had a strong positive correlation with PC2 ($r=0.90$, $p<0.0001$, Spearman’s correlation). $T_{\max }$ and percentile $T_{\max }$ also have strong correlation with PC2 ($r=0.87$ and $r=0.86$, respectively, $p<0.0001$, Spearman’s correlation). BFI, based on the slope of the rising peak, had a negative correlation with PC2 ($r=-0.80$, $p<0.0001$, Spearman’s correlation) [Fig. 4(b)]. We plotted scatter graphs to visualize the correlation between dynamics features and PCs [Figs. 4(c) and 4(d)]. Our results suggest that PC1 is closely related to the original pharmacokinetic curve (MI and AUC parameters), while PC2 seems to contain functional information related to the BFI and MTT parameters.

Fig. 4

Correlation analysis between various vascular parameters and PCs in normal controls. Spearman’s correlation coefficients with (a) PC1 and (b) PC2. Strong positive correlations were observed between MI and PC1 ($r=0.74$, $p<0.0001$), and MTT and PC2 ($r=0.90$, $p<0.0001$). AUC, area under curve; BFI, blood flow index; $T_{\max }$, time-to-max. (c) and (d) Scatter plots of PC score ($x$-axis) and the vascular parameters ($y$-axis) that showed the highest correlation with PCs in (a) and (b).

3.2.

Explained Variance of Principal Components

Proportion of explained variance for PCs was identified as a quantitative value that could differentiate between normal controls and diabetic patients. Each PC contributes to a proportion of the total ICG fluorescence dynamic, based on its eigenvalue. The percentage of explained variance for each PC can be calculated using Eq. (6), which is the corresponding eigenvalue divided by the total variance ($k$ is the component number; $\lambda$ is the eigenvalue).

Table 2 shows the average variance explained by corresponding PCs. For normal controls, the first, second, and third PCs account for 77.71, 17.82, and 2.65% of the average variance, respectively. For diabetic patients, these same PCs account for 92.36, 5.19, and 1.32% of the average variance, respectively. The total variance explained by the first three PCs in controls (98.18%) and diabetic patients (98.87%) was sufficient to represent the high-dimensional ICG fluorescence data. The noticeable difference between diabetic patients and controls was observed with the PC1 and PC2 pattern (Fig. 5), where the proportion of explained variance was significantly lower in diabetic patients than in controls ($p<0.0001$, two-tailed Student’s $t$ test). Figure 5(c) shows the gap between control and diabetes was widened by calculating the ratio of major two variances, PC1 (%) divided by PC2 (%). This finding suggests that explained variance of the first two PCs may be considered a potential imaging-based biomarker that can differentiate between a normal foot and a diabetic foot.

Table 2

Explained variance of PCs (%).

Fig. 5

A dot plot of variance (%). Data are shown as means (red line) with standard error of the mean (light red area) and standard deviation (light blue area). (a) Variance explained of PC1 (%), (b) variance explained of PC2 (%), and (c) the ratio of two major variances, PC1 (%) divided by PC2 (%).

3.3.

Multichromatically Visualized Principal Components Map

To represent the spatial distribution of PCs, PCA results were mapped against projection values of the dynamics of each pixel. A schematic representation of the mapping method is provided in Fig. 6(a). By calculating the inner product of normalized ICG fluorescence dynamics and PCs, the high-dimensional raw data were converted to low-dimensional data that conserve the important elements. The first three PCs became the new basis after the projection. Each pixel has three projection values on PC1, PC2, and PC3 space. Pseudocolored images were created by merging the data with RGB channels (red, green, and blue channels represented PC1, PC2, and PC3, respectively). In PC1–PC2 spaces visualizing all the dynamics on the hallux and dorsal feet, only PC2 was capable of differentiating the hallux and dorsal feet regions [Fig. 6(b)]. The PC2 score on the dorsal feet was higher than the hallux in normal controls, while the reverse pattern was observed in diabetic patients. The distribution of PC2 scores between normal controls and diabetic patients from Fig. 6(b) is shown in histogram form in Fig. 6(c). Figure 6(d) compares the RGB maps of a normal control and a diabetic patient. The projection values of PC2 were mapped to veins on the dorsal feet in normal controls and the hallux in diabetic patients (Fig. 7). These representative PC distribution color-maps show characteristic patterns for normal controls and diabetic patients, and could potentially be used to detect and evaluate vasculopathy in diabetic feet.

Fig. 6

Projection of data onto PCs. (a) A map with PC-matching RGB channels. RGB mapping creates a visual representation of the spatial distribution of the PCs. Fluorescence intensity was projected onto the first three PCs. Each projected value has PC1, PC2, and PC3 scores. The spatial pixels corresponding to the first, second, and third orthogonal temporal components of the image were visualized as a red-green-blue map. (b) Comparison of scatter plots. By projecting raw ICG fluorescence dynamics onto PC1 and PC2, the data could be visualized as a scatter plot. (c) Histograms of PC2 and (d) a case study of PC distribution. Red, green, and blue channels represent PC1, PC2, and PC3, respectively.

Fig. 7

The pseudo color-map of PC2 score. (a) and (b) A comparison of PC2 distribution with green color-map. The projection values of PC2 were mapped to veins on the dorsal foot in a normal subject [green dashed box in (a)] and the hallux in a diabetic patient [green dashed box in (b)]. (c) and (d) A comparison of PC2 distribution with jet color-map. The projection values of PC2 are higher on veins than hallux in a normal subject. Inversely, the score is higher on hallux than veins in diabetic patients.