Experimental tests of indicators for the degree of validness of the diffusion approximation

The diffusion approximation has been one of the central topics in near-infrared spectroscopy (NIRS). When NIRS measurements are analyzed by the diffusion theory, the measurements must be performed in the diffusive regime. However, since most of past researches have focused on theoretical or qualitative nature of the diffusion approximation, it is not easy to know if each measurement is designed in the diffusive regime. In this paper, we consider the diffusion approximation quantitatively and propose indicators that quantify the degree of validness of the diffusion approximation. The difference between the measurement and diffusion theory can be evaluated with the χ 2 value, ℓ 1 and ℓ 2 norms, and Kullback-Leibler divergence. We conduct a liquid phantom experiment to test the proposed χ 2 value. Moreover, the χ 2 value is further investigated by Monte Carlo simulations. We find the χ 2 value becomes significantly large when measurements are performed in the nondiffusive or transport regime. The proposed indicators similarly work. In particular, the χ 2 value is shown to work as an indicator which evaluates the degree of validness of the diffusion approximation. These indicators are general and can be used for different numerical, experimental, and clinical measurements in NIRS.


Introduction
Since near-infrared light is absorbed and scattered in biological tissue, the light detected on the surface has the information on the medium in which light propagates. Hence we can access optical properties of the medium such as the absorption and scattering coefficients by illuminating the medium by near-infrared light. Nearinfrared spectroscopy (NIRS) including optical tomography often relies on the diffusion approximation to the radiative transport equation. When optical properties of biological tissue are retrieved using the diffusion theory, measurements must be conducted in the diffusive regime, in which the source and detector are well separated and the absorption coefficient μ a is sufficiently smaller than the reduced scattering coefficient m ¢ s . Although it has been one of the central topics in NIRS to identify the diffusive regime, most researches have qualitatively studied the diffusion approximation. Hence, in this paper, we revisit the diffusion approximation quantitatively and provide an indicator which helps design measurements in the diffusive regime.
The transition from the transport regime to the diffusive regime takes place asymptotically when the propagation distance of light is large and  m m ¢ a s . That is, there is no clear boundary between the two regimes. However, there is a practical need of drawing the border. In NIRS, optical properties of biological tissue are determined by the comparison between measurements and simulational results. Since the computational cost of the radiative transport equation is quite often too expensive, the diffusion equation should be used if measurements are performed in the diffusive regime. Past researches on the transition from the transport regime to diffusive regime include the following studies.
In [1], the ratio of the fluence rate and its value from the diffusion approximation was plotted against the distance between the source and detector for different μ a values. The plot shows that the diffusion approximation works well when m m ¢ = 1 100 a s and the distance is larger than m ¢ 10 s and the approximation becomes worse when the absorption exceeds m m ¢ = 1 25 a s . Although such tendencies can be read from such qualitative study, it is not easy to capture how satisfactory the diffusion approximation holds for given μ a and source-detector distance. Here we will present an indicator which shows the degree of validness of the diffusion approximation.
In [2], time-resolved measurements were performed in the transmission geometry and measured temporal profiles were compared with the prediction by the diffusion theory. They concluded that the diffusion approximation breaks when the width of the slab becomes smaller than 10ℓ * (ℓ * is the transport mean free path). In [3], solutions to the time-independent diffusion equation were compared with Monte Carlo solutions in the reflection geometry for a semi-infinite medium. From the comparison, the criterion m m < ¢ 4 a s can be read as a condition that the diffusion approximation holds. Moreover they found that the diffusion results are more accurate for small g. By investigating the Brownian motion using a low coherence interferometry, it was experimentally shown that the transition to the diffusive regime depends on the anisotropic factor g even when m ¢ s is unchanged [4]. The breakdown threshold 8ℓ * was suggested in [5]. In their study, this breakdown threshold was obtained through time-resolved experiments in the transmission geometry with a sample of various thicknesses. By observing temporal profiles of transmitted light, slabs of width greater than m ¢ 8 s were defined to be in the diffusive regime [6,7]. In [8], the condition m m < ¢ 10 a s was proposed. The limits of the diffusion approximation was considered with a two-dimensional numerical phantom of the human brain and a threedimensional cubic numerical phantom by the comparison between the diffusion equation and radiative transport equation calculated by finite-difference scheme [1,9]. Their calculation shows that the fluence rate by the diffusion approximation becomes half of its true value when m m = ¢ 2 a s and the distance is 8 in the unit of m ¢ 1 s . It was reported that the deviation of the fluence rate increases when the anisotropic factor g (introduced below) is close to 1, which is typical in biological tissue [10]. The deviation of the dynamic intensity-intensity correlation function from the prediction of the diffusion theory was experimentally observed [11].
The extrapolated boundary condition has often been used to further approximate the diffusion equation. It is known that those approximate solutions are quite precise in the diffusive regime. The reflectance and transmittance were derived as functions of time from the diffusion equation solved using the extrapolated boundary condition [12]. The prediction from those concise formulae agreed well with Monte Carlo simulations. In [13], the validity of the diffusion model was confirmed by experiments and Monte Carlo simulations. In [14], the validity of the diffusion model was studied in frequency domain using the spherical and cylindrical geometries in addition to the slab geometry. The effect of different approximations to the diffusion equation was studied [15,16]. The diffusion approximation with the extrapolated boundary condition was tested using Monte Carlo in the slab geometry [17].
It was found that the quality of the diffusion approximation depends on the refractive index [13,18]. Results for different index mismatches were reported [19]. Moreover, the validity of the diffusion approximation on the surface was studied [20].
The nondiffusive regime was studied by using other approximations to the radiative transport equation. When m m > ¢ a s , a clear difference was observed between solutions to the diffusion equation and telegrapher equation, which takes the ballistic component of photon propagation at short times into account [21]. The diffusion approximation was compared to the P 1 and P 3 approximations [22].
The light propagation is said to be in the diffusive regime if the solution to the radiative transport equation and the solution to the diffusion equation match. However, there is no unique way of setting the border of the diffusive regime because the change from the transport regime to diffusive regime takes place gradually. To draw a line between two regimes, a criterion must be set. In this paper, a constant a (see below) is used as the threshold. The solution to the radiative transport equation can be obtained from measurements since infrared light in biological tissue is governed by the radiative transport equation. In this case, the threshold a must be set in such a way that the determined boundary of the diffusive regime is not affected by measurement noise or error.
In this paper, we propose to use an indicator such as the χ 2 value to identify the diffusive regime. Thus, the degree of validness of diffusion approximation is given as a χ 2 value. Although the change from the nondiffusive or transport regime to the diffusive regime takes place asymptotically, we find that χ 2 values rather rapidly change when measurements shift from one regime to the other. The behavior of the χ 2 value is tested both by liquid phantom experiments and Monte Carlo simulations. By Monte Carlo simulation, it is confirmed that other indicators ℓ 1 and ℓ 2 norms, and the Kullback-Leibler (KL) divergence work as well.

Diffusion theory
Here, the diffusion theory is developed for a three-dimensional semi-infinite medium. In addition to the reflection on the boundary, we take into account the instrument response function, diameters of the optical fibers, and numerical aperture in the diffusion approximation. The approximate form Φ DA (t) of the detected light Φ(ρ d , t) (defined below) will be derived.

The radiative transport equation
Let us consider the three-dimensional half space . We consider the specific intensity ----(ˆ) I t r s , , of near-infrared light at position r = (x, y, z) in directionˆ( J = s sin ) j J j J cos , sin sin , cos at time t. Here, the polar angle ϑ and azimuthal angle j vary in 0 ϑ π and 0 j < 2π, respectively. When absorption coefficient μ a and scattering coefficient μ s are positive constants, the specific intensity (ˆ) I t r s , , obeys the following radiative transport equation. where c > 0 is the speed of light in the medium and (ˆˆ) ¢ p s, s is the scattering phase function, which is normalized as s , s . The anisotropic factor g is introduced as The reduced scattering coefficient is given by There is the source term (ˆ) q t r, s, on the righthand side of (1).
In addition to the initial condition (ˆ) = I r s , , 0 0, the boundary condition is imposed as . Let r s = (ρ s , 0 + ), where ρ s = (x s , y s ), be the position of the optical fiber for the incident beam on the surface. We set ρ s = 0. We write the source term (ˆ) q t r, s, as is the temporal profile of the incident beam. In our measurement in the reflection geometry, the detector is placed on the boundary at r  t denote the response function. Considering the diameter ( ) r  , r d and numerical aperture (ˆ)  s of the optical fiber for the detector at r d , the measured quantity Φ(ρ d , t) is expressed as If the detector radius ρ fiber is small, we can put ( The function ( )  t below corresponds to the instrument response function.
r, s; r , z; s RTE becomes the Green's function for the radiative transport equation in the half space. Let us define Therefore we obtain s sin , r , 0, s; r , z; . The above Φ(ρ d , t) corresponds to the measurement by the TRS-20, which we denote by Φ TRS (ρ d , t).

Diffusion approximation
Let us calculate Φ(ρ d , t) using the diffusion approximation. In the spirit of the diffusion approximation, we write the specific intensity as By substituting the above expression (13) in (1) after some calculations [23], we obtain where the diffusion coefficient is given by [24] ( ) The dependence of the diffusion coefficient on parameters in the radiative transport equation has been studied. Indeed, a is derived from the straightforward calculation. However, the form ( ) m ¢ 1 3 s was suggested by more detailed studies [25][26][27]24]. In Appendix A, we develop an alternative approach to the diffusion approximation. The diffusion coefficient (15 The source term is obtained as Here we introduced The current J(r, t) is obtained as Let us consider the boundary condition for the diffusion equation (14). With the diffuse surface reflection, we obtain the following Robin boundary condition.
Moreover the initial condition is given by u(r, 0) = 0. Let us calculate Then we have within the diffusion approximation, , Let us begin by considering the diffusion approximation to the Green's function G RTE . In the radiative transport equation (1), we set , and the source term for the diffusion equation is obtained as  where we used integration by parts and the reciprocal property ( r , r; . Hence by making use of the Robin boundary condition, Therefore we arrive at We can obtain the Green's function as  . We note that G = 0 when t < 0. The expression (29) is obtained by straightforward calculation [30].
We can write is not required because the constant factors on the right-hand side of (28) can be absorbed by the constant ¢ C 1 . That is, the following relation holds in the diffusive regime.  Figure 1 shows the experimental setup. A liquid phantom was prepared in a bucket of radius 8-10 cm and height 20 cm. The 10% Intralipos (Otsuka, Tokyo, Japan) and green brown ink (Chugai Kasei, Tokyo, Japan) were gradually added to water (initially 3 L) to control the scattering and absorption, respectively. The incident and detecting optical fibers from the TRS-20 were attached to a holder. The separation of the fibers was 2 cm, i.e., ρ d = 20 mm. We set the holder on the surface of the liquid phantom and supported optical fibers by using two arm-type clamps. The liquid phantom was stirred with a magnetic stirrer (SW-030, Nissin Science, Tokyo, Japan) except for during measurement. We confirmed that the bucket is deep enough and the liquid phantom can be regarded as a semi-infinite medium even for smallest m m ¢ , a s used in the experiments. In the diffusive regime for large m ¢ s , optical properties of the liquid phantom were estimated by the TRS-20, which obtains μ a and m ¢ s using the Levenberg-Marquardt method by minimizing the difference between the measured temporal profile and the time-resolved reflectance calculated from the diffusion theory (see [31] for details). We confirmed that estimated μ a values by the TRS-20 were consistent with values from spectrophotometric measurements (UV-3100, Shimadzu, Kyoto, Japan). The absorption coefficient μ a of the liquid phantom can be kept the same for different m ¢ s by adjusting the ink so that the ratio of the ink to the total amount of the liquid phantom is fixed.

Monte Carlo simulations
In addition to the liquid phantom experiment, we employ Monte Carlo simulation, in which optical properties of the medium can be precisely specified. The validity of Monte Carlo simulations for a liquid phantom was confirmed in [32]. Our Monte Carlo simulation was implemented using the variance reduction technique [33,34]. While photons propagate in the medium, they are scattered according to the scattering phase function at each scattering point. The Henyey-Greenstein model [35] was used. The anisotropy factor was set to g = 0.9 for μ a = 0.06 mm −1 and μ a = 0.1 mm −1 , and g = 0.33 otherwise. The refractive index was set to n=1.33. On the boundary at z = 0, the reflection and refraction caused by the refractive index mismatch were considered.
The diameter and numerical aperture of the incident fiber at the origin are 0.2 mm and 0.250, respectively. The detecting optical fiber was placed at a distance of 20 mm from the source (that is, ρ d = 20 mm). For the detector, the diameter and numerical aperture were set to 3 mm and 0.260, respectively.
The weight of detected photon packets was accumulated into time-resolved arrays for absorption coefficients. The duration of an array element was 10 ps, corresponding to the temporal resolution of the TRS instrument used in our experiment. In each Monte Carlo simulation, 10 8 photons were launched.
In Monte Carlo simulations, The constant C 2 is obtained as  There are different ways of measuring the difference between two functions. To evaluate the difference between the measured time-resolved reflectance, which is the solution to the radiative transport equation, and the temporal profile by the diffusion approximation, we introduce the χ 2 value as The lower and upper limits k a , k b are chosen such that ( )  F t k DA for k = 0, K, k a − 1 and k = k b + 1, K, M − 1 are less than one-tenth of the peak height of ( ) . The integers k a , k b will be chosen such that values which are more than 10% of the peak height are considered in the sum in (39).
In addition to the χ 2 value, we propose the ℓ 1 norm, ℓ 2 norm, and KL divergence. They are given by , respectively. Here, . These indicators behave similarly if scaled properly (see figure 8 below). Hence we will focus on the χ 2 indicator in this paper.
A few remarks are necessary: (i) We note that the χ 2 indicator introduced above has nothing to do with the χ 2 test in statistics. The χ 2 value in (39) is used to measure the distance between the solution to the radiative transport equation and the solution to the diffusion equation. In this paper, we admit that experimentally obtained data obey the radiative transport equation and treat ( ) F t EXP as the solution to the radiative transport equation. This means that we do not argue how precisely measured experimental data are described by the radiative transport equation. (ii) The KL divergence sometimes becomes negative whereas the χ 2 value, ℓ 1 norm, and ℓ 2 norm are always nonnegative.
We separate obtained data for different μ a and m ¢ s into the diffusive and nondiffusive regimes using the χ 2 value (39). To this end, we set a threshold a. The measurement is identified as in the diffusive regime if χ 2 < a. Otherwise the measurement is considered to be in the nondiffusive or transport regime.  Figure 2 illustrates that the indicator takes large values when μ a is large.

Results
The diffusive and nondiffusive regimes are segmented in figure 3. The threshold a is set to 5. Since the χ 2 value rapidly grows as shown in figure 2 once performed measurements move from the diffusive regime to the transport regime, the same result is obtained if a = 4 or a = 6. Figure 4 shows time-resolved reflectances for μ a = 0.021 mm −1 (the top curve in figure 3). In Let us try other values of the threshold a and see how the defined diffusive regime is affected by a. In the left panel of figure 6, we set a = 20. Since the criterion a = 20 is loose, the defined diffusive regime is wider than the region defined in figure 5. That is, there are more blue open squares in the left panel of figure 6. In the right panel of figure 6, the a value is set to 5. This a = 5 is severe and measurements are classified in the nondiffusive regime unless ρ d and μ a are very small. When measurements in the non-diffusive regime are of interest, a should be taken rather large (e.g., a = 20), so that the points marked in red are surely in the non-diffisive regime. On the other hand, if measurements in the diffusive regime are of interest, a should be set to small values (e.g., a = 5), so that blue open squares certainly belong to the diffisive regime. Figure 7 shows results for Monte Carlo simulations with μ a = 0.01 mm −1 . All curves in figure 7 are plotted together in figure 8 after scaling. The ℓ 1 norm, ℓ 2 norm, and KL divergence are scaled as where r m = ¢ x d s and c y 2 is the χ 2 value. Here, y is either ℓ y 1, ℓ y 2, or y KL , which are defined in the same way as c y 2. Table 1 shows derivatives dy/dx (central differences) for the indicators χ 2 , ℓ 1 , ℓ 2 , and KL divergence, which are calculated from numerical values in figure 8. It is observed that derivatives take negative large values at r m ¢ = 6

Discussion
Let us compare  MC and  DA . We write

DA
We found A = 2.785 and C 3 = 1.3. The result that C 3 is slightly greater than 1 attributes to the fact that the diffusion approximation breaks on the boundary. In section 2, the source was placed on the surface when Φ DA (t) was calculated. In the diffusion approximation, quite often the source S(r, t) is placed inside the medium (z > 0)   about one transport-mean-free-path away although the source term (ˆ) q t r, s, represents the optical fiber on the boundary. This has an effect of improving the value of C 3 . To see the effect, let us consider the solution with the extrapolated boundary condition. In this case, the Green's function for the diffusion equation in the half space is obtained by using ( ) ¢ g z z t , ; EBC instead of ( ) ¢ g z z t , ; in (29). Here, . 42 Suppose the source is placed at r s = (x s , y s , d), where d is small (i.e., 2Ad = ct and 4DA 2 = ct). Then we have g EBC (0, d; t) = [1 + d/(2DA)]g EBC (0, 0; t). We find that m = ¢ d 0.871 s results in C 3 = 1.0 (the parameters, ρ d = 20 mm, ρ fiber = 1.5 mm, J = 0.26 max , and n = 1.33 are used). Thus, the optimal depth is about the transport mean free path.
In addition to values shown in figure 8, we further tested the symmetric χ 2 value, which is given by The obtained values were close to the χ 2 values. As is explained in section 2.2, values of m ¢ s for the liquid phantom were determined using about several measured points for which the source-detector distance is large enough ( > 20ℓ * ). However, since m ¢ s obtained by the TRS-20 has only about two significant digits, the error in estimating m ¢ s in the nondiffusive regime may not be negligible. To find clearer boundaries between the diffusive and nondiffusive regimes, parameter estimation based on the radiative transport equation is necessary.
In NIRS, the reduced χ 2 value was used to discuss the accuracy of optical parameters calculated by TRS [36]. Instead of evaluating statistical errors, the aim of this paper is to compare the solutions of the radiative transport equation and diffusion equation. In this paper, the χ 2 function was introduced solely as an indicator which takes large values in the nondiffusive regime.
The χ 2 value (39) depends on the choice of k a , k b . In this paper we took into account the values which are more than 10% of the peak height. In [37] and [38], the range between 80% of the peak on the leading edge and 20% on the falling edge was used. In [39], the range from 80% on the leading edge and 1% on the falling edge was used. In addition to 10%-10% cut for k a , k b , we tried 10%-50%, 50%-10%, 50%-50%, 40%-10%, 80%-20%, and 80%-20%, and confirmed that the qualitative conclusion remains the same. In order to capture the difference between ( ) F t EXP and Φ DA (t), we chose the range 10%-10%, which is rather wide. When Φ DA (t) is used  for parameter identification, however, a large value such as 80% should preferably be taken on the leading edge because the diffusion approximation easily breaks at early times.

Conclusion
In this paper, we have proposed indicators such as the χ 2 value (39) to judge if the measurement is in the diffusive regime or in the nondiffusive regime. By using the χ 2 value, the validness of the diffusion approximation can be quantitatively studied. Once the criterion a is given, the diffusive regime is defined by the χ 2 value.
The usefulness of these indicators can be read from figure 2. As is shown in Appendix A, the diffusion equation is obtained from the radiative transport equation when the small parameter ò goes to zero. In this sense, the change from the transport regime to diffusive regime is not abrupt but takes place asymptotically. However, figure 2 shows that the change is rather sharp when the degree of validness of the diffusion approximation is measured by the χ 2 indicator function.
As is seen in section 3, the defined diffusive regime depends on the value of the threshold a. When measurements require that the diffusion approximation severely holds, a should be set to a small value. Suppose we want to conduct measurements in the diffusive regime. In order to draw a line between the diffusive and nondiffusive regimes, a must be defined operationally. First we try a relatively large a, so that many points in the r m ¢  The left-hand side of (A.9) can be expressed as  (14) is derived. We note that the above equation does not have the source term because we dropped the source term in the radiative transport equation in (A.2) by focusing on large space and time.
In our derivation, the diffusion coefficient D in (15), which is independent of μ a , was obtained.  Table B2. TRS measurements for μ a = 0.021 mm −1 (ρ d = 20 mm). The unit of m¢ s is mm −1 . m¢ s r m¢