2.1.

Tissue Model

To generate realistic training data for the ANN, MC simulations of the light-tissue interaction in a three-layer skin tissue model were performed. This process has been described in detail elsewhere^6,16 and will only be summarized here in the context of MELSCI. Optical properties used in the model are relevant for the laser wavelength 780 nm. The first layer in the tissue model is a bloodless epidermis layer with variable thickness ($t_{\mathrm{epi}}$) and melanin concentration. The absorption coefficient of the epidermis layer (${\mu }_{\mathrm{a},\mathrm{epi}}$) is dependent on melanin and is generally low at the used laser wavelength. The second and third layers are dermis layers with variable blood tissue fractions, blood speed distributions, scattering coefficients, and vessel diameters. The upper dermis layer has a fixed thickness of 0.2 mm, and the lower dermis layer is infinite. The reduced scattering coefficient (${\mu }_{\mathrm{s}}^{\prime }$) is the same in all layers. The absorption coefficient of the dermis layers (${\mu }_{\mathrm{a},\text{dermis}}$) is given by their blood tissue fraction ($c_{\text{blood}}$), oxygen saturations ($S_{\mathrm{oxy}}$), and vessel diameters ($d_{\text{vessels}}$), where the latter affects the absorption through the vessel packaging effect.^19,20 The parameters $c_{\text{blood}}$, $S_{\mathrm{oxy}}$, and $d_{\text{vessels}}$ are modeled to be different in the two dermis layers, where the actual values are calculated from a randomized average value and a randomized difference parameter, to ensure the parameters do not differ too much between the layers.

The parameters of interest in the 100,000 randomized training models are summarized in Fig. 1. The average values for $c_{\text{blood}}$, $S_{\mathrm{oxy}}$, and $d_{\text{vessels}}$ are given in the histograms. The last histogram shows the average speed $⟨v_{\mathrm{RBC}}⟩$ of the RBC speed distributions, described in more detail below. The choice of the distributions was based on three principles; (1) the range of the distributions should cover essentially all expected values in skin tissue for all skin sites and all skin types; (2) allow for less probability of parameter values that are rarely expected, such as ${\mu }_{\mathrm{s}}^{\prime }<0.5{\mathrm{mm}}^{-1}$²¹; (3) have higher probability for ranges where a small change has a large effect on the contrast data, thus exponential decay in the histograms of many of the parameters as Beer–Lambert’s law predicts lower absolute change for high ${\mu }_{\mathrm{a}}$ or large path-lengths (related for example to high $t_{\mathrm{epi}}$).

Fig. 1

Normalized histograms of skin parameter distributions in the tissue models. The $x$-axes are scaled to remove the 1% upper tail.

When modeling speckle contrast, it is common to assume either a Gaussian or Lorentzian distribution for the electric field autocorrelation function,^9,22 but the validity of both assumptions has been questioned.¹¹ Here, we instead model the speed distribution of RBCs, $c_{\mathrm{RBC}}(v)$, directly as a weighted sum of 10 uniform distributions from $0\mathrm{mm}/\mathrm{s}$ to exponentially increasing top speeds, representing vessels with parabolic flow profiles with 10 different average speeds. The weights are randomized for each tissue model to create variations in the shape of the distribution. Examples of speed distributions are presented in Fig. 2. Additional details can be found in Ref. 6.

Fig. 2

Examples of RBC speed distributions $c_{\mathrm{RBC}}(v)$, for three different mean speeds as presented in the legend. The distributions were normalized to have a total probability of 1.

2.2.

Doppler Histogram Calculation

Tissue model simulations were performed for calculating Doppler histograms (the optical Doppler spectrum) as previously described in detail.¹³ The simulations and calculations are summarized in this section.

MC simulations were performed with varying epidermis thickness and tissue scattering covering the ranges found in Fig. 1 for these parameters. In each MC simulation, the total path length of each photon in each of the three layers was stored and used to calculate path-length distributions in each layer. For a given set of the model parameters $t_{\mathrm{epi}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, the path-length distribution in each layer was interpolated from the simulated path-length distributions. Using the path-length distribution in each layer, Beer–Lambert’s law was applied to add the absorption effect of ${\mu }_{\mathrm{a},\mathrm{epi}}$ and ${\mu }_{\mathrm{a},\text{dermis}}$, modifying the path-length distributions. These modified path-length distributions were the first input to the calculation of the Doppler histogram of the model.

A single-shifted Doppler histogram was calculated for light Doppler shifted once by a red blood cell with a certain speed. This was done analytically as described in the appendix of Ref. 13, assuming a random angle between the light and the direction of the red blood cell and a scattering angle based on the Gegenbauer kernel phase function.²³ The phase function parameters were set to $g_{\mathrm{Gk}}=0.948$ and ${\alpha }_{\mathrm{Gk}}=1.0$, resulting in an anisotropy factor of 0.991.²⁴ The Doppler histogram for light shifted $n$ times was calculated as the cross-correlation of the single-shifted Doppler histogram with itself $n-1$ times. Probability distributions of the number of Doppler shifts were given by Poisson distributions based on the pathlengths and fraction of red blood cells. The fact that blood is confined in distinct blood vessels rather than being evenly distributed in the layer, i.e. the vessel packaging effect, was also accounted for when calculating the distribution of the number of Doppler shifts.

Finally, the speed distribution, the path-length distributions, the distributions of the amount of Doppler shifts, and the Doppler histograms for various numbers of Doppler shifts and speeds, were used to calculate a final Doppler histogram $H(f)$. Although the entire process is complex, the calculation scheme is efficient, and it takes less than a millisecond to calculate a Doppler histogram for a specific set of model parameters.

2.3.

Speckle Contrast Calculation

The Doppler histogram $H(f)$ obtained as described in the previous section was converted into multi-exposure contrast $K(T)$ using a mathematical framework that is thoroughly described in Refs. 6 and 25. In short, the electric field autocorrelation function $g^{(1)}(\tau )$ was obtained from the Doppler histogram using the Wiener–Khinchine theorem

The contrast computed by this method is noise-free, making it inappropriate for training models that should be applied to real data containing noise. Therefore, a noise component $K_{\text{noise}}(T)$ was added to the simulated contrast data before training the neural network. The noise component model has been described and validated in Ref. 2. In summary, speckle simulations showed that the noise for each exposure time $T$ could be modeled as a sequence of random contrast offsets depending on the average contrast over all exposure times ${⟨K(T)⟩}_T$ and the noise value at the previous exposure time, i.e.,

2.4.

Biological Zero

When using the laser Doppler technique in measurements at low-flow or zero-flow conditions in skin tissue, for example during brachial occlusion, we have observed a residual high frequency energy in the Doppler power spectrum (ongoing work), corresponding to the well-known biological zero (BZ) effect.^27,28 During occlusion, BZ partially originates from Brownian motion of RBCs and from the redistribution of blood between vessels, possibly including reverse flow in the capillaries.²⁹ In addition, movements of macromolecules in the interstitial space contribute to the BZ.³⁰ Due to the small size of these macromolecules, the scattering angles are large compared to scattering from RBCs, resulting in large Doppler shifts. This addition to the power spectrum cannot be explained by Doppler shifts from RBCs, regardless of which speed distributions of the RBCs are used. Furthermore, the residual energy is not present when measuring on flow of phantoms perfused with whole blood,¹⁴ making us confident in assuming that it does not originate from the blood itself, but from some kind of other tissue movements, well in line with the movement of macromolecules. We observe the same effect when using the MELSCI technique. Similar to the residual we observe in the Doppler spectrum, the shape of the contrast curve in measured data during occlusion cannot be explained by any realistic speed distribution of the RBCs. However, since we observe the same phenomenon using both LDF and MELSCI, we are confident that it does not originate from faulty calibrations or instrument noise.

We must account for this macromolecule effect on BZ when generating training data from the skin models, otherwise, the training data will not be valid for low-flow conditions. This can be done by adding a BZ component $H_{\mathrm{BZ}}$ to the Doppler histogram before computing the contrast. This was empirically modeled as

2.5.

Perfusion

Speed-resolved perfusion can be computed from the modeled speed distribution as

When calculating the modeled perfusion $P_{\text{true}}$, the contribution from each skin layer is considered so that $P_{\text{true}}$ reflects the sampling volume. For example, if the epidermis thickness is increased but blood tissue fractions and speed distributions are kept constant in the dermis layers, the total perfusion will slightly decrease since the epidermis layer (with zero perfusion) will constitute a larger part of the sampling volume. See Ref. 13 for further details.

For comparison, we also compute a perfusion estimate from single-exposure contrast as

There are other single exposure models where $P$ is proportional to $K^{-2}$, e.g. Ref. 32. We have done the analysis using this model as well, but the performance was quantitatively worse than for Eq. (6), and the conclusions were unchanged by this choice. Furthermore, Eq. (6) is relatively common in clinical settings as it is used in the PeriCam PSI (Perimed AB, Järfälla, Stockholm, Sweden), and thus a good comparison for the potential benefits of our new model.

2.6.

Artificial Neural Network Perfusion Model

A neural network was trained on the data from the tissue model to estimate speed-resolved perfusion in a single pixel of the multi-exposure contrast images. Thus, the network had seven inputs representing squared multi-exposure contrast $K^2(T)$ calculated at $T=[1,2,4,8,16,32,64]$ ms exposure time and three outputs representing speed-resolved perfusion in the ranges 0 to 1, 1 to 10, and $>10\mathrm{mm}/\mathrm{s}$. The neural network was a fully connected network with a single hidden layer with 25 nodes and hyperbolic tangent (tanh) activation. Linear activation was used for the output layer. To obtain a perfusion image, the model was applied to each individual pixel of the contrast image. The architecture was determined by the process described in Sec. 2.6.1. The dataset consisted of 100,000 tissue models, divided into 70/15/15% for training, validation, and test, respectively. The validation set was used for early stopping and in-training hyperparameter tuning, and the test set was used to get an unbiased performance metric for comparison of multiple networks. 25 networks were trained from random initializations, and the best performing on the 15% test-split was selected as the final model. This model was then evaluated on a separate dataset of 100,000 tissue models, to get the final unbiased performance metrics. The training was done using the Deep Learning toolbox in MATLAB 2021b (MathWorks Inc., MA).

The model was trained using a weighted mean square error (wMSE), where the individual errors of each speed-component were weighted by their contribution to the total perfusion. Specifically,

For each speed component $j$, this metric informs the prediction error relative to the total perfusion. This accounts for the speed distribution and does not inflate small relative errors at low perfusion values, unless the total perfusion is also low. An example is warranted.

Suppose we have the true speed-resolved perfusion $P_{\text{true}}=[0.01,0.1,1]$, and the prediction $P_{\mathrm{ANN}}=P_{\text{true}}+0.01=[0.02,0.11,1.01]$. The unweighted absolute percentage error would be $\mathrm{MAPE}=100|(P_{\mathrm{ANN},j}-P_{\text{true},j})/P_{\text{true},j}|=[100,10,1]\%$, whereas the weighted absolute percentage error would be $\mathrm{wMAPE}=100|(P_{\mathrm{ANN},j}-P_{\text{true},j})/P_{\text{true},\mathrm{tot}}|\approx [0.9,0.9,0.9]\%$, better reflecting the fact that the low speed-component is small in this case and thus should not be given as much weight in the total error.

Empirical results confirmed that using wMSE for training resulted in better overall performance compared to unweighted MSE. Note that wMAPE is only an applicable metric for the speed components. For total perfusion and single-exposure perfusion, unweighted MAPE is used instead.

We also present the coefficient of determination, $R^2$, between estimated and true perfusion, which represents the proportion of the variance in perfusion explained by the model.

2.7.

In-Vivo Measurements and Measurement System

The MELSCI instrument used in this study has been thoroughly described in previous publications.^18,33 In summary, the technique is based on the continuous acquisition of images with 1 ms exposure time from a 1000 frames-per-second high-speed CMOS sensor. The images are successively accumulated to create longer synthetic exposure times at 2, 4, 8, 16, 32, and 64 ms.³³ The use of synthetic exposure times was first proposed and validated in Ref. 34. Further, the accumulation and speckle contrast algorithms are implemented in a field-programmable gate array (FPGA) closely integrated with the CMOS sensor. This enables continuous real-time acquisition and processing of multi-exposure contrast images, without any loss of information due to the negligible interframe time.¹⁸ Each set of 64 speckle images results in one multi-exposure contrast image, enabling a framerate for perfusion images of 15.6 frames-per-second.

An arterial occlusion-release provocation of the forearm was performed in a healthy 24-year-old male volunteer. The subject did not consume any caffeine or nicotine for 6 h prior to the experiment and rested for 15 min before the measurement in a 23°C temperature-controlled room. The protocol followed $5+5+5\min$ of baseline, brachial occlusion at 200 mmHg, and reperfusion, and was approved by the Regional Ethical Review Board in Linköping, Sweden (D.no. 2018/282-31). Written informed consent was given by the volunteer before the experiment.

3.1.

ANN Simulation Performance

The performance of the main and restricted ANN models on the simulated evaluation sets is presented in Fig. 3. The three speed components were evaluated using the wMAPE metric [Eq. (9)] and total perfusion was evaluated using MAPE. The coefficient of determination $R^2$ between predicted and target perfusion was also computed for all four perfusion estimates. The data was divided into 50 bins with 2000 points each, according to the true perfusion ($x$-axis). The mean predicted perfusion in each bin is shown by the black line, with the blue area representing the average deviation from the mean. The ideal prediction is indicated by the red dotted line. The top row shows the performance of the restricted model with only perfusion parameter variation and no noise, and the second row is the main model. The histograms in the last row show the distribution of true perfusion in the dataset.

Fig. 3

The top row shows perfusion predicted by the ANN trained on the restricted tissue model compared to true perfusion [%RBC × mm/s]. The second row shows the same for the main model. The data were divided into 50 bins, each with 2000 data points, based on the true perfusion ($x$-axis). The black lines are the mean predicted perfusion in the bins, and the blue areas indicate the average deviation from the mean in each bin. The red dashed line is the ideal theoretical prediction. Summary metrics weighted mean absolute percentage error (wMAPE) for speed-resolved perfusion, mean absolute percentage error (MAPE) for total perfusion, as well as $R^2$, is presented for both models. The histograms show the distribution of true perfusion in the dataset.

The second row of Fig. 3 demonstrates the combined effect on perfusion from all skin parameters in many randomized models. The relatively low $R^2$ for the lowest speed component for the main model ($R^2=0.41$; 0 to $1\mathrm{mm}/\mathrm{s}$), increased with higher speed and for total perfusion. The same is true for the restricted model presented in the first row, where all metrics are significantly better than the main model case.

To put these values in context, Fig. 4 shows similar performance plots as Fig. 3, using the main model for generating evaluation data, but for single-exposure perfusion estimated with Eq. (6). To allow comparison between the arbitrary perfusion units from the single-exposure model and the absolute units of true perfusion, the estimated perfusion was normalized to have the same mean as true perfusion in the range 0 to 2.5%RBC × mm/s. This normalization was not done for the ANN model but would only have had a small effect on the current values. Regardless, the single-exposure model performs worse than the ANN model for both metrics (MAPE and $R^2$) at all exposure times. Note that the plots in Fig. 4 compare with total perfusion and thus should only be compared to the fourth column in Fig. 3.

Fig. 4

Predicted single-exposure perfusion $P_{\mathrm{SE}}(T)$ [PU] evaluated against true total perfusion (%RBC × mm/s). To allow comparison between the arbitrary perfusion units and absolute perfusion units, the predicted values were normalized to have the same average as true perfusion in the displayed range. As in Fig. 3, the histograms show the distribution of true perfusion in the dataset.

To further evaluate how specific skin parameters influence the perfusion, one skin parameter at a time was changed over a wide range of values while the rest of the parameters were unchanged, as presented in Figs. 5 and 6. The relative change in four perfusion estimates, $P_{\mathrm{ANN},\mathrm{tot}}$, $P_{\mathrm{SE}}(1\mathrm{ms})$, $P_{\mathrm{SE}}(8\mathrm{ms})$, and $P_{\mathrm{SE}}(64\mathrm{ms})$ was compared to the relative change in true perfusion $P_{\text{true},\mathrm{tot}}$. Fig. 4(a) shows the median relative change in perfusion in 5000 tissue models due to changes in blood tissue fraction. Both $P_{\mathrm{ANN},\mathrm{tot}}$ and $P_{\mathrm{SE}}(8\mathrm{ms})$ had close to linear responses. Figure 5(b) shows similar changes when varying the mean blood speed. Here, $P_{\mathrm{ANN},\mathrm{tot}}$ had a significantly more linear response compared to all single-exposure perfusion estimates.

Fig. 5

Relative change in perfusion estimates due to a change in (a) blood tissue fraction and (b) mean blood speed. Changes are relative to perfusion estimates at (a) 0.55% and (b) 1 mm/s, respectively. Each point is the median change in 5000 tissue models randomly selected from the test dataset. An ideal method should follow the dashed line. Note the log-log scales.

Fig. 6

Relative change in perfusion estimates due to a change in (a) reduced scattering coefficient, (b) vessel diameter, and (c) epidermis thickness, indicating the largest sources of error in the absolute value of the estimated perfusion. Changes are relative to (a) $1.6{\mathrm{mm}}^{-1}$, (b) 0.055 mm, and (c) 0.205 mm, respectively. Each point is the median change in 5000 tissue models randomly selected from the test dataset. An ideal method should follow the dashed line.

Figure 6 shows the median relative change in perfusion due to changes in reduced scattering coefficient, vessel diameter, and epidermis thickness, where almost no change is expected in the ideal case (some minor changes due to changes in sample volume are expected in the ideal case). $P_{\mathrm{ANN},\mathrm{tot}}$ was highly sensitive to changes in scattering coefficient [Fig. 6(a)], although less sensitive than the single-exposure methods. While $P_{\mathrm{ANN},\mathrm{tot}}$ was less sensitive to vessel diameter than the single-exposure methods, it still deviated by up to 40% for very small vessel diameters [Fig. 6(b)]. $P_{\mathrm{ANN},\mathrm{tot}}$ was almost insensitive to changes in the epidermis thickness, whereas the single-exposure methods were increasingly sensitive with higher exposure time.

3.2.

In-Vivo Measurements

Figure 7 shows an intensity image of the forearm in the occlusion-release experiment and the region-of-interests (ROI) used in the analysis. The large black ROI was used to automatically determine the color scale for the perfusion images in Figs. 8 and 9 and for studying spatial variations in the perfusion. The small ROIs were used to extract time-resolved perfusion for Figs. 10 and 11 and were placed on high-flow (red) and low-flow (blue) regions outside any visible large vessels, to highlight the spatiotemporal heterogeneity of the microcirculation.

Fig. 7

Intensity image of the forearm presented in the occlusion-release experiment. The large black ROI was used to automatically determine the color scale for the perfusion images in Figs. 8 and 9, and for studying spatial variations in the perfusion. The red and blue ROIs were used to extract the time traces in Figs. 10 and 11 in a high-flow and low-flow region, respectively.

Fig. 8

In-vivo measurement of a forearm during an occlusion-release provocation. Columns represent the three phases of the provocation; baseline, occlusion, and reperfusion. Rows represent the three speed components and total perfusion. The color scale in each image was selected for all images to be visually comparable, using three times the mean perfusion in the large ROI in Fig. 7. Perfusion scales in the unit %RBC × mm/s are displayed in the colorbar below each image. The high perfusion area in the upper left corner of the images is an artifact due to a low intensity. Similarly, the visually high perfusion at the edges of the arm in the occlusion images is due to the low perfusion scale in combination with the low intensity.

Fig. 9

Perfusion estimates based on single-exposure LSCI using the inverse contrast model described in Eq. (6), for exposure times 1, 8, and 64 ms. The color scale in each image was selected for all images to have the same apparent mean, based on the perfusion in the large ROI in Fig. 7. Perfusion values in arbitrary units are presented in the colorbars below each image.

Fig. 10

Speed-resolved perfusion in (a) a low-flow ROI and (b) high-flow ROI during the occlusion-release provocation (baseline 0 to 5 min, occlusion 5 to 10 min, reperfusion 10 to 15 min). Care was taken to place the ROIs (see Fig. 7) outside any large vessels visible in the perfusion and intensity images. Zoom plots show 12 s during each of the three phases of the provocation.

Fig. 11

Single-exposure perfusion $P_{\mathrm{SE}}(T)$ at 1, 8, and 64 ms exposure time, for the same ROIs as Fig. 9. The data were baseline-normalized to enable comparison between the different exposure times.

Figure 8 presents speed-resolved perfusion images during an arterial occlusion-release provocation of the forearm of a healthy volunteer. The images are time averages over 2 s (32 images), selected at 2:00 min (baseline), 9:00 min (occlusion), and 10:06 min (early after release). The absolute value of perfusion is very different for the different speed components and the three phases of the provocation, making direct comparison difficult when using the same color scale for all images. To visually enhance the similarities and differences in the perfusion measures, the color scale of each image was scaled as three times the corresponding mean value computed in the region marked by the large ROI in Fig. 7 (the factor 3 was chosen empirically). The perfusion scales in the unit %RBC × mm/s are shown in the colorbars below each image. Furthermore, to quantify the heterogeneity of the microcirculation, Table 1 summarizes the coefficient of variation (CV) in the large ROI marked in Fig. 7, for each perfusion image in Fig. 8. The CV increases for the higher speed components, indicating a more spatially heterogeneous distribution of large vessels, assuming that higher speeds are more prominent in larger vessels. This is consistent with the morphology of the microvasculature.³⁵ The visible vessels in Fig. 8 are easily observed in the occlusion images for the higher speed components.

Table 1

Coefficient of variation (CV) [%] for the perfusion in Fig. 8, computed in the large ROI in Fig. 7.

Figure 9 presents images at the same time-points as in Fig. 8, with perfusion estimated by single-exposure contrast at 1, 8, and 64 ms, according to Eq. (6). Color scales were selected in the same way as described for Fig. 8. The vascular structures observed during occlusion in the ANN perfusion (see Fig. 8; 1 to $10\mathrm{mm}/\mathrm{s}$, $>10\mathrm{mm}/\mathrm{s}$ and total perfusion), are not seen in the corresponding single-exposure images (Fig. 9). The subjective heterogeneity in Fig. 9 is lower than that in Fig. 8, especially for the longer exposure times, but difficult to directly compare due to the increased noise.

Figure 10 presents time- and speed-resolved perfusion of a low-flow region (top) and a high-flow region (bottom), according to the small ROIs shown in Fig. 7. A 0.5-s mean filter was applied to the graphs in the full timescale to suppress the heart-beat variations and increase visibility of long-term changes. The zoomed-in plots show unfiltered data in 12-second sections during baseline, occlusion, and reperfusion. Gridlines in the zoom plots indicate 1-s intervals. Note that for the low-flow ROI, all speed components are of approximately the same magnitude in baseline and at reperfusion, whereas for the high-flow ROI, the high-speed component is dominant at both baseline and at reperfusion. Furthermore, for both the low-flow and the high-flow ROIs, both the 1 to $10\mathrm{mm}/\mathrm{s}$ and the $>10\mathrm{mm}/\mathrm{s}$ speed components show heart-synchronous pulsations at reperfusion.

Figure 11 presents single-exposure perfusion [Eq. (6)] in the same two ROIs, at 1, 8, and 64 ms exposure time. The perfusion at each exposure time was normalized to baseline (first 5 min) to allow for comparison between the exposure times.