5320
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
Assessing the reliability of diffuse correlation spectroscopy models on noise-free analytical Monte Carlo data TIZIANO BINZONI1,2,*
AND
FABRIZIO MARTELLI3
1
Département de Neurosciences Fondamentales, University of Geneva, Switzerland Département de l’Imagerie et des Sciences de l’Information Médicale, University Hospital, Geneva, Switzerland 3 Dipartimento di Fisica e Astronomia dell’Università degli Studi di Firenze, Sesto Fiorentino, Firenze, Italy *Corresponding author: [email protected] 2
Received 13 March 2015; revised 30 April 2015; accepted 18 May 2015; posted 19 May 2015 (Doc. ID 236202); published 5 June 2015
It is shown that an analytical noise-free implementation of Monte Carlo simulations [Appl. Opt. 54, 2400 (2015)] for diffuse correlation spectroscopy (DCS) may be successfully used to check the ability of a given DCS model to generate a reliable estimator of tissue blood flow. As an example, four different DCS models often found in the scientific literature are tested on a simulated tissue (semi-infinite geometry) with a Maxwell– Boltzmann probability distribution function for red blood cell speed. It is shown that the random model is the best model for the chosen speed distribution but that (1) some inaccuracies in the DCS model in taking into account red blood cell concentration and (2) some inaccuracies, probably due to a low-order approximation of the DCS model, are still observed. The method can be easily generalized for other speed/flow probability distribution functions of the red blood cells. © 2015 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (170.3660) Light propagation in tissues; (170.3340) Laser Doppler velocimetry; (290.5825) Scattering theory. http://dx.doi.org/10.1364/AO.54.005320
1. INTRODUCTION Diffuse correlation spectroscopy (DCS) is a well-known optical technique allowing the assessment of mean tissue blood flow (TBF), noninvasively, in humans (for a review, see e.g. [1,2]). Besides the performance of the hardware, the reliability of DCS depends on the right choice made for the algorithms utilized to derive TBF. In general, TBF is obtained by fitting a suitable analytical model to the data generated by the DCS hardware. The fitting parameters are expected to be related in some way to the TBF. In human tissues, the vascular tree has a complex and partially unknown geometrical structure. The speed distribution of blood covers a large range of values, even inside tissue volumes of a few millimeters cubed. To develop reliable DCS algorithms, this physiological/morphological information should be taken into account as best as possible [3–5]. Moreover, DCS models are derived from an approximation of the correlations transport equation [6], introducing further inaccuracies in the final model. For these reasons, and also for the difficulty in finding reliable tissue flow phantoms to test the DCS instrumentation [7], it remains quite difficult in practice to judge the quality of a chosen DCS analytical model. 1559-128X/15/175320-07$15/0$15.00 © 2015 Optical Society of America
Given the high complexity of the whole problem from the point of view of modeling, it may be useful to study the different DCS models in a controlled situation where all parameters are a priori perfectly known, and in these conditions operate the retrieval of the TBF with the models. With this approach, measurement errors deriving from hardware specifications must not be taken into account, because the aim is to investigate more fundamental properties of the DCS modeling. Although this fact may give the impression of limitation in the present investigation, it actually allows one to have intrinsic performance tests of the different DCS models from which it may be possible to retrieve unique basic information on the modeling of TBF. This study may also facilitate the comprehension of the reasons why, when DCS is applied to real tissues, some DCS models may be preferable to others in terms of agreement with the measurements. Thus, the aim of the present contribution is to make a step further by proposing a new tool facilitating the testing of DCS models. This will be realized by generating sets of synthetic DCS data, where the underlying blood speed/flow probability distribution function is perfectly known. The data will be obtained by exploiting a previously developed noise-free analytical Monte Carlo (MCan ) [8]. As a demonstrative example, the
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
resulting sets of data will be used to test the reliability of four different DCS models classically found in the scientific literature dealing with TBF in humans.
2. METHODS The fact that a mathematical model perfectly fits the experimental data does not imply that the estimated parameters reliably reproduce the corresponding real physical/physiological quantities of interest. For this reason, it is fundamental to test the model on a set of synthetic/experimental data where the actual values of physical/physiological quantities are known. This allows one to check if the model is able to estimate the unknown quantities (e.g., TBF in our case). In Section 2.A we will show how it is possible to generate DCS synthetic data, where the velocity and flow probability density distribution of the moving scatterers (e.g., red blood cells) is perfectly known. In Section 2.B four analytical models, usually used to derive TBF by DCS on humans, are defined. In Section 2.C the optical and physiological parameters for the different simulations are given, and the procedure utilized to fit the four models to the synthetic data is explained.
pv ‖⃗v ‖
The first term on the right-hand side of the equality allows one to obtain the expected value of 1 for g 2MC ∞. Different realizations of the function g 2MC τ, for different physiological/ optical parameter values (see Section 2.C) and valid for a semi-infinite medium, have been generated to test the four analytical DCS models described in Section 2.B. The percentage of photons that have interacted at least with one moving scatterer before reaching the detector, P int , has also been estimated by MCan as P int 1 − p¯ nms 0 × 100;
(2)
where p¯ nms 0 {see Eq. (26) in [8]} is the probability for a photon to undergo no scattering events with a moving scatterer. This information will be utilized for the discussion in Section 4. The probability density function, pv ‖⃗v ‖, for the scatterer speed, ‖⃗v ‖, simulating percolating blood through a randomly directed network of capillaries has been chosen as already proposed elsewhere [10–12]:
1 2
3 − 2 ‖⃗v ‖−‖⃗vt ‖2 − 32 ‖⃗v‖‖⃗vt ‖2 ‖⃗v ‖ 2v 2v B B × e −e ; v B ‖⃗v t ‖ (3)
where v⃗ is a random variable representing the scatters velocity, v⃗t a constant velocity vector representing a bulk translation of the scatterers, v B the parameter defining a Gaussian distribution for the v⃗ coordinates, and as a consequence, generates the typical Maxwell–Boltzmann distribution (when v⃗t 0) for ‖⃗v ‖. In vivo, vB is intuitively linked to the disordered/ random movement of the red blood cells inside the extremely complex and contorted capillary network (with no net movement). The larger the v B value, the larger the shape of the ‖⃗v ‖ distribution, pv ‖⃗v ‖. The velocity v⃗t is linked to the global movement (net movement) of the blood that by definition must enter and then leave a given tissue volume (red blood cells must return to the lung to be reoxygenated before passing again through the same tissue volume). The mean speed of the scatterers, h‖⃗v ‖i, can be obtained as Z ∞ ‖⃗v ‖pv ‖⃗v ‖d‖⃗v ‖ h‖⃗v ‖i 0
A. Analytical Monte Carlo: Reference Data
The main signal generated by the DCS hardware is the autocorrelation function of the detected photoelectric current [it, where t is the time]. In the present context, the abovementioned (normalized) intensity autocorrelation function [g 2MC τ, where τ is the lag time] has been simulated by exploiting the results produced by the MCan approach. For reasons of limited space, the reader is invited to refer to the original literature for complete information on the MCan method [8]. The power density spectrum [P MC Δω, where Δω is the angular frequency] of it has been obtained as previously shown by MCan [8], and g 2MC τ has been derived from the well-known relationship [9] Z 1 ∞ P MC Δωe iΔωτ dΔω: (1) g 2MC τ 1 2π −∞
3 2π
5321
v2B
3‖⃗v t ‖2 erf 3‖⃗v t ‖
pffiffiffipffiffiffi pffiffiffipffiffiffi 2 3 2 ‖⃗v t ‖ 3 2v B −32‖⃗vvt2‖ e B ; 3π vB 2 (4)
and the second moment as Z ∞ 2 h‖⃗v ‖ i ‖⃗v ‖2 pv ‖⃗v ‖d‖⃗v ‖ v 2B ‖⃗v t ‖2 :
(5)
0
By knowing the probability for a photon to scatter with a moving scatterer, P m (see Section 2.C), it is possible to assess the reference scatterers flow, Φ (e.g., blood flow), as Φ ∝ P m h‖⃗v‖i:
(6)
The parameter Φ is the quantity that the experimentalist or the medical doctor want to obtain by DCS. In vivo, the parameter P m is linked to the number of moving red blood cells, and this explains why P m h‖⃗v ‖i in Eq. (6) is related to a blood flow. In the following sections, we will investigate if the different analytical models (see Section 2.B) are actually able to estimate Φ when a real vascular network has morphological/physiological properties represented, e.g., by Eq. (3). To give information on the g 2MC τ relaxation time, the constant τ1.05 has been defined as (see Fig 3) g 2 τ1.05 1.05:
(7)
[When τ increases from 0 to ∞, g 2MC τ always decreases from 2 to 1.] B. Four Analytical DCS Models
The four representative DCS models for g 2 τ usually utilized to fit the experimental DCS data [i.e., g 2MC τ] are derived from the field autocorrelation function, g 1 τ, as [13] g 2 τ 1 jg 1 τj2 ;
(8)
where g 1 τ
G 1 τ ; G 1 0
(9)
5322
and G 1 τ and where
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
3μs0 4π
e
−K r 1
r1
−
e
−K r 2
r2
;
2 2 0 02 2πn hΔr 2 τi: K 3μa μs P m μs λ
(10)
(11)
The parameters λ, n, μa , and μs0 are the laser wavelength, the refractive index, and the absorption and reduced scattering coefficients of the medium, respectively. The parameter ρ is the interoptode distance, r 1 ρ2 z 20 1∕2 , r 2 ρ2 z 0 2z e 2 1∕2 , z 0 1∕μs0 , and z e A∕μs0 , where A is computed using Ref. [14]. The four models, for a the semi-infinite medium, have been derived by suitably defining the mean squared displacement of the moving scatterers, hΔr 2 τi (see Sections 2.B.1–2.B.4). 1. Brownian Model
In this model, the moving scatterers are considered to follow a Brownian movement and for this reason hΔr 2 τi is expressed as [2] hΔr 2 τi 6DB τ:
(12)
In this case the moving scatterers flow, ΦB , is considered to be proportional to [2] ΦB ∝ P m DB ;
(13)
where DB is the diffusion coefficient related to the scatterers movement. In vivo, DB describes the rate at which red blood cells move through a particular tissue with given physiological/ physical conditions. Thus, considering that blood evolves inside the vessels, if its speed increases (blood flow increases) DB must also increase. This gives an intuitive explanation of the relationship existing between ΦB and DB . 2. Random Model
Another well-known and general model is when the moving scatterers follow a random movement. In this case, hΔr 2 τi is expressed as [2] hΔr 2 τi hV 2R iτ2 ;
(14)
and the flow, ΦR , as ΦR ∝ P m hV 2R i; where
hV 2R i
(15)
is the second moment of the scatterers speed.
where DL is the effective diffusion coefficient and τL is the time scale of the randomization of the scatterers velocity. 4. Mixed Model
Another model trying to better reproduce the physical/ physiological reality and allowing one to simultaneously take into account Brownian and random flow has been proposed in Ref. [17] and is expressed as hΔr 2 τi 6DM τ hV 2M iτ2 : In this case the flow, ΦM , is estimated by ΦM ∝ P m D M ; where DM is the diffusion coefficient and moment of the scatterers speed.
(19) hV 2M i
In Ref. [15] it is explained that the Brownian model (Section 2.B.1) assumes that the ballistic to random walk transition in the diffusion of the moving scatterers occurs at time scales shorter than that probed by DCS. Considering that there are no data to support this assumption, a more general formulation releasing the assumption has been proposed for hΔr 2 τi, i.e., [15,16] −ττ 2 ; (16) hΔr τi 6DL τ − τL 1 − e L and the flow, ΦL , is expressed as ΦL ∝ P m D L ;
(17)
is the second
5. “Exact” Model
The previous four models are typically utilized in real experimental conditions, where it is not known if the underlining microvascular geometry and blood flow speed distribution actually correspond to the chosen model. In the present case, the exact speed distribution is known [Eq. (3)]. Thus, one can derive an “exact” model corresponding to this specific distribution. The aim here is to further test the common mathematical techniques utilized to derive all the previous four models. By exactly following the same mathematical procedure utilized to derive the previous four models, it is straightforward to demonstrate [18] that if a semi-infinite medium has a speed distribution for the scatterers equal to Eq. (3), Eq. (10) always holds and in this case hΔr 2 τi v2B ‖⃗v t ‖2 τ2 h‖⃗v ‖2 iτ2 :
(20)
Equation (20) is similar to Eq. (14), and for this reason we expect hV 2R i and h‖⃗v ‖2 i to be equal if fitted to the present MCan data. It is essential to note again at this point that models typically utilized in DCS are P 1 approximations of the correlation transport equation [6,18]. This means that some parameters appearing in the final equations must satisfy some mathematical constraints. This is the case for τ, which cannot be larger than the critical constant τc (if not applicable, the model brakes). For the “exact” model, τc is equal to −1∕2 1 ; (21) τc k2o v 2B ‖⃗v t ‖2 3 where k o 2πn∕λ; and the constraint may be expressed as τ1.05 < τc :
3. Langevin Model
(18)
(22)
Intuitively, this means that the g 2 τ curve (“exact” model) must roughly reach the value 1 in a time τc (see Fig. 3), determined by the parameters v B and ‖⃗v t ‖. If this is not the case, the model brakes. Note also that Eqs. (14) and (20) are formally equivalents. For this reason, in the concluding Section 4 we will use the terms random or “exact” indistinctly. C. Simulations and Fitting Procedure
The reference data [i.e., g 2MC τ curves] have been generated for a semi-infinite medium with refractive index n 1.4. This n value implies A ≈ 1.96; precisely calculated “on flight” during the simulations (see Section 2.B). The optical fibers (normal to
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
5323
the medium surface) spacing was set to ρ 30 mm and the CW laser wavelength to λ 785 nm. Two sets of optical parameters were utilized for the simulations, denoted data set 1 and data set 2, generating the reference g 2MC τ curves. To facilitate the reading, consider that the set of parameters containing μa 0.005 mm−1 and μs0 1 mm−1 (see the following subsections) are reasonably within the range of the typical values for a bone tissue; they have been called “bone” [19]. For the same reason, the set of parameters containing μa 0.025 mm−1 and μs0 0.5 mm−1 have been called “muscle.” 1. Data Set 1: “Bone”
All combinations of the following parameters were considered: μa 0.005 mm−1 , μs0 1 mm−1 , P m ∈ f0.01; 0.025; 0.05; 0.075; 0.1g [8,20], v B ∈ f0.025; 0.05; 0.075; 1; 1.25; 1.5; 1.75; 2g mm s−1 , ‖⃗v t ‖ ∈ f0; 0.025; 0.05; 0.075; 1; 1.25; 1.5; 1.75; 2g mm s−1 , g 0.9, and λ 785 nm. This gave a total of 360 reference g 2MC τ curves. 2. Data Set 2: “Muscle”
The same as in Section 3.C.1 with the exception that μa 0.025 mm−1 and μs0 0.5 mm−1 . This gave a total of 360 reference g 2MC τ curves. 3. Fitting Procedure and Flow Assessment
The reference g 2MC τ curves were fitted using the four investigated DCS models, represented by g 2 τ and the relative Eqs. (12), (14), (16), and (18). A trust-region-reflective algorithm (MATLAB language function), minimizing the squared l 2 -norm of the residuals, e 2 min ‖g 2MC τ − g 2 τ‖22 ; C
(23)
was applied to perform the fitting; where C refers as usual to a particular set of fitting parameters, depending on the chosen DCS model. For the Brownian model C ∈ fP m DB g (to be seen as a single fitting parameter), for the random model C ∈ fP m hV 2R ig (one fitting parameter), for the Langevin model C ∈ fP m DL ; τL g (two fitting parameters), and for the mixed model C ∈ fP m DM ; P m hV 2M ig (two fitting parameters). It is important to note that when analyzing a real experiment, P m can never be isolated as a single fitting parameter. The above fitting parameters gave the ΦB , ΦR , ΦL , and ΦM estimators of Φ (in arbitrary units) for different physiological conditions and tissues. It is clear that the resulting estimators, obtained using Eqs. (12), (14), (16), and (18), do not strictly correspond to a flow, but this is actually the way the flow is estimated in DCS. This is a further reason explaining why it is so important to test the models on situations where the flow values are known. For explanatory purposes, the “exact” model was also fit to the reference g 2MC τ curves, but in this case P m was considered as known (this is not possible in a real experiment), and h‖⃗v ‖2 i was taken as a fitting parameter. By knowing the v B and ‖⃗v t ‖ values of the simulations, this approach allowed us to investigate if the “exact” model actually reproduces the relationship v 2B ‖⃗v t ‖2 h‖⃗v ‖2 i [see Eq. (20)], and thus demonstrate its reliability.
Fig. 1. Number of simulations, nsim , that have a percentage P int of photons that have interacted at least once with a moving scatterer before reaching the photodetector.
3. RESULTS In the present contribution, we have simulated 720 reference g 2MC τ curves, i.e., 360 simulations for “bone” and 360 simulations for “muscle.” From Fig. 1 one can see that all the 360 simulations for “bone” tissue have P int ≈ 100% of the photons reaching the photodetector that have interacted at least once with a moving scatterer. This, independently of the chosen simulation parameters. However, 144 simulations for “muscle” tissue have P int values slightly below 100%. This means that some photons have not yet interacted with a moving scatterer when reaching the photodetector. For the sake of completeness, in Fig. 2 is shown an example comparing P MC Δω and g 2MC τ obtained by MCan with the relative simulations obtained by a classical Monte Carlo (MC) method [21]. As already shown in a previous study, MCan and MC produce the same results [8]. However, it is important to note that the MC computing time for 720 P MC Δω spectra would be of the order of 1340 days while for MCan only 24 h were necessary.
Fig. 2. Example of P MC Δω and g 2MC τ computed both by classical MC and MCan . The simulation parameters were: μa 0.025 mm−1 , μs0 0.5 mm−1 , P m 0.05, v B 1, mm s−1 , ‖⃗v t ‖ 2 mm s−1 , g 0.9, n 1.4, and λ 785 nm.
5324
Vol. 54, No. 17 / June 10 2015 / Applied Optics
Research Article
Fig. 3. Example of g 2MC τ curve (black line) with the relative g 2 τ fittings (remaining colors) for the four different models. The dashed vertical lines represent the position of τ1.05 [Eq. (7)] and τc [Eq. (21)] on the abscissa.
In Fig. 3, a typical g 2MC τ curve with the relative g 2 τ fitting results for the four models (see Section 2.B) are presented. By visual inspection it clearly appears that the Brownian model is the worst one. To give a more quantitative estimate of the quality of the fittings for the four investigated models, in Fig. 4 are reported the mean and standard deviation (720 curves) of the squared 2-norm of the residuals, e 2 [Eq. (23)], for each model. As expected (see Section 4), the Brownian model is the worst model and, in general (not shown), it reproduces fitted curves similar to the one shown in Fig. 3 (blue dotted line). In Fig. 5 are shown the different flow estimators ΦB , ΦR , ΦL , and ΦM , compared to the exact value Φ [Eq. (6)]. It immediately appears that none of the models correctly takes into account the dependence on P m . A methodical P m -induced shift in the curves is observed. It is worth remarking that the Brownian model, ΦB , curiously generates reasonably “linear” data for a given tissue. However, considering that the Brownian fitted curves appear in general as in Fig. 3, i.e., the model does not fit at all the curves, this result must be considered as a pure mathematical artifact. Moreover, the
Fig. 4. Squared 2-norm of the residuals, e 2 [Eq. (23)], for the four fitted models (360 simulations for “bone” and 360 simulations for “muscle”). Vertical bars are standard deviations.
Fig. 5. Estimation of the scatterers’ flow, Φ, assessed by means of the four investigated models. P m 0.01 (black), P m 0.025 (red), P m 0.05 (blue), P m 0.075 (green), and P m 0.1 (magenta).
Brownian model is not able to correct for μa and μs0 differences [compare ordinate values, for the same Φ, in Figs. 5(a) and 5(b)]. Clearly, ΦL and ΦM are also not good estimators, even if they have the smallest e 2 , because the data are scattered on too large a region compared with the others models (good fitting, bad Φ estimator). In conclusion, the best estimator is represented by the random model, even if ΦR is not perfectly linearly related to Φ, and that the influence of P m is not correctly taken into account (shift in the colored curves). This is not astonishing because the random model Eq. (14) is in practice formally equal to the “exact” model Eq. (20). Let us now investigate if the “exact” model itself has not some unwanted hidden defect. While in a real experimental setup hV 2R i cannot be computed independently of P m , in the present study this has been done because P m was obviously known because it was utilized to generate the reference MCan curves. In this manner it has been possible to compare h‖⃗v ‖2 i obtained with the simulation parameters for MCan [Eq. (5)] and h‖⃗v ‖2 i estimated by fitting; theoretically, they should give exactly the same values (see Section 2.B.5). In Fig. 6 one can see that hV 2R i values for “bone” follows the identity line quite well, and that the error depending on P m is not visible. Thus, the “exact” model seems to be able to reasonably generate for this set of data the exact estimators. However, the results for “muscle,” even if linear, are less precise and the P m -induced shifts are always present for small P m values. It is worth noting that in this case the deviation from the identity line observed in Fig. 6(b) can be in part explained by plotting only the points
Research Article
Fig. 6. Panels (a) and (b) comparison of hV 2R i, assessed by means of the random model, with the exact value h‖⃗v ‖2 i Eq. (5), and supposing P m is known. This case does not correspond to the experimental reality where P m and hV 2R i cannot be separated (see manuscript). P m 0.01 (black), P m 0.025 (red), P m 0.05 (blue), P m 0.075 (green), and P m 0.1 (magenta). The dashed line is the identity line. Panels (c) and (d) similar to (a) and (b), but where only the points satisfying the constraint τ1.05 < τc (see Section 4) are represented. The data point represented by the red star is computed from the data in Fig. 3 (random model).
that are obtained from g 2MC τ curves satisfying the inequality τ1.05 < τc (see Section 2.B.5); a necessary condition allowing one to derive the model represented by Eq. (10). Figures 6(c) and 6(d) clearly show that if τ1.05 < τc , then the fitted hV 2R i are better distributed on the identity line. The red star appearing in Fig. 6(d) represents a particular hV 2R i value computed from the example (random model) of Fig. 3, where the constraint τ1.05 < τc is also valid. 4. DISCUSSION AND CONCLUSIONS In the present contribution we have shown that the MCan approach allows one to investigate if a classical DCS model is able to reasonably estimate scatters flow, ϕ (e.g., blood flow), in arbitrary units. It has been shown that the random model, among the proposed four models, is the best model allowing one to describe a flow that has a Maxwell–Boltzmann distribution for the moving scatterers speed. Even if the random or “exact” models have been built ad hoc for the Maxwell–Boltzmann speed distribution, this result remains nontrivial. In fact, we have seen (Section 2.B.5) that the derivation of g 2 τ implies a well-known specific approximation (oτ3 , [18]), on the field-autocorrelation function for a single scattering event, that obviously do not appear in the MCan -generated g 2MC τ curves. This means that g 2 τ and g 2MC τ are not necessarily the same, even for the “exact” model. This problem has been highlighted in Figs. 6(a) and 6(b) where, even after having taken into account the τ1.05 < τc constraint, some differences with the identity line still remain. The remaining inaccuracies cannot ensue from the reference data noise, because MCan is noise free. The inaccuracies are probably due, as highlighted above, to the too low order approximation utilized when deriving Eq. (20). In practice, this means that if the v 2B or ‖⃗v t ‖ (or a combination) values are too
Vol. 54, No. 17 / June 10 2015 / Applied Optics
5325
large, the random (or “exact”) model may be not valid. This point will certainly be a matter for future investigations. In the present work, it has also been shown [Figs. 5(c) and 5(d)] that when estimating ΦR , an inexact dependence on P m is observed, and that the ΦR − Φ relationship is not perfectly linear. One of the reasons might be the fact that the estimator ΦR contains the second moment of the speed and not the speed itself [see Eq. (15)]. However, the positive point is that for given real tissue (in vivo), the parameter ΦR (proportional to the number of moving scatterers in the observed tissue) is not subject to such large P m changes. Thus, ΦR probably moves in a neighborhood along one given colored line in Figs. 5(c) and 5(d). This possibly explains also why DCS appears to give reasonable results in real experimental setups, even in the presence of this mathematical problem. In the past it has been noted that, experimentally, the photons detected by the DCS instrumentation, scatter not only on moving scatterers, but also on a fixed background [18]. This fact is known to introduce a distortion in the g 2 τ function (a phenomenon not described by Monte Carlo simulations). It has been hypothesized that the good correlation of the DCS models with the experimental measurements, even in the presence of fixed background scattering, may be due to the fact that most of the photons interact at least once with a moving body [2]. By inspecting Fig. 1 this seems to be true because of the high P int values. However, this remain a necessary but not sufficient condition, and other investigations must be performed on this topic. Summarizing the main concepts, we can say that the proposed method allows one to test DCS analytical models on synthetic reference data, g 2MC τ, precisely reproducing a perfectly known moving scatterers speed/flow distribution. In fact, we have seen that DCS analytical models typically utilized in biomedical optics are only approximated solutions. This means that DCS analytical models may not exactly reproduce a fortiori the synthetic data, g 2MC τ, or generate the correct desired parameter values. Thus, the proposed method allows one to study, e.g., the domain of validity of the DCS models (e.g., speed range where they remain valid) or their robustness on synthetic data where the speed distribution is not exactly the one considered by the DCS model. Conditions that may easily be present in the real experimental setting. Thus, the proposed approach is to be considered before any other test such as the study of the influence of noise or hardware artifacts [3,22,23]. If the proposed test on g 2MC τ reference data fails, e.g., for a particular range of speed parameters, then it is theoretically pointless to apply the model in this situation. This, because there is an obvious intrinsic problem in the DCS model, is independent from noise or hardware artifacts. In this case, the possibility of a further improvement of the model must be considered. It must be noted that in vivo the Brownian model better fits the experimental data than the random model [2,24–27]. For demonstrative purposes in the present contribution, the model utilized for the moving scatterers speed was a Maxwell– Boltzmann distribution (with a bulk drift). Thus, it is mathematically expected that the best DCS model (without considering the observed approximation errors) is in this case
5326
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
the (“exact”) random model (see Section 2.B.5). Thus, there is no inconsistency between the present results and the experimental one. These results only show that the classical Brownian model is not able to extract the right physiological parameters if the speed distribution is Maxwell–Boltzmann, and that the Brownian model is not robust enough in this case. Considering that the probability density function of the speed in real tissue is not exactly known, future tests of the robustness of the Brownian model for other speed distributions will certainly be of interest. The Langevin model deserves a last comment. In fact, it represent a clear example where a perfect fitting does not mean that the derived parameters of interest are necessarily exact (see Figs. 4 and 5). This is probably due to the cross talk existing between DL and τL . This motivates again the interest of the tests proposed in the present contribution. In conclusion, we have demonstrated that the MCan -based approach may be used to check the reliability of a given DCS model to generate the suitable Φ estimator, when applied to a given moving scatteres speed distribution. Moreover, it must be noted that MCan can be adapted to other speed distributions, simulating in this way other experimental conditions (e.g., satisfying the Brownian DCS model). This will allow, for example, one to test in the future the robustness of a given DCS model on flow distributions that do not exactly represent the conditions for which the model was built (this is probably what happens in a real experiment). We hope that the present work will contribute to gaining a further understanding about the DCS technique. REFERENCES 1. D. A. Boas, L. E. Campbell, and A. G. Yodh, “Scattering and imaging with diffusing temporal field correlations,” Phys. Rev. Lett. 75, 1855– 1858 (1995). 2. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73, 076701 (2010). 3. Y. Shang, T. Li, L. Chen, Y. Lin, M. Toborek, and G. Yu, “Extraction of diffuse correlation spectroscopy flow index by integration of Nth-order linear model with Monte Carlo simulation,” Appl. Phys. Lett. 104, 193703 (2014). 4. Y. Shang and G. Yu, “A Nth-order linear algorithm for extracting diffuse correlation spectroscopy blood flow indices in heterogeneous tissues,” Appl. Phys. Lett. 105, 133702 (2014). 5. W. B. Baker, A. B. Parthasarathy, D. R. Busch, R. C. Mesquita, J. H. Greenberg, and A. G. Yodh, “Modified Beer–Lambert law for blood flow,” Biomed. Opt. Express 5, 4053–4075 (2014). 6. B. J. Ackerson, R. L. Dougherty, N. M. Reguigui, and U. Nobbman, “Correlation transfer: application of radiative transfer solution methods to photon correlation problems,” J. Thermophys. Heat Transfer 6, 577–588 (1992). 7. T. Binzoni, A. Torricelli, R. Giust, B. Sanguinetti, P. Bernhard, and L. Spinelli, “Bone tissue phantoms for optical flowmeters at large interoptode spacing generated by 3D-stereolithography,” Biomed. Opt. Express 5, 2715–2725 (2014). 8. T. Binzoni and F. Martelli, “Random numbers free analytical implementation of Monte Carlo for laser-Doppler flowmetry at large
9.
10. 11.
12.
13.
14.
15.
16. 17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
interoptode spacing: application to human bone tissue,” Appl. Opt. 54, 2400–2406 (2015). H. Z. Cummins and H. L. Swinney, “Light beating spectroscopy,” in Progress in Optics, E. Wolf, ed. (North-Holland, 1970), Vol. 8, pp. 133–200. R. Bonner and R. Nossal, “Model for laser Doppler measurements of blood flow in tissue,” Appl. Opt. 20, 2097–2107 (1981). J. Zhong, A. M. Seifalian, G. E. Salerud, and G. E. Nilsson, “A mathematical analysis on the biological zero problem in laser Doppler flowmetry,” IEEE Trans. Biomed. Eng. 45, 354–364 (1998). T. Binzoni and D. Van De Ville, “Full-field laser-Doppler imaging and its physiological significance for tissue blood perfusion,” Phys. Med. Biol. 53, 6673–6694 (2008). C. Cheung, J. P. Culver, K. Takahashi, J. H. Greenberg, and A. G. Yodh, “In vivo cerebrovascular measurement combining diffuse near-infrared absorption and correlation spectroscopies,” Phys. Med. Biol. 46, 2053–2065 (2001). D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). S. A. Carp, N. Roche-Labarbe, M. A. Franceschini, V. J. Srinivasan, S. Sakadžić, and D. A. Boas, “Due to intravascular multiple sequential scattering, diffuse correlation spectroscopy of tissue primarily measures relative red blood cell motion within vessels,” Biomed. Opt. Express 2, 2047–2054 (2011). B. J. Berne and R. Pecora, Dynamic Light Scattering: with Applications to Chemistry, Biology, and Physics (Dover, 2000). H. M. Varma, C. P. Valdes, A. K. Kristoffersen, J. P. Culver, and T. Durduran, “Speckle contrast optical tomography: a new method for deep tissue three-dimensional tomography of blood flow,” Biomed. Opt. Express 5, 1275–1289 (2014). D. A. Boas, “Diffuse photon probes of structural and dynamical properties of turbid media: theory and biomedical applications,” Ph.D. thesis (University of Pennsylvania, 1996). P. Farzam, P. Zirak, T. Binzoni, and T. Durduran, “Pulsatile and steady-state hemodynamics of the human patella bone by diffuse optical spectroscopy,” Physiol. Meas. 34, 839–857 (2013). A. Kienle, “Non-invasive determination of muscle blood flow in the extremities from laser Doppler spectra,” Phys. Med. Biol. 46, 1231– 1244 (2001). A. Sassaroli and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A 29, 2110–2117 (2012). L. Dong, L. He, Y. Lin, Y. Shang, and G. Yu, “Simultaneously extracting multiple parameters via fitting one single autocorrelation function curve in diffuse correlation spectroscopy,” IEEE Trans. Biomed. Eng. 60, 361–368 (2013). C. Zhou, G. Yu, D. Furuya, J. Greenberg, A. Yodh, and T. Durduran, “Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain,” Opt. Express 14, 1125– 1144 (2006). R. C. Mesquita, T. Durduran, G. Yu, E. M. Buckley, M. N. Kim, C. Zhou, R. Choe, U. Sunar, and A. G. Yodh, “Direct measurement of tissue blood flow and metabolism with diffuse optics,” Philos. Trans. R. Soc. A 369, 4390–4406 (2011). G. Yu, T. Durduran, C. Zhou, R. Cheng, and A. G. Yodh, “Nearinfrared diffuse correlation spectroscopy for assessment of tissue blood flow,” in Handbook of Biomedical Optics, D. A. Boas, C. Pitris, and N. Ramanujam, eds. (CRC Press, 2011). J. Dong, R. Bi, J. H. Ho, P. S. Thong, K. C. Soo, and K. Lee, “Diffuse correlation spectroscopy with a fast Fourier transform-based software autocorrelator,” J. Biomed. Opt. 17, 0970041 (2012). T. Durduran and A. G. Yodh, “Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement,” Neuroimage 85, 51–63 (2014).
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
Assessing the reliability of diffuse correlation spectroscopy models on noise-free analytical Monte Carlo data TIZIANO BINZONI1,2,*
AND
FABRIZIO MARTELLI3
1
Département de Neurosciences Fondamentales, University of Geneva, Switzerland Département de l’Imagerie et des Sciences de l’Information Médicale, University Hospital, Geneva, Switzerland 3 Dipartimento di Fisica e Astronomia dell’Università degli Studi di Firenze, Sesto Fiorentino, Firenze, Italy *Corresponding author: [email protected] 2
Received 13 March 2015; revised 30 April 2015; accepted 18 May 2015; posted 19 May 2015 (Doc. ID 236202); published 5 June 2015
It is shown that an analytical noise-free implementation of Monte Carlo simulations [Appl. Opt. 54, 2400 (2015)] for diffuse correlation spectroscopy (DCS) may be successfully used to check the ability of a given DCS model to generate a reliable estimator of tissue blood flow. As an example, four different DCS models often found in the scientific literature are tested on a simulated tissue (semi-infinite geometry) with a Maxwell– Boltzmann probability distribution function for red blood cell speed. It is shown that the random model is the best model for the chosen speed distribution but that (1) some inaccuracies in the DCS model in taking into account red blood cell concentration and (2) some inaccuracies, probably due to a low-order approximation of the DCS model, are still observed. The method can be easily generalized for other speed/flow probability distribution functions of the red blood cells. © 2015 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (170.3660) Light propagation in tissues; (170.3340) Laser Doppler velocimetry; (290.5825) Scattering theory. http://dx.doi.org/10.1364/AO.54.005320
1. INTRODUCTION Diffuse correlation spectroscopy (DCS) is a well-known optical technique allowing the assessment of mean tissue blood flow (TBF), noninvasively, in humans (for a review, see e.g. [1,2]). Besides the performance of the hardware, the reliability of DCS depends on the right choice made for the algorithms utilized to derive TBF. In general, TBF is obtained by fitting a suitable analytical model to the data generated by the DCS hardware. The fitting parameters are expected to be related in some way to the TBF. In human tissues, the vascular tree has a complex and partially unknown geometrical structure. The speed distribution of blood covers a large range of values, even inside tissue volumes of a few millimeters cubed. To develop reliable DCS algorithms, this physiological/morphological information should be taken into account as best as possible [3–5]. Moreover, DCS models are derived from an approximation of the correlations transport equation [6], introducing further inaccuracies in the final model. For these reasons, and also for the difficulty in finding reliable tissue flow phantoms to test the DCS instrumentation [7], it remains quite difficult in practice to judge the quality of a chosen DCS analytical model. 1559-128X/15/175320-07$15/0$15.00 © 2015 Optical Society of America
Given the high complexity of the whole problem from the point of view of modeling, it may be useful to study the different DCS models in a controlled situation where all parameters are a priori perfectly known, and in these conditions operate the retrieval of the TBF with the models. With this approach, measurement errors deriving from hardware specifications must not be taken into account, because the aim is to investigate more fundamental properties of the DCS modeling. Although this fact may give the impression of limitation in the present investigation, it actually allows one to have intrinsic performance tests of the different DCS models from which it may be possible to retrieve unique basic information on the modeling of TBF. This study may also facilitate the comprehension of the reasons why, when DCS is applied to real tissues, some DCS models may be preferable to others in terms of agreement with the measurements. Thus, the aim of the present contribution is to make a step further by proposing a new tool facilitating the testing of DCS models. This will be realized by generating sets of synthetic DCS data, where the underlying blood speed/flow probability distribution function is perfectly known. The data will be obtained by exploiting a previously developed noise-free analytical Monte Carlo (MCan ) [8]. As a demonstrative example, the
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
resulting sets of data will be used to test the reliability of four different DCS models classically found in the scientific literature dealing with TBF in humans.
2. METHODS The fact that a mathematical model perfectly fits the experimental data does not imply that the estimated parameters reliably reproduce the corresponding real physical/physiological quantities of interest. For this reason, it is fundamental to test the model on a set of synthetic/experimental data where the actual values of physical/physiological quantities are known. This allows one to check if the model is able to estimate the unknown quantities (e.g., TBF in our case). In Section 2.A we will show how it is possible to generate DCS synthetic data, where the velocity and flow probability density distribution of the moving scatterers (e.g., red blood cells) is perfectly known. In Section 2.B four analytical models, usually used to derive TBF by DCS on humans, are defined. In Section 2.C the optical and physiological parameters for the different simulations are given, and the procedure utilized to fit the four models to the synthetic data is explained.
pv ‖⃗v ‖
The first term on the right-hand side of the equality allows one to obtain the expected value of 1 for g 2MC ∞. Different realizations of the function g 2MC τ, for different physiological/ optical parameter values (see Section 2.C) and valid for a semi-infinite medium, have been generated to test the four analytical DCS models described in Section 2.B. The percentage of photons that have interacted at least with one moving scatterer before reaching the detector, P int , has also been estimated by MCan as P int 1 − p¯ nms 0 × 100;
(2)
where p¯ nms 0 {see Eq. (26) in [8]} is the probability for a photon to undergo no scattering events with a moving scatterer. This information will be utilized for the discussion in Section 4. The probability density function, pv ‖⃗v ‖, for the scatterer speed, ‖⃗v ‖, simulating percolating blood through a randomly directed network of capillaries has been chosen as already proposed elsewhere [10–12]:
1 2
3 − 2 ‖⃗v ‖−‖⃗vt ‖2 − 32 ‖⃗v‖‖⃗vt ‖2 ‖⃗v ‖ 2v 2v B B × e −e ; v B ‖⃗v t ‖ (3)
where v⃗ is a random variable representing the scatters velocity, v⃗t a constant velocity vector representing a bulk translation of the scatterers, v B the parameter defining a Gaussian distribution for the v⃗ coordinates, and as a consequence, generates the typical Maxwell–Boltzmann distribution (when v⃗t 0) for ‖⃗v ‖. In vivo, vB is intuitively linked to the disordered/ random movement of the red blood cells inside the extremely complex and contorted capillary network (with no net movement). The larger the v B value, the larger the shape of the ‖⃗v ‖ distribution, pv ‖⃗v ‖. The velocity v⃗t is linked to the global movement (net movement) of the blood that by definition must enter and then leave a given tissue volume (red blood cells must return to the lung to be reoxygenated before passing again through the same tissue volume). The mean speed of the scatterers, h‖⃗v ‖i, can be obtained as Z ∞ ‖⃗v ‖pv ‖⃗v ‖d‖⃗v ‖ h‖⃗v ‖i 0
A. Analytical Monte Carlo: Reference Data
The main signal generated by the DCS hardware is the autocorrelation function of the detected photoelectric current [it, where t is the time]. In the present context, the abovementioned (normalized) intensity autocorrelation function [g 2MC τ, where τ is the lag time] has been simulated by exploiting the results produced by the MCan approach. For reasons of limited space, the reader is invited to refer to the original literature for complete information on the MCan method [8]. The power density spectrum [P MC Δω, where Δω is the angular frequency] of it has been obtained as previously shown by MCan [8], and g 2MC τ has been derived from the well-known relationship [9] Z 1 ∞ P MC Δωe iΔωτ dΔω: (1) g 2MC τ 1 2π −∞
3 2π
5321
v2B
3‖⃗v t ‖2 erf 3‖⃗v t ‖
pffiffiffipffiffiffi pffiffiffipffiffiffi 2 3 2 ‖⃗v t ‖ 3 2v B −32‖⃗vvt2‖ e B ; 3π vB 2 (4)
and the second moment as Z ∞ 2 h‖⃗v ‖ i ‖⃗v ‖2 pv ‖⃗v ‖d‖⃗v ‖ v 2B ‖⃗v t ‖2 :
(5)
0
By knowing the probability for a photon to scatter with a moving scatterer, P m (see Section 2.C), it is possible to assess the reference scatterers flow, Φ (e.g., blood flow), as Φ ∝ P m h‖⃗v‖i:
(6)
The parameter Φ is the quantity that the experimentalist or the medical doctor want to obtain by DCS. In vivo, the parameter P m is linked to the number of moving red blood cells, and this explains why P m h‖⃗v ‖i in Eq. (6) is related to a blood flow. In the following sections, we will investigate if the different analytical models (see Section 2.B) are actually able to estimate Φ when a real vascular network has morphological/physiological properties represented, e.g., by Eq. (3). To give information on the g 2MC τ relaxation time, the constant τ1.05 has been defined as (see Fig 3) g 2 τ1.05 1.05:
(7)
[When τ increases from 0 to ∞, g 2MC τ always decreases from 2 to 1.] B. Four Analytical DCS Models
The four representative DCS models for g 2 τ usually utilized to fit the experimental DCS data [i.e., g 2MC τ] are derived from the field autocorrelation function, g 1 τ, as [13] g 2 τ 1 jg 1 τj2 ;
(8)
where g 1 τ
G 1 τ ; G 1 0
(9)
5322
and G 1 τ and where
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
3μs0 4π
e
−K r 1
r1
−
e
−K r 2
r2
;
2 2 0 02 2πn hΔr 2 τi: K 3μa μs P m μs λ
(10)
(11)
The parameters λ, n, μa , and μs0 are the laser wavelength, the refractive index, and the absorption and reduced scattering coefficients of the medium, respectively. The parameter ρ is the interoptode distance, r 1 ρ2 z 20 1∕2 , r 2 ρ2 z 0 2z e 2 1∕2 , z 0 1∕μs0 , and z e A∕μs0 , where A is computed using Ref. [14]. The four models, for a the semi-infinite medium, have been derived by suitably defining the mean squared displacement of the moving scatterers, hΔr 2 τi (see Sections 2.B.1–2.B.4). 1. Brownian Model
In this model, the moving scatterers are considered to follow a Brownian movement and for this reason hΔr 2 τi is expressed as [2] hΔr 2 τi 6DB τ:
(12)
In this case the moving scatterers flow, ΦB , is considered to be proportional to [2] ΦB ∝ P m DB ;
(13)
where DB is the diffusion coefficient related to the scatterers movement. In vivo, DB describes the rate at which red blood cells move through a particular tissue with given physiological/ physical conditions. Thus, considering that blood evolves inside the vessels, if its speed increases (blood flow increases) DB must also increase. This gives an intuitive explanation of the relationship existing between ΦB and DB . 2. Random Model
Another well-known and general model is when the moving scatterers follow a random movement. In this case, hΔr 2 τi is expressed as [2] hΔr 2 τi hV 2R iτ2 ;
(14)
and the flow, ΦR , as ΦR ∝ P m hV 2R i; where
hV 2R i
(15)
is the second moment of the scatterers speed.
where DL is the effective diffusion coefficient and τL is the time scale of the randomization of the scatterers velocity. 4. Mixed Model
Another model trying to better reproduce the physical/ physiological reality and allowing one to simultaneously take into account Brownian and random flow has been proposed in Ref. [17] and is expressed as hΔr 2 τi 6DM τ hV 2M iτ2 : In this case the flow, ΦM , is estimated by ΦM ∝ P m D M ; where DM is the diffusion coefficient and moment of the scatterers speed.
(19) hV 2M i
In Ref. [15] it is explained that the Brownian model (Section 2.B.1) assumes that the ballistic to random walk transition in the diffusion of the moving scatterers occurs at time scales shorter than that probed by DCS. Considering that there are no data to support this assumption, a more general formulation releasing the assumption has been proposed for hΔr 2 τi, i.e., [15,16] −ττ 2 ; (16) hΔr τi 6DL τ − τL 1 − e L and the flow, ΦL , is expressed as ΦL ∝ P m D L ;
(17)
is the second
5. “Exact” Model
The previous four models are typically utilized in real experimental conditions, where it is not known if the underlining microvascular geometry and blood flow speed distribution actually correspond to the chosen model. In the present case, the exact speed distribution is known [Eq. (3)]. Thus, one can derive an “exact” model corresponding to this specific distribution. The aim here is to further test the common mathematical techniques utilized to derive all the previous four models. By exactly following the same mathematical procedure utilized to derive the previous four models, it is straightforward to demonstrate [18] that if a semi-infinite medium has a speed distribution for the scatterers equal to Eq. (3), Eq. (10) always holds and in this case hΔr 2 τi v2B ‖⃗v t ‖2 τ2 h‖⃗v ‖2 iτ2 :
(20)
Equation (20) is similar to Eq. (14), and for this reason we expect hV 2R i and h‖⃗v ‖2 i to be equal if fitted to the present MCan data. It is essential to note again at this point that models typically utilized in DCS are P 1 approximations of the correlation transport equation [6,18]. This means that some parameters appearing in the final equations must satisfy some mathematical constraints. This is the case for τ, which cannot be larger than the critical constant τc (if not applicable, the model brakes). For the “exact” model, τc is equal to −1∕2 1 ; (21) τc k2o v 2B ‖⃗v t ‖2 3 where k o 2πn∕λ; and the constraint may be expressed as τ1.05 < τc :
3. Langevin Model
(18)
(22)
Intuitively, this means that the g 2 τ curve (“exact” model) must roughly reach the value 1 in a time τc (see Fig. 3), determined by the parameters v B and ‖⃗v t ‖. If this is not the case, the model brakes. Note also that Eqs. (14) and (20) are formally equivalents. For this reason, in the concluding Section 4 we will use the terms random or “exact” indistinctly. C. Simulations and Fitting Procedure
The reference data [i.e., g 2MC τ curves] have been generated for a semi-infinite medium with refractive index n 1.4. This n value implies A ≈ 1.96; precisely calculated “on flight” during the simulations (see Section 2.B). The optical fibers (normal to
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
5323
the medium surface) spacing was set to ρ 30 mm and the CW laser wavelength to λ 785 nm. Two sets of optical parameters were utilized for the simulations, denoted data set 1 and data set 2, generating the reference g 2MC τ curves. To facilitate the reading, consider that the set of parameters containing μa 0.005 mm−1 and μs0 1 mm−1 (see the following subsections) are reasonably within the range of the typical values for a bone tissue; they have been called “bone” [19]. For the same reason, the set of parameters containing μa 0.025 mm−1 and μs0 0.5 mm−1 have been called “muscle.” 1. Data Set 1: “Bone”
All combinations of the following parameters were considered: μa 0.005 mm−1 , μs0 1 mm−1 , P m ∈ f0.01; 0.025; 0.05; 0.075; 0.1g [8,20], v B ∈ f0.025; 0.05; 0.075; 1; 1.25; 1.5; 1.75; 2g mm s−1 , ‖⃗v t ‖ ∈ f0; 0.025; 0.05; 0.075; 1; 1.25; 1.5; 1.75; 2g mm s−1 , g 0.9, and λ 785 nm. This gave a total of 360 reference g 2MC τ curves. 2. Data Set 2: “Muscle”
The same as in Section 3.C.1 with the exception that μa 0.025 mm−1 and μs0 0.5 mm−1 . This gave a total of 360 reference g 2MC τ curves. 3. Fitting Procedure and Flow Assessment
The reference g 2MC τ curves were fitted using the four investigated DCS models, represented by g 2 τ and the relative Eqs. (12), (14), (16), and (18). A trust-region-reflective algorithm (MATLAB language function), minimizing the squared l 2 -norm of the residuals, e 2 min ‖g 2MC τ − g 2 τ‖22 ; C
(23)
was applied to perform the fitting; where C refers as usual to a particular set of fitting parameters, depending on the chosen DCS model. For the Brownian model C ∈ fP m DB g (to be seen as a single fitting parameter), for the random model C ∈ fP m hV 2R ig (one fitting parameter), for the Langevin model C ∈ fP m DL ; τL g (two fitting parameters), and for the mixed model C ∈ fP m DM ; P m hV 2M ig (two fitting parameters). It is important to note that when analyzing a real experiment, P m can never be isolated as a single fitting parameter. The above fitting parameters gave the ΦB , ΦR , ΦL , and ΦM estimators of Φ (in arbitrary units) for different physiological conditions and tissues. It is clear that the resulting estimators, obtained using Eqs. (12), (14), (16), and (18), do not strictly correspond to a flow, but this is actually the way the flow is estimated in DCS. This is a further reason explaining why it is so important to test the models on situations where the flow values are known. For explanatory purposes, the “exact” model was also fit to the reference g 2MC τ curves, but in this case P m was considered as known (this is not possible in a real experiment), and h‖⃗v ‖2 i was taken as a fitting parameter. By knowing the v B and ‖⃗v t ‖ values of the simulations, this approach allowed us to investigate if the “exact” model actually reproduces the relationship v 2B ‖⃗v t ‖2 h‖⃗v ‖2 i [see Eq. (20)], and thus demonstrate its reliability.
Fig. 1. Number of simulations, nsim , that have a percentage P int of photons that have interacted at least once with a moving scatterer before reaching the photodetector.
3. RESULTS In the present contribution, we have simulated 720 reference g 2MC τ curves, i.e., 360 simulations for “bone” and 360 simulations for “muscle.” From Fig. 1 one can see that all the 360 simulations for “bone” tissue have P int ≈ 100% of the photons reaching the photodetector that have interacted at least once with a moving scatterer. This, independently of the chosen simulation parameters. However, 144 simulations for “muscle” tissue have P int values slightly below 100%. This means that some photons have not yet interacted with a moving scatterer when reaching the photodetector. For the sake of completeness, in Fig. 2 is shown an example comparing P MC Δω and g 2MC τ obtained by MCan with the relative simulations obtained by a classical Monte Carlo (MC) method [21]. As already shown in a previous study, MCan and MC produce the same results [8]. However, it is important to note that the MC computing time for 720 P MC Δω spectra would be of the order of 1340 days while for MCan only 24 h were necessary.
Fig. 2. Example of P MC Δω and g 2MC τ computed both by classical MC and MCan . The simulation parameters were: μa 0.025 mm−1 , μs0 0.5 mm−1 , P m 0.05, v B 1, mm s−1 , ‖⃗v t ‖ 2 mm s−1 , g 0.9, n 1.4, and λ 785 nm.
5324
Vol. 54, No. 17 / June 10 2015 / Applied Optics
Research Article
Fig. 3. Example of g 2MC τ curve (black line) with the relative g 2 τ fittings (remaining colors) for the four different models. The dashed vertical lines represent the position of τ1.05 [Eq. (7)] and τc [Eq. (21)] on the abscissa.
In Fig. 3, a typical g 2MC τ curve with the relative g 2 τ fitting results for the four models (see Section 2.B) are presented. By visual inspection it clearly appears that the Brownian model is the worst one. To give a more quantitative estimate of the quality of the fittings for the four investigated models, in Fig. 4 are reported the mean and standard deviation (720 curves) of the squared 2-norm of the residuals, e 2 [Eq. (23)], for each model. As expected (see Section 4), the Brownian model is the worst model and, in general (not shown), it reproduces fitted curves similar to the one shown in Fig. 3 (blue dotted line). In Fig. 5 are shown the different flow estimators ΦB , ΦR , ΦL , and ΦM , compared to the exact value Φ [Eq. (6)]. It immediately appears that none of the models correctly takes into account the dependence on P m . A methodical P m -induced shift in the curves is observed. It is worth remarking that the Brownian model, ΦB , curiously generates reasonably “linear” data for a given tissue. However, considering that the Brownian fitted curves appear in general as in Fig. 3, i.e., the model does not fit at all the curves, this result must be considered as a pure mathematical artifact. Moreover, the
Fig. 4. Squared 2-norm of the residuals, e 2 [Eq. (23)], for the four fitted models (360 simulations for “bone” and 360 simulations for “muscle”). Vertical bars are standard deviations.
Fig. 5. Estimation of the scatterers’ flow, Φ, assessed by means of the four investigated models. P m 0.01 (black), P m 0.025 (red), P m 0.05 (blue), P m 0.075 (green), and P m 0.1 (magenta).
Brownian model is not able to correct for μa and μs0 differences [compare ordinate values, for the same Φ, in Figs. 5(a) and 5(b)]. Clearly, ΦL and ΦM are also not good estimators, even if they have the smallest e 2 , because the data are scattered on too large a region compared with the others models (good fitting, bad Φ estimator). In conclusion, the best estimator is represented by the random model, even if ΦR is not perfectly linearly related to Φ, and that the influence of P m is not correctly taken into account (shift in the colored curves). This is not astonishing because the random model Eq. (14) is in practice formally equal to the “exact” model Eq. (20). Let us now investigate if the “exact” model itself has not some unwanted hidden defect. While in a real experimental setup hV 2R i cannot be computed independently of P m , in the present study this has been done because P m was obviously known because it was utilized to generate the reference MCan curves. In this manner it has been possible to compare h‖⃗v ‖2 i obtained with the simulation parameters for MCan [Eq. (5)] and h‖⃗v ‖2 i estimated by fitting; theoretically, they should give exactly the same values (see Section 2.B.5). In Fig. 6 one can see that hV 2R i values for “bone” follows the identity line quite well, and that the error depending on P m is not visible. Thus, the “exact” model seems to be able to reasonably generate for this set of data the exact estimators. However, the results for “muscle,” even if linear, are less precise and the P m -induced shifts are always present for small P m values. It is worth noting that in this case the deviation from the identity line observed in Fig. 6(b) can be in part explained by plotting only the points
Research Article
Fig. 6. Panels (a) and (b) comparison of hV 2R i, assessed by means of the random model, with the exact value h‖⃗v ‖2 i Eq. (5), and supposing P m is known. This case does not correspond to the experimental reality where P m and hV 2R i cannot be separated (see manuscript). P m 0.01 (black), P m 0.025 (red), P m 0.05 (blue), P m 0.075 (green), and P m 0.1 (magenta). The dashed line is the identity line. Panels (c) and (d) similar to (a) and (b), but where only the points satisfying the constraint τ1.05 < τc (see Section 4) are represented. The data point represented by the red star is computed from the data in Fig. 3 (random model).
that are obtained from g 2MC τ curves satisfying the inequality τ1.05 < τc (see Section 2.B.5); a necessary condition allowing one to derive the model represented by Eq. (10). Figures 6(c) and 6(d) clearly show that if τ1.05 < τc , then the fitted hV 2R i are better distributed on the identity line. The red star appearing in Fig. 6(d) represents a particular hV 2R i value computed from the example (random model) of Fig. 3, where the constraint τ1.05 < τc is also valid. 4. DISCUSSION AND CONCLUSIONS In the present contribution we have shown that the MCan approach allows one to investigate if a classical DCS model is able to reasonably estimate scatters flow, ϕ (e.g., blood flow), in arbitrary units. It has been shown that the random model, among the proposed four models, is the best model allowing one to describe a flow that has a Maxwell–Boltzmann distribution for the moving scatterers speed. Even if the random or “exact” models have been built ad hoc for the Maxwell–Boltzmann speed distribution, this result remains nontrivial. In fact, we have seen (Section 2.B.5) that the derivation of g 2 τ implies a well-known specific approximation (oτ3 , [18]), on the field-autocorrelation function for a single scattering event, that obviously do not appear in the MCan -generated g 2MC τ curves. This means that g 2 τ and g 2MC τ are not necessarily the same, even for the “exact” model. This problem has been highlighted in Figs. 6(a) and 6(b) where, even after having taken into account the τ1.05 < τc constraint, some differences with the identity line still remain. The remaining inaccuracies cannot ensue from the reference data noise, because MCan is noise free. The inaccuracies are probably due, as highlighted above, to the too low order approximation utilized when deriving Eq. (20). In practice, this means that if the v 2B or ‖⃗v t ‖ (or a combination) values are too
Vol. 54, No. 17 / June 10 2015 / Applied Optics
5325
large, the random (or “exact”) model may be not valid. This point will certainly be a matter for future investigations. In the present work, it has also been shown [Figs. 5(c) and 5(d)] that when estimating ΦR , an inexact dependence on P m is observed, and that the ΦR − Φ relationship is not perfectly linear. One of the reasons might be the fact that the estimator ΦR contains the second moment of the speed and not the speed itself [see Eq. (15)]. However, the positive point is that for given real tissue (in vivo), the parameter ΦR (proportional to the number of moving scatterers in the observed tissue) is not subject to such large P m changes. Thus, ΦR probably moves in a neighborhood along one given colored line in Figs. 5(c) and 5(d). This possibly explains also why DCS appears to give reasonable results in real experimental setups, even in the presence of this mathematical problem. In the past it has been noted that, experimentally, the photons detected by the DCS instrumentation, scatter not only on moving scatterers, but also on a fixed background [18]. This fact is known to introduce a distortion in the g 2 τ function (a phenomenon not described by Monte Carlo simulations). It has been hypothesized that the good correlation of the DCS models with the experimental measurements, even in the presence of fixed background scattering, may be due to the fact that most of the photons interact at least once with a moving body [2]. By inspecting Fig. 1 this seems to be true because of the high P int values. However, this remain a necessary but not sufficient condition, and other investigations must be performed on this topic. Summarizing the main concepts, we can say that the proposed method allows one to test DCS analytical models on synthetic reference data, g 2MC τ, precisely reproducing a perfectly known moving scatterers speed/flow distribution. In fact, we have seen that DCS analytical models typically utilized in biomedical optics are only approximated solutions. This means that DCS analytical models may not exactly reproduce a fortiori the synthetic data, g 2MC τ, or generate the correct desired parameter values. Thus, the proposed method allows one to study, e.g., the domain of validity of the DCS models (e.g., speed range where they remain valid) or their robustness on synthetic data where the speed distribution is not exactly the one considered by the DCS model. Conditions that may easily be present in the real experimental setting. Thus, the proposed approach is to be considered before any other test such as the study of the influence of noise or hardware artifacts [3,22,23]. If the proposed test on g 2MC τ reference data fails, e.g., for a particular range of speed parameters, then it is theoretically pointless to apply the model in this situation. This, because there is an obvious intrinsic problem in the DCS model, is independent from noise or hardware artifacts. In this case, the possibility of a further improvement of the model must be considered. It must be noted that in vivo the Brownian model better fits the experimental data than the random model [2,24–27]. For demonstrative purposes in the present contribution, the model utilized for the moving scatterers speed was a Maxwell– Boltzmann distribution (with a bulk drift). Thus, it is mathematically expected that the best DCS model (without considering the observed approximation errors) is in this case
5326
Research Article
Vol. 54, No. 17 / June 10 2015 / Applied Optics
the (“exact”) random model (see Section 2.B.5). Thus, there is no inconsistency between the present results and the experimental one. These results only show that the classical Brownian model is not able to extract the right physiological parameters if the speed distribution is Maxwell–Boltzmann, and that the Brownian model is not robust enough in this case. Considering that the probability density function of the speed in real tissue is not exactly known, future tests of the robustness of the Brownian model for other speed distributions will certainly be of interest. The Langevin model deserves a last comment. In fact, it represent a clear example where a perfect fitting does not mean that the derived parameters of interest are necessarily exact (see Figs. 4 and 5). This is probably due to the cross talk existing between DL and τL . This motivates again the interest of the tests proposed in the present contribution. In conclusion, we have demonstrated that the MCan -based approach may be used to check the reliability of a given DCS model to generate the suitable Φ estimator, when applied to a given moving scatteres speed distribution. Moreover, it must be noted that MCan can be adapted to other speed distributions, simulating in this way other experimental conditions (e.g., satisfying the Brownian DCS model). This will allow, for example, one to test in the future the robustness of a given DCS model on flow distributions that do not exactly represent the conditions for which the model was built (this is probably what happens in a real experiment). We hope that the present work will contribute to gaining a further understanding about the DCS technique. REFERENCES 1. D. A. Boas, L. E. Campbell, and A. G. Yodh, “Scattering and imaging with diffusing temporal field correlations,” Phys. Rev. Lett. 75, 1855– 1858 (1995). 2. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73, 076701 (2010). 3. Y. Shang, T. Li, L. Chen, Y. Lin, M. Toborek, and G. Yu, “Extraction of diffuse correlation spectroscopy flow index by integration of Nth-order linear model with Monte Carlo simulation,” Appl. Phys. Lett. 104, 193703 (2014). 4. Y. Shang and G. Yu, “A Nth-order linear algorithm for extracting diffuse correlation spectroscopy blood flow indices in heterogeneous tissues,” Appl. Phys. Lett. 105, 133702 (2014). 5. W. B. Baker, A. B. Parthasarathy, D. R. Busch, R. C. Mesquita, J. H. Greenberg, and A. G. Yodh, “Modified Beer–Lambert law for blood flow,” Biomed. Opt. Express 5, 4053–4075 (2014). 6. B. J. Ackerson, R. L. Dougherty, N. M. Reguigui, and U. Nobbman, “Correlation transfer: application of radiative transfer solution methods to photon correlation problems,” J. Thermophys. Heat Transfer 6, 577–588 (1992). 7. T. Binzoni, A. Torricelli, R. Giust, B. Sanguinetti, P. Bernhard, and L. Spinelli, “Bone tissue phantoms for optical flowmeters at large interoptode spacing generated by 3D-stereolithography,” Biomed. Opt. Express 5, 2715–2725 (2014). 8. T. Binzoni and F. Martelli, “Random numbers free analytical implementation of Monte Carlo for laser-Doppler flowmetry at large
9.
10. 11.
12.
13.
14.
15.
16. 17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
interoptode spacing: application to human bone tissue,” Appl. Opt. 54, 2400–2406 (2015). H. Z. Cummins and H. L. Swinney, “Light beating spectroscopy,” in Progress in Optics, E. Wolf, ed. (North-Holland, 1970), Vol. 8, pp. 133–200. R. Bonner and R. Nossal, “Model for laser Doppler measurements of blood flow in tissue,” Appl. Opt. 20, 2097–2107 (1981). J. Zhong, A. M. Seifalian, G. E. Salerud, and G. E. Nilsson, “A mathematical analysis on the biological zero problem in laser Doppler flowmetry,” IEEE Trans. Biomed. Eng. 45, 354–364 (1998). T. Binzoni and D. Van De Ville, “Full-field laser-Doppler imaging and its physiological significance for tissue blood perfusion,” Phys. Med. Biol. 53, 6673–6694 (2008). C. Cheung, J. P. Culver, K. Takahashi, J. H. Greenberg, and A. G. Yodh, “In vivo cerebrovascular measurement combining diffuse near-infrared absorption and correlation spectroscopies,” Phys. Med. Biol. 46, 2053–2065 (2001). D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). S. A. Carp, N. Roche-Labarbe, M. A. Franceschini, V. J. Srinivasan, S. Sakadžić, and D. A. Boas, “Due to intravascular multiple sequential scattering, diffuse correlation spectroscopy of tissue primarily measures relative red blood cell motion within vessels,” Biomed. Opt. Express 2, 2047–2054 (2011). B. J. Berne and R. Pecora, Dynamic Light Scattering: with Applications to Chemistry, Biology, and Physics (Dover, 2000). H. M. Varma, C. P. Valdes, A. K. Kristoffersen, J. P. Culver, and T. Durduran, “Speckle contrast optical tomography: a new method for deep tissue three-dimensional tomography of blood flow,” Biomed. Opt. Express 5, 1275–1289 (2014). D. A. Boas, “Diffuse photon probes of structural and dynamical properties of turbid media: theory and biomedical applications,” Ph.D. thesis (University of Pennsylvania, 1996). P. Farzam, P. Zirak, T. Binzoni, and T. Durduran, “Pulsatile and steady-state hemodynamics of the human patella bone by diffuse optical spectroscopy,” Physiol. Meas. 34, 839–857 (2013). A. Kienle, “Non-invasive determination of muscle blood flow in the extremities from laser Doppler spectra,” Phys. Med. Biol. 46, 1231– 1244 (2001). A. Sassaroli and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A 29, 2110–2117 (2012). L. Dong, L. He, Y. Lin, Y. Shang, and G. Yu, “Simultaneously extracting multiple parameters via fitting one single autocorrelation function curve in diffuse correlation spectroscopy,” IEEE Trans. Biomed. Eng. 60, 361–368 (2013). C. Zhou, G. Yu, D. Furuya, J. Greenberg, A. Yodh, and T. Durduran, “Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain,” Opt. Express 14, 1125– 1144 (2006). R. C. Mesquita, T. Durduran, G. Yu, E. M. Buckley, M. N. Kim, C. Zhou, R. Choe, U. Sunar, and A. G. Yodh, “Direct measurement of tissue blood flow and metabolism with diffuse optics,” Philos. Trans. R. Soc. A 369, 4390–4406 (2011). G. Yu, T. Durduran, C. Zhou, R. Cheng, and A. G. Yodh, “Nearinfrared diffuse correlation spectroscopy for assessment of tissue blood flow,” in Handbook of Biomedical Optics, D. A. Boas, C. Pitris, and N. Ramanujam, eds. (CRC Press, 2011). J. Dong, R. Bi, J. H. Ho, P. S. Thong, K. C. Soo, and K. Lee, “Diffuse correlation spectroscopy with a fast Fourier transform-based software autocorrelator,” J. Biomed. Opt. 17, 0970041 (2012). T. Durduran and A. G. Yodh, “Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement,” Neuroimage 85, 51–63 (2014).
