I. Introduction

The rapidly increasing use of light in diagnostic and therapeutic medicine has created the need to determine noninvasively the optical properties of living tissue. These properties may be of intrinsic interest; for example, changes in the optical absorption coefficient can be related to blood oxygenation and tissue metabolism.[1],[2] Alternatively, the optical properties may be required to calculate the distribution of light in the tissue during, for example, photodynamic therapy[3] or laser surgery.[4] Many papers have been devoted to the measurement of the optical properties of excised tissue specimens, and the detection of diffusely reflected or transmitted light has been studied as a means of monitoring blood oxygenation in vivo.[5] More recently, Groenhuis et al.,[6] Bonner et al.,[7] and Steinke and Shepherd[8] have studied the relationship between the absorption and scattering coefficients of tissue and the spatial dependence of diffuse reflectance near a finite light source. Such measurements may be useful for noninvasive optical characterization, but absolute measurements of reflectance must be made at several locations.

The temporal spreading of a short light pulse as it propagates through a scattering medium also contains information about the optical interaction coefficients. Such measurements have long been of interest in atmospheric research; for example, Weinman and Shipley[9] used the time dependence of a transmitted pulse to deduce the optical thickness of clouds. Theoretical studies of light pulse propagation in multiple scattering media based on the diffusion approximation have been published by Ishimaru[10] and Furutsu.[11] Shimizu et al.[12] measured the time resolved reflectance of a plane wave from suspensions of microspheres and suggested that this technique might be used to determine the optical properties of random media such as tissue.

We propose a practical in vivo technique in which a short light pulse is produced on the tissue surface by a small source, such as a laser beam or optical fiber, and a small detector some distance away is used to measure the time resolved pulse. Two geometries are considered: a semi-infinite medium which represents measurements made on the surface of a large volume and a finite slab representing reflectance or transmittance measurements on smaller organs. A simple model based on the time dependent diffusion equation is introduced which allows calculation of the pulse shape at any point on either face of the slab. The predictions of the model are shown to be in good agreement with Monte Carlo simulations performed by us and other investigators. We show that in certain conditions the absorption and scattering coefficients can be expressed as functions of simple descriptors of the pulse shape: the time of maximum signal and the decay constant at long times. The pulse shape recorded from the calf muscle of a human volunteer is accurately predicted by the model, and the derived absorption and scattering coefficients are similar to published values for excised muscle tissue. If economical pulsed light sources can be used, this method offers the potential for fast simple noninvasive determination of tissue optical properties.

II. Theory

The geometry of the problem is illustrated in Figs. 1(a) and 1(b). A narrow collimated pulsed light beam is normally incident on the surface of a semi-infinite [Fig. 1(a)] or finite [Fig. 1(b)] homogeneous tissue slab. We will assume that the diffuse photon fluence rate ϕ(r,t) satisfies the diffusion equation

where c is the speed of light in the tissue, D is the diffusion coefficient,

μ_a is the linear absorption coefficient, μ_s is the linear scattering coefficient, g is the mean cosine of the scattering angle, and S(r,t) is the photon source. This equation can be derived from the radiative transfer equation, and its application to light propagation in turbid media has been discussed by many investigators including Ishimaru.[13] The fluence rate can be accurately calculated using Eq. (1) if μ_a ≪ (1 − g)μ_s and if the point of interest is far from sources or boundaries. The first condition is generally true for soft tissues in the 650–1300 nm range. Strictly speaking, the second condition is violated in calculating the reflectance and transmittance, but even here we will show that accurate estimates of the relative (rather than absolute) quantities can be obtained.

For a short pulse from an isotropic point source, S(r,t) = δ(0,0), it may be shown[14] that in an infinite medium the solution of Eq. (1) is

We can use this Green’s function to solve the problem posed in Fig. 1(a) by making two further assumptions. First, we assume that all the incident photons are initially scattered at a depth

so that the actual source term becomes the simple delta function described above. This localization of the first interactions will not produce inaccuracies if we are interested in the fluence rate far from the source or at times long after pulse incidence. We must also specify a boundary condition at the tissue surface. Duderstadt and Hamilton[15] have shown that a useful approach is to set the diffuse fluence rate to zero at an extrapolated boundary some distance beyond the actual surface. For example, if the interface is between tissue and a nonscattering medium with the same index of refraction, this extrapolated boundary is located at z = −(0.7104)3 ·D.[13] Hemenger[16] has shown that a mismatch in the index of refraction at the surface, for example, an air–tissue interface, can be accounted for by changing this extrapolation length. For our application, where observations are made at a large distance compared with the extrapolation length, we have found that the pulse shape is insensitive to the exact location of the extrapolated boundary. Hence, for simplicity, our second assumption is that ϕ(r,t) = 0 on the physical boundary z = 0. As discussed by Eason et al.,[17] this boundary condition can be met by adding a negative or image source of photons to the infinite medium problem as shown in Fig. 1(a). The fluence rate per incident photon can then be written in cylindrical coordinates as the sum of contributions from the two sources:

We wish to know the number of photons reaching the surface per unit area per unit time |J(ρ,0,t)|, which can be calculated from Fick’s law[15]:

which leads to the final expression for the reflectance R(ρ,t):

For the case where ${\rho }^2\gg z_0^2$ we also note that

The observation that

leads to the suggestion recently made by Chance et al.[18] that the absorption coefficient of the tissue can be determined from the asymptotic slope of the log_eR(ρ,t) vs t curve. (The speed of light depends on the index of refraction n of the tissue which is known to a few percent. We assume that n = 1.4,[19] so that c = 0.214 mm ps⁻¹.)

The transport scattering coefficient (1 − g)μ_s can also be determined from the log_eR(ρ,t) vs t curve by noting that at t_max, the time of maximum detected signal, the slope is zero. Solving Eq. (8) yields the expression

Thus we see from Eqs. (9) and (10) that the optical properties of a semi-infinite slab of tissue could in principle be obtained by measuring the diffusely reflected light some distance from the source as a function of time. A superior signal-to-noise would be obtained by integrating the reflected light over some larger area. As an illustration of the potential of this method, we can integrate R(ρ,t) over the entire surface to obtain R(t):

This expression agrees with the prediction of Ito and Furutsu[20] that for a nonabsorbing medium the total diffuse reflectance should depend on t^−3/2.

An important question is how measurements of the diffusely reflected light are affected by a finite tissue volume. A related question is whether a simple model can yield results for the time resolved transmittance through a tissue volume. Such measurements have been made for rat brains by Delpy et al.[19] and have been shown to carry information about the tissue optical properties. We will consider the simplest case: a finite tissue slab of thickness d illustrated in Fig. 1(b). Here we have an additional boundary where the condition ϕ(ρ,d,t) = 0 must apply. This condition can be met by adding two sources near z = 2d as shown, but then the boundary condition at z = 0 is violated for t > 2d/c. Both boundary conditions can be met at all times only by adding an infinite number of dipole sources as shown in Fig. 1(b). In practice, the number of sources required depends on the optical properties of the slab and the maximum time at which the reflectance or transmittance must be calculated.

Following the same development as above, it can be shown that, where three dipoles are retained, the reflectance R(ρ,d,t) is given by

The spatial integral analogous to Eq. (11) can also be computed, and we have

The transmittance T(ρ,d,t) calculated by retaining four dipoles is

and the spatially integrated transmittance is just

A quantity of interest in in vivo transmission spectroscopy is the mean distance traversed by photons before exit.[19] An analytical expression for this quantity can be derived from Eq. (15) for the slab geometry. The mean path length 〈ct〉 is given by

where the integrals have been evaluated using the solutions of Gradsteyn and Ryzhik.[21]

III. Materials and Methods

The predictions of the diffusion model were compared with data from three sources: experimental results reported by Chance et al.,[18] Monte Carlo simulation performed by the authors using a code developed by Flock et al.,[22] and Monte Carlo results published by Delpy et al.[19]

The experimental arrangement has been described in detail elsewhere.[18] Briefly, a pulsed dye laser generated 6-ps pulses at 760 nm, which were incident on the calf of a human volunteer. Scattered light was collected at a distance of 40 mm from the source using a 3-mm diam fiber optic bundle in contact with the skin. The detector was a two-stage microchannel plate tube, and the pulse shape was recorded using standard time resolved photon counting methods. The time resolution of the system, as determined by the full width at half-maximum of the directly detected pulse, was 150 ps.

Results for spatially integrated time resolved reflectance were obtained using a Monte Carlo code described by Flock et al.[22] This program allows anisotropic light scattering where the angular dependence of the scattering is described by the analytical Henyey-Greenstein phase function.[23] The cumulative path length traveled by each photon before escape is scored so that the time resolved reflectance can be examined.

For time resolved transmittance we used the results of Delpy et al.[19] Their Monte Carlo code computes the spatially integrated transmittance through a slab. The program incorporates an experimentally determined phase function measured for rat brain tissue at 783 nm.

IV. Results and Discussion

Figure 2 shows the logarithm of the signal detected in vivo vs time. The optical properties of the muscle could be obtained by fitting the data to Eq. (7), but Eqs. (9) and (10) allow a simple determination of the absorption and transport scattering coefficients. From the asymptotic slope we calculated that μ_a = 0.0176 mm⁻¹, and, using this value and the time of maximum signal, Eq. (10) gave the value (1 − g)μ_s = 0.85 mm⁻¹. These results are in good agreement with the values of μ_a = 0.023 ± 0.004 and (1 − g)μ_s = 0.85 ± 0.08 published by Wilson et al.[24] for excised bovine muscle. Given that our experiment was performed in vivo and that two different species were involved, this agreement may be fortuitous and further studies are clearly required. Nonetheless, this preliminary result shows that the derived coefficients for muscle tissue are reasonable. Using the calculated values, Eq. (7) was used to generate the complete pulse shape, and, as shown in Fig. 2, the agreement with experiment is quite good.

A comparison of Monte Carlo simulations of spatially integrated time resolved reflectance with Eq. (11) is presented in Fig. 3. The tissue was a semi-infinite slab with the same properties as those deduced for the calf muscle discussed above. In the simulation we used a value of g = 0.95 as reported for bovine muscle at 630 nm by Flock et al.[25] The diffusion model predictions have been normalized to the Monte Carlo data at 50 ps, and the agreement between the two is excellent out to 550 ps where the statistical uncertainties in the Monte Carlo data become large. At times shorter than 25 ps there is considerable discrepancy between the two. This is not surprising as photons exiting at these early times have not undergone multiple large angle scattering and their behavior is inadequately described by diffusion theory. Because of the uncertainty due to the finite number of photon histories the data are only shown out to 550 ps. At this time the asymptotic slope −μ_ac has not been attained. Substitution in Eq. (11) shows that at 4 ns the slope would be within 10 % of its terminal value for this particular sample. While μ_a can thus be readily determined from spatially integrated measurements, a close inspection of Eq. (11) shows that the shape of such curves is relatively insensitive to the scattering properties of the medium. Spatially resolved measurements, such as that presented in Fig. 2, are superior if information about (1 − g)μ_s is desired. Of course, if μ_a is determined from a time resolved measurement, (1 − g)μ_s could be calculated from some other noninvasive probe, such as the total (spatially integrated, temporally integrated) reflectance which in the diffusion limit depends only on μ_a/(1 − g)μ_s and the index of refraction.[26]

Also shown in Fig. 3 is the reflectance predicted from Eq. (12) for a slab of identical tissue 10 mm thick along with the results of the Monte Carlo simulation for this geometry. As expected, the loss of photons from the back surface causes the reflectance observed at the front surface to decrease more rapidly with time, although the diffusion model predicts a larger effect than observed in the Monte Carlo simulation. This simple illustration shows that finite tissue geometries can have a significant effect on the observed signal if the observation time is long enough. Caution must, therefore, be exercised in applying simple relations like Eq. (9) to observations made on smaller tissue volumes such as the organs of small animals.

Finally, in Fig. 4 we show that the simple diffusion model is also quite successful in predicting the total time resolved transmittance through a tissue slab as calculated by a Monte Carlo simulation of photon transport. The Monte Carlo results were taken from Delpy et al.[19] and apply to a 10-mm slab of tissue for which μ_a = 0.00434 mm⁻¹ and μ_s = 6 mm⁻¹. An experimental phase function[27] which has an anisotropy parameter, g = 0.72 was used in their simulation. A parameter of considerable interest is the mean path length traveled by photons before leaving the rear surface. From the results of Delpy et al.[19] shown in Fig. 4, we calculated this to be 80.6 mm. The analytic expression of Eq. (16) gives a result of 77.6 mm, an error of only 3.7%. The shape of the transmitted pulse is strongly influenced not only by μ_a but also by (1 − g)μ_s. Such integrated measurements thus have the potential for deducing both optical coefficients.

V. Conclusions

The purpose of this paper has not been to test exhaustively the capability of time resolved reflectance and transmittance measurements for the deduction of tissue optical properties. Rather, the potential of this novel technique has been illustrated by Monte Carlo simulations and preliminary in vivo measurements. A simple model based on the diffusion approximation to the radiative transfer equation has been used to derive analytic expressions for the temporal and spatial dependence of diffusely transmitted and reflected light. These expressions have been shown to predict the detected pulse shape quite accurately for optical properties typical of tissue in the so-called therapeutic window of 650–1300 nm. A limitation of the model is the remaining ambiguity between g and μ_s, as only the product (1 − g)μ_s can be obtained. Some of this ambiguity might be resolved by examining the behavior of the reflectance at very early times where we have shown diffusion theory to be inadequate. The technical demands of such measurements are higher, but time resolution of the order of a picosecond is within current capabilities. The sort of measurements described in this paper require resolution of the order of 100 ps where off-the-shelf hardware is quite adequate. At some wavelengths, practical systems could incorporate pulsed laser diodes which would greatly reduce the cost and complexity of the light source. Further modeling and experimentation are also required to study the influence of more realistic geometries on the observed signals.

The full potential of time resolved measurements has not been completely explored. Besides the possibility of picosecond work, the interesting problem of spatial location of an optical inhomogeneity in, say, the brain or breast may be aided by studying the time behavior of scattered light. Here numerical solutions of the diffusion equation may provide an alternative to time-consuming Monte Carlo simulations. Much more work needs to be done to explore fully the dimension added by time resolution, but we have shown that such measurements are already capable of noninvasively probing the optical properties of living tissue.

The authors are grateful to Stephen Flock for providing the Monte Carlo results. This work was financially supported by the National Cancer Institute of Canada.

Michael Patterson and B. C. Wilson also work in the Departments of Radiology and Physics of McMaster University.

Figures

 

Fig. 1 (a) Geometry for the calculation of R(ρ,t) for a semi-infinite homogeneous medium. The incident pencil beam is assumed to create an isotropic photon source at depth z₀, indicated by the filled circle. The boundary condition ϕ(ρ,0,t) = 0 can be met by adding a negative source indicated by the open circle. (b) Geometry for the calculation of the time resolved reflectance and transmittance from a homogeneous slab. The boundary conditions ϕ(ρ,0,t) = 0 and ϕ(ρ,d,t) = 0 can be met by adding an infinite series of dipole photon sources. The first four are shown in this illustration.

Download Full Size | PDF

 

Fig. 2 Experimentally recorded reflectance signal from the calf muscle of a human volunteer at 760 nm is indicated by the solid line. The detector fiber bundle was positioned 40 mm from the pulsed laser source. The open circles show the pulse shape predicted by the diffusion model expression of Eq. (7) normalized to the peak value.

Download Full Size | PDF

 

Fig. 3 Monte Carlo (symbols) and diffusion model (smooth curves) calculation of R(t), the spatially integrated time resolved reflectance, from a semi-infinite tissue slab and a slab 10 mm thick. The optical properties of the slab are typical of human skeletal muscle at 760 nm. The diffusion model results have been matched to the Monte Carlo data at 50 ps. Note that the units of R(t) are ps⁻¹ as R(t) is the number of photons reaching the surface per unit time per incident photon.

Download Full Size | PDF

 

Fig. 4 Comparison of the predictions of the diffusion model (solid line) and Monte Carlo (dashed line) results of Delpy et al.[19] for the spatially integrated time-resolved transmittance through a 10-mm tissue slab. Each curve has been normalized at its peak value.

Download Full Size | PDF

References

1. M. Tamura, O. Hazeki, S. Nioka, B. Chance, and D. Smith, “Simultaneous Measurements of Tissue Oxygen Concentration and Energy State by Near-Infrared and Nuclear Magnetic Resonance Spectroscopy,” in Chemoreceptors and Reflexes in Breathing, S. Lahiri, Ed. (Oxford U.P., New York, in press 1989).

2. M. Cope and D. T. Delpy, “System for Long-Term Measurement of Cerebral Blood and Tissue Oxygenation on Newborn Infants by Near Infrared Transillumination,” Med. Biol. Eng. Comput. 26, 289–294 (1988).

3. B. C. Wilson and M. S. Patterson, “The Physics of Photodynamic Therapy,” Phys. Med. Biol. 31, 327–360 (1986).

4. S. L. Jacques and S. A. Prahl, “Modeling Optical and Thermal Distributions in Tissue During Laser Irradiation,” Lasers Surg. Med. 6, 494–503 (1987).

5. F. F. Jobsis, “Noninvasive, Infrared Monitoring of Cerebral and Myocardial Oxygen Sufficiency and Circulatory Parameters,” Science 198, 1264–1267 (1977).

6. R. A. J. Groenhuis, J. J. Ten Bosch, and H. A. Ferwerda, “Scattering and Absorption of Turbid Materials Determined from Reflection Measurements. 2: Measuring Method and Calibration,” Appl. Opt. 22, 2463–2467 (1983).

7. R. F. Bonner, R. Nossal, S. Havlin, and G. H. Weiss, “Model for Photon Migration in Turbid Biological Media,” J. Opt. Soc. Am. A 4, 423–432 (1987).

8. J. M. Steinke and A. P. Shepherd, “Diffuse Reflectance of Whole Blood: Model for a Diverging Light Beam,” IEEE Trans. Biomed. Eng. BME-34, 826–834 (1987).

9. J. A. Weinman and S. T. Shipley, “Effects of Multiple Scattering on Laser Pulses Transmitted Through Clouds,” J. Geophys. Res. 77, 7123–7128 (1972).

10. A. Ishimaru, “Diffusion of a Pulse in Densely Distributed Scattered,” J. Opt. Soc. Am. 68, 1045–1050 (1978).

11. K. Furutsu, “On the Diffusion Equation Derived from the Space–Time Transport Equation,” J. Opt. Soc. Am. 70, 360–366 (1980).

12. K. Shimizu, A. Ishimaru, L. Reynolds, and A. P. Bruchner, “Backscattering of a Picosecond Pulse from Densely Distributed Scatterers,” Appl. Opt. 18, 3484–3488 (1979).

13. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic, New York, 1978), pp. 175–190.

14. S. Chandrasekhar, “Stochastic Problems in Physics and Astronomy,” Rev. Mod. Phys. 15, 1–88 (1943).

15. J. J. Duderstadt and L. J. Hamilton, Nuclear Reactor Analysis (Wiley, New York, 1976), pp. 140–144.

16. R. P. Hemenger, “Optical Properties of Turbid Media with Specularly Reflecting Boundaries: Applications to Biological Problems,” Appl. Opt. 16, 2007–2012 (1977).

17. G. Eason, A. Veitch, R. Nisbet, and F. Turnbull, “The Theory of the Backscattering of Light by Blood,” J. Phys. D 11, 1463–1479 (1978).

18. B. Chance, S. Nioka, J. Kent, K. McCully, M. Fountain, R. Greenfeld, and G. Holtom, “Time Resolved Spectroscopy of Hemoglobin and Myoglobin in Resting and Ischemic Muscle,” Anal. Biochem. 174, 698–707 (1988).

19. D. T. Delpy, M. Cope, P. van der Zee, S. Arridge, S. Wray, and J. Wyatt, “Estimation of Optical Pathlength Through Tissue from Direct Time-of-Flight Measurements,” Phys. Med. Biol. 33, 1433–1442 (1988).

20. S. Ito and K. Furutsu, “Theory of Light Pulse Propagation Through Thick Clouds,” J. Opt. Soc. Am. 70, 366–374 (1980).

21. I. S. Gradsteyn and I. M. Ryzhik, Table of Integrals, Series and Products (Academic, New York, 1980), p. 340.

22. S. T. Flock, M. S. Patterson, and B. C. Wilson, “Monte Carlo Modelling of Light Propagation in Highly Scattering Tissues. I. Model Predictions and Comparison with Diffusion Theory,” IEEE Trans. Biomed. Eng.1989, in press.

23. L. G. Henyey and J. L. Greenstein, “Diffuse Radiation in the Galaxy,” Astrophys. J. 93, 70–83 (1941).

24. B. C. Wilson, M. S. Patterson, S. T. Flock, and J. D. Moulton, “The Optical Absorption and Scattering Properties of Tissues in the Visible and Near-Infrared Wavelength Range,” in Light in Biology and Medicine, Vol. 1R. H. Douglas, J. Moan, and F. Doll’Acqua eds. (Plenum, New York, 1988) pp. 45–52.

25. S. T. Flock, B. C. Wilson, and M. S. Patterson, “Total Attenuation Coefficients and Scattering Phase Functions of Tissues and Phantom Materials at 633 nm,” Med. Phys. 14, 835–841 (1987).

26. B. C. Wilson and M. S. Patterson, “The Determination of Light Fluence Distributions in Photodynamic Therapy,” in Photodynamic Therapy of Neoplastic Disease, D. Kessel, Ed. (CRC Press, Cleveland, 1989), in press.

27. P. van der Zee and D. T. Delpy, “Computed Point Spread Functions for Light in Tissue Using a Measured Volume Scattering Function,” Adv. Exp. Med. Biol. 222, 191–197 (1988).