3.1.

Data Acquisition Principle

The performance was controlled by LabVIEW (National Instruments), which is one of the most popular and powerful program tools and suitable for interactive data acquisition and signal processing. The TTL output of the detector was counted over sampling time $\tau$ by the counter/timer board. The PCI-6602 has a timebase of 80 MHz, which means the fastest input pulses that can be counted are 12.5 ns ($1/80,000,000\mathrm{s}$) apart. The PCI-6602 is the 32-bit counter and can count up to 4,294,967,295 ($2^{32}-1$) before it rolls over. The counter can count the events continuously and store the data in the buffer. Due to the limitation of buffer size, time resolution of anything less than 2 μs would fill up the buffer in less than 1 s. For the purpose of stable data acquisition, we conducted experiments using 2.5 μs. In the data acquisition program, we applied “queue” structure to link the “producer” loop and the “consumer” loop, which were in charge of reading photon counts and writing results to a file, respectively. By using this kind of structure, the two loops were able to work in parallel so that a faster speed can be achieved. Simultaneously, another LabVIEW program detected the latest written photon-counting file and calculated the autocorrelation function by FFT algorithm. The result was then displayed and saved to a specified folder. Running on our current system, one core of CPU was fully occupied to acquire the data, while the other core was used for the autocorrelation calculation.

3.2.

Data Processing

In this section, the working principle of the hardware correlator board will be briefly described, based on the instruction manual of a most commonly used commercial product (Flex99R-12) from Correlator.com. The algorithm behind the hardware correlator will be simulated by the software approach, and its performance will be compared against our new FFT-based approach in a later section.

3.2.1.

Simulating hardware autocorrelator board

Inside the hardware autocorrelator board

The typical hardware correlators such as the ones manufactured by ALV and Correlator.com are based on a multitau scheme, and have a multiple tier structure which contains a certain amount of registers as shown in Fig. 1. If we take the one from Correlator.com for instance, the first tier has 32 registers, and higher tiers have 16 registers each. The incoming photon count starts filling up the registers in the first tier from the left with a given delay time $T_0$ (typically 160 ns). When it reaches the right end, the second tier fills up from the left at a twice as slow rate, since the sum of two register values in the right end of the first tier constructs the first register value of the second tier. Likewise, the first register value of nth tier is constructed by sum of the last two register values in ($n-1$)th tier, which means the registers of $n$’th tier is updated every $2^{n-1}{T}_0$. The hardware correlator can have upto 31 tiers, which means there are $512[=(1\times 32)+(30\times 16)]$ registers with 512 different delay times. However, registers at the 31st tier updates every $2^{30}{T}_0$, which is too long with a practical $T_0$ value of $\sim 100\mathrm{ns}$, so later tiers are hardly used in practice.

Fig. 1

The diagram of hardware autocorrelator inner structure, with several tiers of registers where the delay times are the power of two times of initial delay time $T_0$.

For each shift that occurs in the register, the unnormalized intensity autocorrelation coefficient for the ith register, $G_2(t_i)$, is calculated as the average value of $n_i·n_0$, where $n_i$ indicates the photon count in the ith register, $n_0$ is the photon count at zero delay time with the same bin width as the ith register, and $t_i$ is the delay time between $n_i$ and $n_0$.²⁸ The delay time $t_i$ is calculated as the summation of all the bin widths on the left of the ith register. With a few hundred register channels (512 channels for a 31 tier system, for example), a delay time that ranges from hundred of nanoseconds to minutes can be achieved, greatly reducing the computational load compared to the linear autocorrelator.

Software implementation

The above mentioned hardware correlator algorithm is implemented in Matlab, in order to examine its behavior with a known input signal, generated either by simulation or experiment. The number of tiers, $N$, was given as a variable (normally set at around 10), and a cell array variable with a size of $16(N+1)$ [32 for 1st tier, and $16(N-1)$ afterwards] was allocated for register channels. The Matlab script initiates a loop, where the register channels in the first tier shifts to the right by one register per iteration and a new datum is read in the left-most register. When the iteration number reaches $2^{n+1}$, registers in $n$’th tier become shifted to the right by one register, after which a new value [the sum of two register values at the end of ($n-1$)’th tier] is assigned to the left-most register in $n$’th tier. Figure 2 shows a snapshot of the register profile in different tiers, where the register values are represented in a logarithmic scale to better visualize the wide-range of the values.

Fig. 2

A snapshot of the registers of simulated hardware correlator. Color represents the logarithm of the photon count in each register. Ignore the blue portion on the left where there are no registers.

When the register fully fills up, the autocorrelation coefficient between ith register and the 0th register (with the same bin width, as described previously in the hardware correlator section) is calculated. The autocorrelation curve for whole register is calculated for each shift event in the last tier, and they are averaged over hundreds of trials before being plotted.

3.2.2.

FFT-based software autocorrelator

The unnormalized intensity autocorrelation, represented as $G_2(\tau )=\langle I(t)I(t+\tau )\rangle$, resembles the convolution operation, except that in the convolution calculation, the second term being multiplied has to be the time-reversal of the original function $I(t)$, namely $I(-t)$. Therefore, the same autocorrelation can be obtained by convolution between two functions $I(t)$ and $I(-t)$. Then, by the convolution theorem, one can obtain the autocorrelation function based on the Fourier transform as follows:

Eq. (4)

$$G_2(\tau )=(I(t)*I(-t))(\tau )=F^{-1}[F\{I(t)\}F\{I(-t)\}]=F^{-1}\{\tilde{I}(\upsilon )\tilde{I}(\upsilon {}^{*})\}=F^{-1}\{{|\tilde{I}(\upsilon )|}^2\},$$

where $F$ and $F^{-1}$ refer to forward and inverse Fourier transform, respectively, and $\tilde{I}(\upsilon )$ is the Fourier transform of $I(t)$. Note that the last equality of Eq. (4) is also called Wiener-Khinchin theorem and real valued $I(t)$ always results in real values of $G_2$.

Since the Fourier transform can be very efficiently implemented by the fast Fourier transform (FFT), one can obtain the linear autocorrelation function in a very fast and stable manner using above equation.

Figure 3 shows the comparison between the two normalized intensity autocorrelation $g_2(\tau )$ calculations, one from the multi-tau approach simulating the hardware correlator board, and the other from the FFT-based linear approach. The data used in Fig. 3(a) was from a flow phantom experiment, acquired by the NI counter board (PCI-6602) at the sampling rate of 400 kHz. For the FFT-based approach, the data size used was about $4\times {10}^5$, covering about 1 s. For the hardware simulating approach, the delay time was only calculated upto about 30 ms due to the small number of tiers (10 tiers were used), but it was averaged over 500 autocorrelation calculations, and the total length of data used is about the same. The FFT-based calculation always ends up showing a spike near the highest delay time, but that is due to the cyclic nature of the FFT-based result, and our visualization in Fig. 3 only shows the first half, thereby removing the unwanted peak. Both of them show a nice convergence to the value of 1 for long delay times, and they nicely overlap with each other. For the experimental data in general, we were able to observe a trend that the hardware (multitau) algorithm tends to show higher fluctuation for small delay times, whereas the FFT-based algorithm always show the smoothest behavior for small delay times due to its single-tau calculation.

Fig. 3

Comparison between hardware and software autocorrelation calculations. (a) Autocorrelation calculated from a typical flow phantom data. (b) Autocorrelation calculated from a simulated data. Hardware-simulating algorithm (red squares) was implemented under the assumption of a 10-tier register, averaged over 500 different autocorrelations. The sampling rate was 400 kHz.

The comparison was also performed using simulated data based on a Markov chain with a Gaussian noise, and the result was always overlapping unless the average count value of the data was too low. When the average was low, we were able to regain the equivalence between the two methods by binning the data, thereby increasing the average photon count value. Typical autocorrelation curves for simulated data are shown in Fig. 3(b).

4.1.

In-House Flow Phantom Experiment

Phantom, which is the tissue-simulating object to mimic the properties of human or animal tissues, has played an important role in the evolvement of diagnostic systems, and most physical therapeutic interventions. It can be used for the purposes of initial testing of system design, optimizing signal to noise ratio (SNR) in existing system, and so on.⁴⁴ In order to validate and assess the ability of our system, a flow phantom was built as shown in Fig. 4.

Fig. 4

The in-house silicon flow phantom, in which an acrylic cylindrical tube was embedded to mimic the blood vessel. The tube was filled up with glass beads with a diameter of 3 mm in order to make random flow directions.

The phantom making procedure can be found elsewhere.⁴⁵ Briefly, RTV 12A was the base compound, RTV 12C (Momentive Performance Materials) was the curing agent, while carbon black and ${\mathrm{TiO}}_2$ were mixed to adjust the absorption coefficient and scattering coefficient, respectively. The optical properties of the solid flow phantom were ${\mu }_s^{\prime }=8{\mathrm{cm}}^{-1}$ and ${\mu }_a=0.03{\mathrm{cm}}^{-1}$ (measured by OxiplexTS, ISSInc, at 830 nm), which satisfied the ${\mu }_s^{\prime }\gg {\mu }_a$ criterion, and the photon diffusion equation remains valid.

The schematic of the DCS experimental setup is shown in Fig. 5. The Lipofundin N 20% (B.BraunMelsungen AG, Germany) with a concentration of 0.6% was pumped through the cylindrical tube, with flow speeds of $0\mathrm{ml}/\mathrm{s}$, $0.05\mathrm{ml}/\mathrm{s}$, $0.125\mathrm{ml}/\mathrm{s}$, $0.2\mathrm{ml}/\mathrm{s}$, and $0.325\mathrm{ml}/\mathrm{s}$. The source-detector pair was centered over the solid phantom surface at a distance of 1 cm from the flow tube surface. Reflection geometry was adopted in which incident light was injected into the phantom by the source fiber and detected 3.0 cm away along with the cylindrical tube direction.

Fig. 5

(a) Schematic of the DCS experimental setup using in-house flow phantom. (b) Screen shot of the LabVIEW interface of the FFT-based autocorrelation calculation.

4.2.

In Vivo Cuff Occlusion Experiment on a Human Arm

After the flow phantom experiment, we conducted the validation study on a human subject. The in vivo cuff occlusion experiment was done on the arm of a 28 year-old healthy male subject. In order to vary the blood flow speed, a cuff on the subject’s upper arm was constricted for 15 s for temporary arterial occlusion. The source-detector separation was 1.5 cm and the probes were located on the forearm. The cuff inflation pressure was 200 mmHg.

4.3.

Photodynamic Therapy Response Monitoring

Besides the flow phantom and in vivo human arm cuff occlusion experiments, we also applied our instrument to a lab-based PDT experiment for evaluation. Since the antivascular effects of PDT is a key indicator of tissue response to treatment,²⁹ we aim to provide instantaneous hemodynamic feedback on treatment response. The in vivo PDT experiment was approved by the Institutional Animal Care and Use Committee of Singapore Health Services Pte Ltd and carried out at the National Cancer Centre Singapore. A Balb/c mouse bearing xenograft tumors on the flanks was used as an animal model. The mouse received chlorin e6 (Ce6) as the photosensitizer at a dose of $5\mathrm{mg}/\mathrm{kg}$. Approximately 4 h after injection, the tumors were irradiated with laser light. The irradiation lasted 2450 s, consisting of alternating 250-s light exposures and 300-s dark intervals. The treatment light for PDT was delivered by a 665-nm medical-grade laser (Biolitec, Germany). The light dose used was $100\mathrm{J}/{\mathrm{cm}}^2$, delivered at a rate of $80\mathrm{mW}/{\mathrm{cm}}^2$ via an optical fiber. The treatment light was expanded to a 2.5-cm diameter on the tumor surface. Our DCS measurement was carried out before and immediately after each exposure, as well as at 25 and 50 min after the last treatment. Source-detector separation was 1 cm, which fitted the size of the tumors.

5.1.

In-House Flow Phantom Experiment

Figure 6 shows the normalized intensity autocorrelation functions in a semilogarithm plot with varying flow speed. It can be observed that the normalized intensity autocorrelation function $g_2(\tau )$ showed an initial constant plateau, but decayed exponentially with increasing delay time. As speed increased, the decay rate of the autocorrelation curve increased. If we recall the definition of autocorrelation, the finite data length will always cause the copy of the data shift out of the data window, thus eventually the autocorrelation function will go down to zero naturally. However, by employing an FFT, one can actually get the cyclic autocorrelation, thus it avoids the edge effect caused by data finiteness and offers a faster calculating speed. Ideally, the zero flow ($0\mathrm{ml}/\mathrm{s}$) was not supposed to show any decay, but we observed a decay that is mainly attributed to the intrinsic fluctuation of the laser source and the room temperature Brownian motion of the Lipofund in remaining in the tubing. The autocorrelation function calculating speed was less than 20 milliseconds for 800,000 data points with the type of double-precision integer, which is faster than the data writing speed (new file generating speed) and enables the display of the autocorrelation curve in almost real-time. We note from the autocorrelation curves, that increasing the flow causes the plateau region of the curves to be shorter and makes it difficult to determine the $\beta$ value. Therefore, it might cause inaccuracy in measuring high flow, but one can obtain stable plateaus for most of the relevant flow ranges in vivo.

Fig. 6

Autocorrelation curves of the phantom with varied speed of flow ($0\mathrm{ml}/\mathrm{s}$, $0.05\mathrm{ml}/\mathrm{s}$, $0.125\mathrm{ml}/\mathrm{s}$, $0.2\mathrm{ml}/\mathrm{s}$, and $0.325\mathrm{ml}/\mathrm{s}$) as a function of delay time, $\tau$.

Since FFT-based autocorrelation calculation provided smoother starting and ending plateaus, it was obviously able to improve curve fitting. For better comparison between the Brownian motion and random flow model, we performed curve fitting by using both models. The Levenberg-Marquardt algorithm was applied in Matlab (Mathwork, Inc., nlinfit function). As indicated by Fig. 7, which was the autocorrelation function of the flow with speed of $0.325\mathrm{ml}/\mathrm{s}$, the Brownian motion model fit the experimental data better than the random flow model did. This result also suggested that our self-designed phantom is able to mimic the capillary flow of the human body. Therefore, we conducted curve fitting by using the Brownian motion model, and the blood flow index ($\mathrm{BFI}=\alpha D_B$) was obtained. The fitting result in general depends highly on the methods of data sampling. In order to obtain better fitting with a uniform weightage in the logarithmic time scale, we sampled the entire autocorrelation function equally on 2-base logarithm scale as shown in Fig. 7. This subsampling helps with faster fitting speed without compromising its accuracy.

Fig. 7

Comparison between Brownian motion and the random flow model. The autocorrelation function was calculated from the flow phantom data with a speed of $0.325\mathrm{ml}/\mathrm{s}$. The entire autocorrelation data was sampled and shown as dots.

In order to compare the goodness of fitting of the FFT-based software autocorrelator with that of the simulated hardware autocorrelator, we calculated a 95% confidence interval of fitted parameters by using an nlparci function in Matlab (Mathwork, Inc.,). As shown in Fig. 8, by using a speed of $0\mathrm{ml}/\mathrm{s}$ as a base, relative $\alpha D_B$ and $\sqrt{\alpha \langle V^2\rangle }$ with their standard deviations at different flow speeds were calculated in both software and hardware cases. It can be found that in either the Brownian motion or random flow model, the standard deviations of the fitted parameter in FFT-based software autocorrelator were always smaller than in the hardware one. Concurrently, in the same autocorrelator, the parameter fitted for the Brownian motion had a smaller standard deviation than that for the random flow model and it also can be observed from individual fitting as shown in Fig. 7. Moreover, we performed a linear regression of the relative $\alpha D_B$ and $\sqrt{\alpha \langle V^2\rangle }$. It showed higher linearity to flow speed in the software autocorrelator case than in the hardware one. In addition, as mentioned previously, the autocorrelation calculating speed is faster than the data writing speed, our instrument holds the potential to make full use of this time to accomplish data fitting and provide the rBF values in real-time.

Fig. 8

Linearity test of relative BFI (either $\alpha D_B$ or $\sqrt{\alpha \langle V^2\rangle }$) versus flow speed. (a) and (c) are obtained using a software autocorrelator (shown in blue), while (b) and (d) are obtained using a hardware autocorrelator (shown in red). All rBFIs are normalized to the BFI value at zero flow.

5.2.

In Vivo Human Arm Cuff Occlusion Experiment

Figure 9(a) shows the normalized intensity autocorrelation function $g_2(\tau )$ of three stages of the in vivo cuff occlusion experiment, namely baseline, cuff inflation, and cuff deflation. The baseline curve was averaged over the first 15 s, the inflation curve was averaged over the 15th to 35th seconds, while the deflation curve was the average of only the 35th to 45th seconds. The decay rate of $g_2(\tau )$ shows significant change. Figure 9(b) reveals rBF, which was obtained by fitting $g_2(\tau )$ to the solution of the photon diffusion correlation equation. For the first 15 s, the blood flow was constant as no constriction was applied. Once the constriction was applied to the cuff, the blood flow dropped down sharply and when constriction was released, the blood flow immediately showed an overshoot followed by a slow decrease to the baseline due to homeostasis.

Fig. 9

(a) Autocorrelation functions of cuff occlusion stages, namely baseline (blue solid line), cuff inflation (red dash-dot line), and cuff deflation (green dashed line). (b) Representative time response of relative blood flow (rBF) during cuff occlusion experiment, including baseline recording, cuff inflation (cuff occlusion period on the figure), and cuff deflation stages.

5.3.

PDT Monitoring Experiment

We plotted the DCS monitoring result during the PDT treatment in Fig. 10. The color bars in the figure represent PDT light intervals.

Fig. 10

The time behavior of relative blood flow (rBF) before, during, and after PDT treatment. The color bars represent the light exposures during treatment; there was no data recorded during light exposure intervals.

The mean tumor rBF decreased by $\sim 50\%$ immediately after the whole PDT, and decreased to $\sim 45\%$ of the baseline 60 min after treatment. During each PDT dark interval, blood flow changes showed minor fluctuations. They indicate impermanent vessel occlusion during PDT light intervals. These temporary vessel occlusions can be reversible, and it allows flow recovery during the dark intervals. A similar trend was observed using laser Doppler measurements, although different animal models were used in each case.⁴⁶ However, the results after 60 min suggest that permanent vessel occlusion had occurred due to PDT.