Quantification of blood flow index in diffuse correlation spectroscopy using long short-term memory architecture

: Diffuse correlation spectroscopy (DCS) is a noninvasive technique that derives blood flow information from measurements of the temporal intensity fluctuations of multiply scattered light. Blood flow index (BFI) and especially its variation was demonstrated to be approximately proportional to absolute blood flow. We investigated and assessed the utility of a long short-term memory (LSTM) architecture for quantification of BFI in DCS. Phantom and in vivo experiments were established to measure normalized intensity autocorrelation function data. Improved accuracy and faster computational time were gained by the proposed LSTM architecture. The results support the notion of using proposed LSTM architecture for quantification of BFI in DCS. This approach would be especially useful for continuous real-time monitoring of blood flow.


Introduction
As an important clinical biomarker, blood flow (BF) delivers nutrients such as oxygen to tissue and removes metabolic by-products from tissue. The dynamic changes of blood flow in tissues provide useful information for diagnosis of diseases and monitoring therapeutic effects. For example, blood flow measurements in skeletal muscle hold special significance for patients with peripheral artery disease, wherein insufficient flow can lead to nutrient deficits, tissue damage, and muscle weakness [1][2][3]. Many other diseases, such as stroke [4][5][6], head trauma [7,8], cardiovascular disease and caners [9] are also associated with abnormal blood flow. Therefore, measurement of blood flow is an effective way to inform clinicians about oxygen delivery and metabolism.
Several approaches have been developed for tracking blood flow, such as computed tomography (CT) [10], positron emission tomography (PET), arterial-spin labeled MRI (ASL-MRI) [11,12], doppler ultrasound (DU) [13,14]. However, standard radiological modalities (i.e., CT, PET, and MRI) cannot provide continuous measures of blood flow. Transcranial Doppler (TCD) ultrasonography is limited to observations of large vessel flow velocities, which do not necessarily reflect microvascular perfusion in patients with cerebrovascular disease, and laser Doppler flowmetries is invasive [10]. Diffuse correlation spectroscopy (DCS) has become more attractive due to its advantages of noninvasive, no ionizing radiation, portability, inexpensive cost and potentiality of continuous bedside monitoring [15][16][17]. To date, DCS has been extensively validated to measure blood flow in vivo in deep tissues and widely applied in clinic for early diagnosis and therapeutic monitoring various diseases, such as stroke, tumor, peripheral artery disease, obstructive sleep apnea and etc.
In theory, DCS uses near-infrared light within the physiological window (∼650-950 nm) to measure blood flow, wherein photons could travel deep in tissue, as a result of the relatively small absorption of water and hemoglobin. Specifically, the temporal intensity autocorrelation function of detected light is measured, and correlation diffusion theory is employed to fit the measured curve using a tissue blood flow index (BFI) as a fitting parameter. Although previous studies have provided calibration approaches for continuous monitoring of absolute blood flow in skeletal muscle [18] and brain [19], sometimes it is difficult to achieve these schemes in current clinical practice. Thus, relative blood flow (rBF) are typically derived by the blood flow changes relative to baseline period. It has been verified in many studies that the BFI and especially its variation were demonstrated to be approximately proportional to absolute blood flow. The current approaches for quantification of BFI is to fit measured intensity autocorrelation function with either the analytical or the Monte Carlo model [20]. These approaches, however, are based on iteratively fitting which is susceptible to the data noise and quite time-consuming. To overcome these limitations of traditional fitting methods, several approaches were proposed such as Nth-order linear (NL) algorithm [21,22], least-absolute minimization (L1 norm) and support vector regression (SVR) algorithm [23]. Very recently, deep learning method was also applied for quantification of BFI on a human subject [24]. For general-purpose sequence modeling, long short-term memory (LSTM) [25,26], as a typical RNN structure, has proven stable and powerful for modeling time sequences in various previous studies [26,27], such as detection of epilepsy seizures [28], emotion recognition [29] based on EEG signals.
Here, we investigated and assessed the utility of proposed LSTM architecture for quantification of blood flow index in DCS. Phantom and in vivo experiments were established to measure normalized intensity autocorrelation function g 2 (τ) data. To verify our hypothesis, the g 2 (τ) datasets in both phantom (n=24000) and in vivo (n=67200) experiments were trained to develop LSTM models for BFI quantification. Performances of the LSTM models were calculated and compared with traditional fitting method using prediction groups.

Methods
We present a LSTM-based deep learning model as an alternative method to extract the blood flow index comparing to traditional fitting method in diffuse correlation spectroscopy technique. Phantom experiment and in vivo experiment were established in this study to verify the proposed LSTM architecture.

Diffuse correlation spectroscopy
DCS estimates blood flow by quantifying rapid speckle intensity fluctuations of multiply scattered coherent near-infrared light induced by red blood cell motion. These speckle intensity fluctuations are characterized by computing the normalized intensity autocorrelation function, i.e., at multiple delay-times, τ, where I(t) is the detected light intensity at time t, and the brackets <> represent time-averages (for experiments) or ensemble averages (for calculations). To quantitatively relate the measurements of g 2 (τ) to blood flow (BF), a correlation diffusion approach is employed to calculate the electric field autocorrelation function, i.e., as a function of blood flow index (BFI) describing the dynamics of red blood cells. The blood flow index is ascertained by fitting the calculated g 1 (τ) to the experimentally measured g 2 (τ) using the Siegert relation [15,30]: Here, β is a constant, y-intercept on the g 2 axis, depends on the number of spatial modes detected. Hence, the aim is to minimize the sum of squared differences (SSD) between the calculated and measured autocorrelation functions g 2 (τ). The minimization of SSD is performed by using Nelder-Mead simplex algorithm [31] (fminsearch function) in MATLAB (Mathwork, Inc., USA) [32]. Although the BFI does not have absolute blood flow units, it is directly proportional to tissue blood flow [18,19], i.e., BF = γBFI (4) and has been successfully validated against a plethora of gold-standard techniques, such as ASL-MRI [11,12], transcranial Doppler ultrasound [14] and so on [10]. Here, γ is a proportionality constant of which the unit is [(ml · 100ml −1 · min −1 )/(cm 2 /s)]. It is apparent that relative changes of tissue blood flow (rBF) compared to a baseline blood flow (BF 0 ) can be derived from corresponding changes in BFI, i.e., The proportional relationship between BFI and blood flow has been introduced by several review papers [33,34] based on many validation studies of DCS blood flow indices in human brain and muscle.

Instrumentation
The DCS system employed in this study was described in our previous study [35]. The system consists of a continuous-wave laser (785 nm, CrystaLaser Inc., USA) with a coherence length longer than 10 meters, a single-photon-counting module (APD) (SPCM-AQRH, Excelitas Inc., Canada), and real-time software correlator based on PC-based data acquisition board (PCIe-6612, National Instruments Inc., USA) [36]. The laser source illuminates the sample via a multimode fiber (200 um diameter, Chunhui Science & Technology Inc., China), and remitted light that travelled through the sample is detected by a single mode fiber (5 um diameter, Chunhui Science & Technology Inc., China) 1.5 cm away from the source. The single-photon-counting module (APD) converts the detected photons to digital TTL pulses, which are sent to pulse counters in the data acquisition board. The counters continuously accumulate photon counts, and the computation of the normalized intensity autocorrelation function g 2 (τ) is performed by software implementing a multi-tau scheme. Traditionally, correlators employ embedded programming and a multi-tau algorithm to compute the autocorrelation functions over a large range of delay times (from ∼ 1µs up to as much as 1 ∼ 2s). Briefly, as we know that deep tissue g 2 (τ) curve most often decays to 1 at delay times of τ ∼ 250 µs or less. Hence, correlation data at delay times greater than ∼ 250 µs does not offer significant information about tissue dynamics. Based on the previous study [36], 64 delay-times were set in this study, which were ranged from 3.3333 µs to 3.4 ms at logarithmic rule. The number of delay-times were reduced comparing with the traditional hardware correlators, hence the data acquisition can be made more efficient and faster. The frequency of g 2 (τ) measurement was set as 10 Hz in both phantom and in vivo experiment.

Phantom experiment
Liquid phantom serves as the gold standard to test our algorithm using a combination of absorbers and scatterers with known spectra. The liquid phantom, contained in a tank (200 mm ×120 mm × 150 mm), was comprised of distilled water, India ink (Royal Talens Inc., Netherlands) and 30% intralipid (Huarui Pharmaceutical Inc., China). The expected µ a is 0.05 cm −1 , and µ' s is 8 cm −1 . A glass cuboid (100 mm × 30 mm × 10 mm) filled with liquid phantom was placed in the tank as shown in Fig. 1. The import of the glass cuboid was connected to a 50 ml syringe by a plastic flexible tube (3 mm in diameter). The syringe filled with liquid phantom was placed in a medical injection pump (SP-2000, Ningbo Annol Medical Device Technology Inc., China), which was employed to control the velocity of liquid phantom in the glass cuboid. The liquid phantom flowing in the glass cuboid was directly output to the tank through the export of the glass cuboid. The DCS optical probe was fixed on the surface of the glass cuboid. After a 1-minute baseline, the velocity of medical injection pump was set to 400 ml/h and maintained for 1.5 minutes. Then, the velocity increased 0.5 ml/h and maintained for another 1.5 minutes. In this way, the velocity increased 50 ml/h per 1.5 minutes until it increased to 6 ml/h. Thus, 5 velocities (i.e., 400 ml/h, 450 ml/h, 500 ml/h, 550 ml/h and 600 ml/h) of the liquid phantom in total were measured in one trial as shown in Fig. 2(a). Twelve trials were performed in phantom experiment. Please note that when we set the velocity of the medical injection pump, the velocity of the liquid phantom in the glass cuboid could not immediately increase to the setting velocity. To ensure that the accuracy of liquid phantom velocity, the 30 seconds data measured after the velocity setting point were excluded. Hence, 1-min data were left for each velocity in one trial, which means that 5-min data in total were used for the following processing in each trial. Example blood flow measurement acquired from the liquid phantom during a single trial is shown in Fig. 2(b). The BFI values, which increased with the velocity, were extracted based on the traditional fitting method. In addition, example g 2 (τ) curves measured from the liquid phantom to quantify BFIs were presented in Fig. 2(c). These example g 2 (τ) curves were selected randomly from five different velocities and baseline period (i.e., the pump velocity is zero), respectively. The g 2 (τ) curves at five different velocities are selected randomly from the example BFI in panel (b). The g 2 (τ) curve in the baseline period marked as yellow pentagram in Fig. 2(c) is selected randomly from 1-min baseline measurement (i.e., the pump velocity is zero) of the example BFI. It can be seen that the g 2 (τ) curve in the baseline period decayed slowly, which is similar as that from the hardware correlator in previous phantom experiment [37]. As the pump velocity increased, the g 2 (τ) curves at five different velocities decayed much faster, which were marked as circular, diamond, triangle, square and star in Fig. 2(c). It was shown that the g 2 (τ) curves decayed faster as the pump velocity increased. Previous studies [38] [39] have demonstrated that the pump speed is proportional to the blood flow measured by DCS. BFIs quantified by LSTM model and traditional fitting method were compared with the known pump velocities of liquid phantom.

In vivo experiment
We also established in vivo experiment with cuff occlusion. Twelve healthy adult volunteers (seven males and five females, age 23 ± 1.6, BMI 22.3 ± 4.2) participated in this study after giving their written informed consent. The subject sat comfortably with his/her right arm extended in a comfortable position throughout the arterial occlusion (AO) protocol (note that all subjects measured in the study were right-handed, hence the right arm was their dominant arm). The DCS optical probe was secured to the subject's right forearm with an elastic neoprene strap indicated. A blood pressure cuff was placed around the right bicep. After a 3-minute baseline, the cuff was rapidly inflated to 200 mmHg to achieve arterial occlusion. Arterial occlusion was maintained for 1 minute. The cuff was then rapidly deflated, and the recovery was monitored for 3 minutes.
The 7-min "AO trial" protocol, as shown in Fig. 3(a), was repeated twice for each subject. Example blood flow measurement acquired from a subject during a single "AO trial" is shown in Fig. 3(b). Similarly, the BFI values were extracted based on the traditional fitting method. In addition, example g 2 (τ) curves measured at four different time points (i.e., baseline, AO, overshoot and recovery) from this subject were presented in Fig. 3(c), respectively. The time points selected for the g 2 (τ) curves were marked as circular, star, triangle and square in Fig. 3(b), respectively.

Dataset
The blood flow index traditionally is ascertained by fitting the calculated g 1 (τ) to the experimentally measured g 2 (τ) using the Siegert relation. The dataset consisted entirely of the measured g 2 (τ) data. Of note, the measured g 2 (τ) curve was computed by software implementing a multi-tau scheme of which contained 64 delay-times. Hence, one g 2 (τ) curve contained 64 data points (i.e., size of g 2 (τ) size: 1×64).
For phantom experiment, 12 trials were performed in total. The phantom experiment data set was divided into training group (8 trials), verification group (2 trials) and prediction group (2 trials) randomly. The effective measurement duration of one trial was 5 minutes (1 minute per one velocity, the first 30 seconds data measured after the velocity setting point were excluded). As mentioned before, the blood flow index data acquisition frequency was 10 Hz. Thus, the amount of g 2 (τ) measured for one trial was 3000. Each g 2 (τ) included 64 data points. The size of g 2 (τ) for training group (8 trials) was 24000×64, the size of g 2 (τ) for verification group (2 trials) and prediction group (2 trials) were both 6000×64.
For in vivo experiment, 12 subjects were participated in this study and each subject performed twice 7-min "AO trial". Hence, 24 trials of 12 subjects were established in total. Similarly, the in vivo experiment data set was also divided into training group (16 trials), verification group (4 trials) and prediction group (4 trials) randomly. The 2 trials of one subject were in the same group. The arterial occlusion (AO) protocol lasted for 7 minutes for a single trial (4200 data points). The blood flow index data acquisition frequency was also 10 Hz. Thus, the size of g 2 (τ) for training group was 67200×64, the size of g 2 (τ) for verification group and prediction group were both 16800×64.
In brief, effective measurement duration of phantom experiment was 60 minutes in total, of which the g 2 (τ) data across five velocities (i.e., 400 ml/h, 450 ml/h, 500 ml/h, 550 ml/h and 600 ml/h). In other words, the size of total g 2 (τ) data in phantom experiment was 36000×64, and the ratio of training group, verification group and prediction group was 4:1:1. The measurement duration of in vivo experiment was 168 minutes in total, of which the g 2 (τ) data across four different time points (i.e., baseline, AO, overshoot and recovery). The size of total g 2 (τ) data of in vivo experiment was 100800×64, and the ratio of training group, verification group and prediction group was also 4:1:1.

LSTM architecture
For general-purpose sequence modeling, LSTM as a typical RNN structure has proven stable and powerful for modeling time sequences in various previous studies. The LSTM architecture varies with the change in the nature of data. The number of memory units in the hidden layers, the activation function and number of stacked layers are variables in any LSTM structure. The proposed LSTM architecture in this study is shown in Fig. 4. The model input is the measured g 2 (τ) data of which the size is 1×64. The model output is the predicted blood flow index data BFI of which the size is 1×1. Note that, the measured g 2 (τ) data were normalized before putting into the model. Specifically, as shown in Fig. 4, input (1×64) is divided into 8 groups, each group is a 1D signal (1×8). These 1D signals are fed to two LSTM layers one by one, i.e., for totally 8 times. Then, the Rectified Linear Unit (ReLU) 38 operation is applied to introduce nonlinearity. At the end of processing, the fully connected layer uses the softmax classifier 46,47 to output the predicted blood flow index data. Table 1 lists all the selected parameters of the proposed LSTM architecture. The mean square error (MSE) loss function was employed as the loss function. Adaptive moment estimation (Adam) 48 was adopted for the training process, which iteratively updates the neural network weights based on the training data. The learning rate was initialized with 1e-4. The model was implemented in python using Pytorch with Intel Core i7-4700MQ CPU.

Evaluation criteria
The normalized intensity autocorrelation functions g 2 (τ), measured in phantom and in vivo experiments, were processed based on both the traditional fitting algorithm and proposed LSTM model to quantify BFIs. Usually, root mean squared error (RMSE) and mean absolute error (MAE) are the common used evaluation criteria to assess the quantification performance. In this study, however, RMSE percentage (RMSE%) and MAE percentage (MAE%) were calculated, as follows: Here, BFI pred represents the BFI values predicted or measured, BFI true represents the BFI true values. The advantage of RMSE% and MAE% is that average of the true values are contained in their definitions. In this study, for the phantom experiments, five velocities of medical injection pump were defined as true values of BFI (i.e., BFI true ), since they were already known. For the in vivo experiments, the BFIs quantified by the traditional fitting method were defined as true values of BFI (i.e., BFI true ).

Results
Results from phantom experiment and in vivo experiment were presented in this section. For phantom experiment, BFIs quantified by LSTM model and traditional fitting method were compared with the known velocities of liquid phantom. For in vivo experiment, BFIs quantified by the proposed LSTM model were compared with BFIs quantified by traditional fitting method.

Phantom experiment results
Outputs of the phantom LSTM model includes 2 trials BFI data based on the measured g 2 (τ) data in prediction group. The measured g 2 (τ) data of these 2 trials were also used to estimate BFI by traditional fitting method. In order to compare with the known velocities of phantom, relative blood flow changes (rBF) were derived from corresponding changes in BFI based on Eq. (5). Here, the baseline blood flow index (BFI 0 ) is defined as the mean value of their first 1-min data. Similarly, relative changes of the phantom velocities varied by medical injection pump were also normalized to the mean value of its first 1-min data. The above three parameters were denoted as rBF LSTM , rBF trad and rBF v , respectively.
Here, BFI LSTM is the BFI calculated by the LSTM model, BFI trad is the BFI calculated by the traditional fitting method, BFI v is the BFI calculated based on the velocity of the medical injection pump. Previous studies have demonstrated that the pump speed is proportional to the blood flow measured by DCS [38,39]. Thus, the relationship between BFI and velocity of the medical injection pump can be defined as BFI = kV. Here, k is the proportionality constant, and V is velocity of the medical injection pump. Note that, subscript 0 indicated the mean value of its first 1-min data. Figure 5 displays the time traces of rBF (3000 data points) during a single trial across five velocities. Figure. 5(a) shows the time traces of rBF LSTM (solid blue line) and rBF v (solid red line), while Fig. 5(b) shows the time traces of rBF trad (solid blue line) and rBF v (solid red line). More careful observation of the data, the 0.3 min extracted time-windows present in Fig. 5(c) and Fig. 5(d). To contrast with rBF v , it can be seen that the rBF LSTM and rBF trad both match well with the known phantom velocities. However, the traditional fitting method causes larger fluctuations than LSTM model. The Pearsons correlation coefficient between rBF LSTM and rBF v is r 2 =0.95 while the Pearsons correlation coefficient between rBF trad and rBF v is r 2 =0.74.  Figure 6 displays RMSE% and MAE% calculated on five phantom velocities. No significant differences were found to be related with velocities. The RMSE% of LSTM model and traditional fitting method ranges from 0.99% to 3.95% and 6.04% to 8.29%, respectively. The MAE % of LSTM model and traditional fitting method ranges from 0.05% to 3.37% and 5.1% to 6.6%, respectively. RMSE% calculated on all the phantom velocities were 3.3% for LSTM and 7.62% for traditional fitting method. MAE% calculated on all the phantom velocities were 2.18% for LSTM and 5.9% for traditional fitting method. Better performance was observed for LSTM model.

In-vivo experiment results
Outputs of the in vivo LSTM model includes 4 trials BFIs data. Figure 7(a) shows the time traces of BFIs (4200 data points) during an example AO trial quantified by LSTM model (BFI LSTM ) and traditional fitting method (BFI trad ). It demonstrates that the two methods had same trend and good correlation in BFIs. Figure 7(b) shows a 0.5 min extract of BFIs in baseline, i.e., yellow shaded region in Fig. 7(a). It can be seen that pulsatile nature of tissue blood flow were recorded with clear resolution. Moreover, RMSE% and MAE% of quantified BFI for this AO trial were 4.39% and 3.38%.    Fig. 8(a); solid red regression line, dotted green one-to-one line] and a Tukey mean-difference (or Bland-Altman) analysis [ Fig. 8(b)] show excellent agreement between the two methods, per extraction of BFIs. The slope between BFI LSTM and BFI trad is 0.97 and the Pearsons correlation coefficient is r 2 =0.99. The observed agreement between the BFIs quantified by the two methods in Fig. 8 further supports the notion of using LSTM model for quantification of BFI in DCS.

Discussion
Our study demonstrates the promise of using LSTM architecture to quantify BFI for blood flow measurement in software-correlator-based DCS. This approach can readily be used for continuous blood flow monitoring in long time, such as all-day monitoring of stroke patients.
LSTM architecture used in this study has been widely applied in various areas of medicine, including early diagnosis of Alzheimer's disease [40,41], ECG signals for detection of epilepsy seizures and emotion recognition. We used LSTM architecture for quantification of BFI from the measured temporal intensity autocorrelation functions. This proposed approach was inspired by the fact that the LSTM is stable and powerful for modeling time sequences, such as ECG data. Our experimental results indicate the prediction performance of the LSTM model exceeded the performance of the traditional fitting method (i.e., Nelder-Mead simplex algorithm). This suggests that the LSTM architecture could be an alternative approach for quantification of BFI in DCS, especially for continuous long-time monitoring of blood flow.
The phantom experiment was established to verify feasibility of the proposed LSTM model in BFI quantification, and compare its results with the traditional fitting method. The proportional relationship between pump velocity and blood flow were validated [38,39]. The recent study [38] conducted the phantom experiment by changing the pump velocity as this study. Linear regression between the pump velocity and BFI was presented with excellent correlations (r 2 =0.998). Different shapes of the anomaly and manipulated pump speed were designed to mimic the various situations inside biological tissues [39]. Three different pump speeds were varied in their phantom experiment and it was observed that the BFI values were positively relevant to the pump speed in their study. It is clear that the pump velocity is proportional to the BFI measured by DCS. In our phantom experiment, BFI of the liquid phantom were known which is directly proportional to velocities controlled by the medical injection pump. Therefore, the LSTM model was trained and verified to the known velocities based on the training and verification group.
As shown in Fig. 5, the time traces of rBF during a single trial across five velocities matched well with the known phantom velocities. The flow of the phantom caused by the medical injection pump was fast and unstable, hence the Fig. 5 demonstrates more fluctuations. Moreover, BFI was measured based on the software correlator and the acquisition frequency was set as 10 Hz, which was at least 10 times faster than the original hardware correlator. The acquisition frequency of BFI based traditional hardware correlator is usually less than 1 Hz. Much more data points (at least 10 times) were gathered in the same time period, thus more fluctuations seemed to be observed. Actually, these fluctuations were the flow of phantom itself. The Pearsons correlation coefficient between the LSTM model and the known velocities was r 2 =0.95, while the RMSE% and MAE% were 3.3% and 2.18%. By constrast, these parameters of traditional fitting method were r 2 =0.74, RMSE%=7.62%, and MAE%=5.9%. This observation suggests that the accuracy of BFI quantification was improved by the proposed LSTM model for measuring blood flow.
For in vivo experiment, we actually could not obtain the true BFI of subjects. Thus, the BFIs estimated by the traditional fitting method were considered as the standard. The LSTM model was trained and verified to the BFIs estimated by the traditional fitting method. We observed a high correlation (r 2 =0.99), good agreement (slope=0.97) and small difference (RMSE%=4.39%, MAE% = 3.38%) in BFI predicted by LSTM model versus the traditional fitting method in Fig. 8. This observation suggests that the proposed LSTM model as an alternative approach for BFI quantification for in vivo blood flow monitoring.
In addition, the computing time of BFI quantification is important for real-time blood flow monitoring, especially for continuous long-time monitoring. We measured temporal intensity autocorrelation functions at 10 Hz based on software correlator in both phantom and in vivo experiment. The in vivo experiment, for instance, took 7 min to measure 4200 g 2 (τ) curves during an AO trial. The proposed LSTM model took 44.88 ms to quantify the 4200 BFI data points vs the traditional fitting method took 2381.86 ms. However, if we only quantified one BFI data, the proposed LSTM model took 19.95 ms vs the traditional fitting method took 176.16 ms. We found that the bigger the amount of BFI data was, the more time would be saved by LSTM model. The computing time (including the time of transferring data and inference time) and the inference time of BFI quantification based on these two approaches were provided in Table 2. As shown in Table 2, the computing time of 4200 BFI data points was less than the twice the computing time of 2100 BFI data points. The inference time of 4200 BFI data points was also less than the twice inference time of 2100 BFI data points. The time of transferring data for LSTM approach was around 12.56 ± 1.52 ms, while the time of transferring data for traditional fitting method was around 162.45 ± 42.32 ms. Thus, the computing time of BFI quantification is mainly determined by the inference time. The inference times per sample of the proposed LSTM approach and traditional fitting method were 0.99 ms and 2.48 ms, respectively. More than half of the time was saved by the proposed LSTM approach. Especially for continuous long-time monitoring, for example, for continuous blood flow measurement of 24 hours, the quantification of BFI (4.32×10 6 data points) would take less than 46.16 s by the proposed LSTM model. Hence, the proposed LSTM model performed good computational capacity in BFI quantification for blood flow monitoring. Therefore, it is most probably to realize the real-time blood flow monitoring with the LSTM architecture. LSTM and the graded recurrent unit (GRU) are two variants of recurrent neural network (RNN) [42,43]. Since GRU normally needs less parameter compared with LSTM, we compared the performance of LSTM model with GRU model using the same in vivo dataset. Table 3 was the evaluation parameters of the same example AO trial shown in Fig. 7. As shown in Table 3, although GRU model had slightly fewer parameters, all the evaluation parameters, including r 2 , RMSE% and MAE%, demonstrate that the LSTM model has better performance than the GRU model. The source-detector separation of 1.5 cm was used in this study for both phantom and in vivo experiments. Normally, 1.5 cm is sufficient to probe the muscle tissue, such as skeletal muscle in forearm or calf. Based on photon diffusion theory in biological tissues, light penetration depth depends on tissue optical properties and the source-detector separation. The penetration depth of NIR light in biological tissues is approximately half of the source-detector separation. The thickness of the superficial layers (skin and adipose tissues) above forearms was 2.8 ± 0.6 mm, which was measured from n=10 volunteers in previous study [44]. Besides, the source-detector separation of 1.5 cm was established in several muscle blood flow measurement studies [45][46][47]. For instance, blood flow in the forearm during artery cuff occlusion was measured using 1.5 cm source-detector (S-D) separation between the two laser diodes and the camera in Ref. [45]. Blood flow in the forearm during handgrip exercise was measured using 1.5 cm source-detector (S-D) while the sampling rate was 1 Hz [46]. Four source-detector separations (10.0, 15.7, 22.8, and 30.0 mm) were used to measure muscle blood flow during the operation at intraoperative flap elevation step and on the three postoperative days [47]. In our in vivo study, the blood flow was also measured in the forearm. Thus, the DCS signals from the 1.5 cm separations derive mainly from the muscle layer, which is sufficient to probe the muscle tissue. Our future work will further validate the performance of our proposed algorithm by measuring cerebral blood flow using the DCS probe with larger source-detector separation (2.5 cm or above).

Conclusion
In this study, we proposed a LSTM architecture for quantification of BFI from the measured normalized intensity autocorrelation function g 2 (τ) in diffuse correlation spectroscopy. The experimental results indicated that our proposed LSTM model observed improved accuracy and faster computational time than traditional fitting method. This approach would be especially useful for continuous real-time monitoring of blood flow. Data availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.