Study of photon migration with various source-detector separations in near-infrared spectroscopic brain imaging based on three-dimensional Monte Carlo modeling

We have simulated the photon migration with various sourcedetector separations based on a three-dimensional Monte Carlo code. The whole brain MRI structure images are introduced in the simulation, and the brain model is more accurate then the previous studies. The brain model consists of scalp, skull, CSF layer, gray matter, and white matter. We demonstrate dynamic propagating movies under different source-detector separations. The multiple backscattered intensity from every layer of the brain model is obtained by marking the deepest layer which every photon can reach. Also the influences of an absorption target on brain cortex are revealed. ©2005 Optical Society of America OCIS codes: 170.3660 (Light Propagation Through Tissue); 170.5280 (Photon Migration) References and Links 1. G. Boas, "Noninvasive Imaging of the Brain," Optics News 15, 52-55 (2004) 2. A. Villringer and B. Chance, "Non-invasive optical spectroscopy and imaging of human brain function," Trends Neurosci. 20, 435-442 (1997) 3. A. Villringer, J. Planck, C. Hock, L. schleinkofer, U. Dirnagl, "Near infrared spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of brain function in human adults," Neurosci. Lett. 154, 101-104 (1993) 4. Y. Fukui, Y. Ajichi, E. Okada, "Monte Carlo prediction of near-infrared light propagation in realistic adult and neonatal head models," Appl. Opt. 42, 2881-2887 (2003) 5. H. Koizumi, T. Yamamoto, A. Maki, Y. Yamashita, H. Sato, H. Kawaguchi, N. Ichikawa, "Optical topography: practical problems and new applications," Appl. Opt. 42, 3054-3062 (2003) 6. E. Okada, D.T. Delpy "Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer" Appl. Opt. 42, 2906-2914 (2003) 7. T. Hayashi, Y. Kashio, E. Okada, "Hybrid Monte Carlo-diffusion method for light propagation in tissue with a low-scattering region," Appl. Opt. 42, 2888-2896 (2003) 8. A. Ishimaru, Wave propagation and scattering in random media, I and II (Academic, New York, 1978) 9. A. Yodh and B. Chance, "Spectroscopy and imaging with diffusing light," Phys. Today 48, 38-40 (1995) 10. S.K. Gayen and R.R. Alfano, "Emerging optical biomedical imaging techniques," Opt. Photon. News 7, 1722 (1996) 11. S. R. Arridge, "Optical tomography in medical imaging," Inverse Problems 15, R41-R93 (1999) 12. P. Bruscaglioni, G. Zaccanti, and Q. Wei, "Transmission of a pulsed polarized light beam through thick turbid media: numerical results," Appl. Opt. 32, 6142-6150 (1993) 13. M. J. Rakovic, G.W. Kattawar, M. Mehrbeolu, B.D. Cameron, L. V. Wang, S. Rastegar, and G. L. Cote, "Light Backscattering Polarization Patterns from Turbid Media: Theory and Experiment," Appl. Opt. 38, 3399-3408 (1999) 14. S. Bartel and A. H. Hielscher, "Monte Carlo Simulations of the Diffuse Backscattering Mueller Matrix for Highly Scattering Media," Appl. Opt. 39, 1580-1588 (2000) (C) 2005 OSA 17 October 2005 / Vol. 13, No. 21 / OPTICS EXPRESS 8339 #8338 $15.00 US Received 1 August 2005; revised 13 September 2005; accepted 26 September 2005 15. M. Moscoso, J. B. Keller, and G. Papanicolaou, "Depolarization and blurring of optical images by biological tissue," J. Opt. Soc. Am. A 18, 948-960 (2001) 16. H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, "Monte Carlo and Multicomponent Approximation Methods for Vector Radiative Transfer by use of Effective Mueller Matrix Calculations," Appl. Opt. 40, 400-412 (2001) 17. B. Kaplan, G. Ledanois, and B. villon, "Mueller Matrix of Dense Polystyrene Latex Sphere Suspensions: Measurements and Monte Carlo Simulation," Appl. Opt. 40, 2769-2777 (2001) 18. X. Wang and L. V. Wang, "Propagation of polarized light in birefringent turbid media: A Monte Carlo study," J. Biomed. Opt. 7, 279-290 (2002) 19. I. Lux and L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations (CRC Press, Boca Ration, Fla., 1991) 20. G. W. Kattawar and G. N. Plass, "Radiance and polarization of multiple scattered light from haze and clouds," Appl. Opt. 7, 1519-1527 (1968) 21. E. Okada, D. T. Delpy "Near-infrared light propagation in an adult head model. II. Effect of superficial tissue thickness on the sensitivity of the near-infrared spectroscopy signal" Appl. Opt. 42, 2915-2922 (2003) 22. D. A. Boas, J. P. Culver, J. J. Stott, A. K. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10, 159170 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-3-159 23. C. F. Bohren and D. R. Huffman, “Absorption and Scattering of Light by Small Particles,” John Wiley & Sons, 1983 24. D. Haensse, P. Szabo, D. Brown, J. Fauchere, P. Niederer, H. Bucher, M. Wolf, “New multichannel near infrared spectrophotometry system for functional studies of the brain in adults and neonates,” Opt. Express 13, 4525-4538(2005), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-12-4525 25. M. A. Franceschini, D. A. Boas, “Noninvasive measurement of neuronal activity with near-infrared optical imaging,” NeuroImage 21, 372-386 (2004) 26. J. C. Ramella-Roman, S. A. Prahl, S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: Part I,” Opt. Express 13, 4420-4438 (2005), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-12-4420 27. X. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: A Monte Carlo study,” J. Biomed. Opt. 7, 279-290 (2002) 28. L. H. Wang, S. L. Jacques, L-Q Zheng, “MCML – Monte Carlo modeling of photon transport in multilayered tissues,” Computer Methods and Programs in Biomedicine 47,131-146 (1995).


Introduction
Generally, brain imaging techniques such as conventional x-ray radiography, magnetic resonance imaging (MRI), computed tomography (CT), and positron emission tomography (PET) have allowed the noninvasive investigation in the human brain.Besides, optical imaging of the brain is a rapidly growing field because it has attracted considerable interest recently due to a number of theoretical advantages in comparison with other brain imaging modalities [1][2][3].Light in the near-infrared range is only moderately absorbed by water, hemoglobin, and other significant body substances and thus can penetrate several centimeters inside the human brain.The near-infrared spectroscopy (NIRS) is such a technique that measures changes in oxy-hemoglobin and deoxy-hemoglobin concentrations with different absorption coefficients of hemoglobin at various wavelengths.Since absorption of the activated region varies by changes in the blood volume and oxygenation, one can measure brain activation by detecting the intensity change of near-infrared light that passes through the brain.Due to its high time resolution (~ms) and good spatial resolution (~cm) offers the option of continuously recording changes of physiological parameters during brain activation.Therefore, the NIRS technique is widely used for functional brain imaging in neurosurgery, neurology, psychiatrics, pediatrics, rehabilitation, cognitive neuroscience, linguistics, and etc [1,4].
Several significant advantages of NIRS brain imaging include high biochemical specificity, compact system setup, and relatively low cost.Near-infrared light in the wavelength range from 650 to 950 nm can penetrate the adult skull and propagate while in the brain.With brain activation, the oxyhemoglobin pours into a specific region that corresponds to neuron activation.The different oxygenation indicates the changes in brain function with neurovascular coupling processes.Based on the modified Beer-Lambert's Law (MBL), the brain activation of a specific region can be measured from the detected light signals.In order to improve the signal-to-noise ratio in NIRS systems, some approaches were adopted, such as varying of the source-detector separation and the time-gating method in previous studies [5][6][7].
As mentioned above, the information regarding the brain was associated with detected optical signals that strongly depend on the separation distance between the source and the detector.In this study, we observed the photon migration with multiple scattering events in a three-dimensional head model of a human based on Monte Carlo simulations.10 8 photons were emitted into the head model and propagated within 10 hours' computing time on a Pentium IV 3.2 GHz CPU.The detectors were located on the surface of the scalp.The received photons of the detectors were traced with various source-detector separations.In this paper, the new contributions are stated as follows: 1.In previous study, the results of Monte-Carlo simulations were usually based on the semiinfinity five-layer structure [6] or two-dimensional head model with a MRI scan slice [4].
In our study, we simulated of photon migration within the human brain on the basis of realistic three-dimensional adult head model obtained by means of an anatomical MRI system.2. Based on our algorithm, the dynamic photon migration analysis can be described in threedimensional model.Although many papers discuss the photon migration in human brain, as our knowledge, it is a first paper to show the dynamic photons traveling form source to several detectors, respectively.3.In this paper, we used determination of contributions of the photons scattered from the different regions within the head model as a function of separation distance between the light source and the detectors locations.Because the three-dimensional model approaches the realistic human head, the signal-to-noise ratio evaluation and optimal choice of source-detector separation may provide more helpful information for NIRS systems design.4. One of the goals of the work is to study the sensitivity of the NIRS measurement to functional hemodynamic signals generated in the gray and white matters which are considered as useful information in the NIRS.For this study, we set an abnormality on the surface of gray matter in simulation.The received intensity distribution on the surface of head illustrates the locations of the target and the relative response.This paper is organized as follows: In Section 2, we describe the Monte Carlo algorithm, the structure of three-dimensional brain model, and optical parameters of brain tissues.Simulation results are discussed in Section 3. Finally, discusses are drawn in Section 4.

Methods
The behavior of photon migration in turbid media is a fundamental research to many practical applications including imaging of colloidal suspensions and biological tissues [8][9][10][11].Because the analytical solution of radiative transport equation within a bounded turbid medium is complicated, Monte Carlo algorithms of light propagation through turbid media have been pursued extensively and applied to characterize particle suspensions and biological tissues [12][13][14][15][16][17][18][19][20][26][27][28].Based on the method, some research groups defined the differential pathlength factor (DPF), partial differential path-length factor (PDPF) and spatial sensitivity profile (SSP) to indicate the sensitivity of the NIRS signal with absorption changes in the tissues [4,[6][7].Also, the effect of the thickness of the skull and CSF layer on light propagation in an adult head model was investigated [21].Three-dimensional Monte Carlo code has been described to calculate the photon migration through highly scattering media with varying optical properties and arbitrary boundary conditions [22].

Monte Carlo algorithm
Monte Carlo modeling offers a flexible approach to analyze light transport in turbid media.Our method is based on a set of rules that describe the characteristics of a photon behavior in turbid tissues.The propagation behavior of a photon in tissues can be decided with two parameters: (1) the mean free path of a scattering or absorption event, and (2) the scattering angle of a scattering event.In addition, Snell's law and Fresnel reflection formulas were applied at each boundary.When a photon is encounter with a scattering event, a new step size, deflection angle, and azimuthal angle will be assigned.The formal solution, Mie theory, describes absorption and/or scattering event with a sphere that has been available in previously study [23].After a sphere was illuminated by a beam of light with specified characteristics, scattered photon propagates in all directions.Then, the photon will be scattered by the other scatters again and again.It seems an intractable task to trace all the scattered waves in all directions simultaneously when two or more scatters are arranged randomly.However, we adopted the "photon-appearance" light with proper energy weighting for simulating of wave behaviors.
In our simulation, before the photons start to emit, we created a matrix to record the spatial arrangement of tissue constituents.Each voxel in the matrix is assigned an integer to represent the type of tissue in this voxel (e.g. 1 = air, 2 = scalp, etc.).And the scattering properties governed by the Mueller matrix are calculated with the optical parameters in Table 1.The initial position and direction of photon were defined in each photon emission.In our case, the point source was used that means all the photons start to emit at the same direction.Then the pass length was determined.The free path selection of the photon is illustrated in Eq. ( 2).Besides, if the photon propagates through a voxel boundary, the step size should be modified as the way in Eq. ( 3).If the refraction index is different between two adjacent voxels, the effect of refraction or reflection may occur alternatively.During the propagation, all the voxels that the photon has passed through are recorded.Therefore, the movie of photon migration in human brain can be made.
By the definition of scattering, absorption, and extinction coefficients (μ s, μ a, and μ t ), we have known that μ t Δs denotes the probability of interaction with the turbid medium when a photon travels an infinitesimal distance Δs.The probability density function of free path s 1 is: (1) We can obtain a sampled value s 1 : where ξ is a sampled value of a uniform random variable within the interval [0,1].Therefore, Eq. ( 2) was used to sample the step size of a photon movement.However, before the photon is scattered by a particle, the photon may pass through the boundary of different voxels.So the decision of the step size of photon should be modified.In this condition, Eq. ( 2) should be: where i and μ ti are the index and the extinction coefficient of the voxel, respectively.s i is the path length of photon in the i voxel.Some concepts were concerned with μ t should be remarked again.The value of μ t represents the degree of scattering and absorption when light transports in a turbid medium.Intuitively, if the concentration of the medium is dense, μ t will be large.Besides the concentration of the medium, the other parameters, including the particle size, incident wavelength, material inside the particles and the composition of solvent, etc., will also influence the value of μ t as well.In other words, μ t is a macroscopic parameter of a turbid medium.Total attenuation coefficient, μ t , not only abate the intensity of light but also decrease the polarization degree of light.

Brain Model
A three-dimensional image of brain structure contains 256x256x50 voxels was adopted in simulations.Each voxel indicates a 1x1x1 mm 3 cube in size.The brain consists of five layers including scalp, skull, cerebrospinal fluid layer (CSF layer), gray matter, and white matter.For the wavelength of 800 nm, the scattering coefficients of the five tissues are 1.9, 1.6, 0.24, 2.2, and 9.1 mm -1 , respectively.Also, the absorption coefficients of the five tissues are 0.018, 0.016, 0.004, 0.036, and 0.014 mm -1 , respectively [6].The optical properties of each layer of brain that were used in the simulation are shown in Table 1.First, we calculated the Mueller matrix of the five layers and save the row data in the flash memory.Then, based on the threedimensional MRI image, we assign each voxel as air, scalp, skull, CSF layer, gray matter, and white matter.The numbers of voxels, the size of cube, the type of tissue, the light source position and the detector positions were adjusted to various simulations.Here, the voxel size was chosen as 1x1x1 mm 3 for fitting of the MRI resolution.The left figures of Fig. 1 indicate tomograms with different depths of the brain.The 50 two-dimensional MRI images were used from the top head to down, and each image contains 256x256 pixels.Thus, each voxel of the space in our simulations can be well defined.Also, the right figure of Fig. 1 illustrates the structure of five brain layers.

Results
In our study, the region of visual cortex was monitored with migrated photons in Monte Carlo simulations.The light source was located at the 2.2 cm from the top of rear head.In Fig. 2(a), the movie shows the migration of all the 10 8 photons in the horizontal cross section with excitation.The intensity of light was revealed with pseudo-color mapping in the movie.The energy propagation of light diffuse with temporal evolution.In the movie, we can see the penetration depth of light is about 4 cm that reach the white matter after 400 to 600 psec light propagation.Besides, we can observe that a part of photons was guided along the CSF layer in the movie.

Scalp Skull CSF Gray matter White matter
Figure 2 (b) demonstrates the photon migration in the vertical cross section at the center at the same time with Fig. 2(a).The penetration depth of light is the same with Fig. 2(a).Again, photons propagate longer distance along the CSF layer can be observed clearly.
In practical situation of NIRS experiments, there are two questions we interested in this study.First, how many photons can be backscattered outward through the scalp and be collected by the detectors?Furthermore, how many received photons were backscattered from the gray matter and the white matter in a NIRS measurement?Figure 3 demonstrates the paths of detected photons with various source-detector separations.The movies show the trajectories of the photons received by the detectors in the horizontal cross section.Figure 3 (a), (b), and (c) are the cases of source-detector separation with 1 cm, 2 cm, and 3 cm, respectively.The yellow arrows indicate the light source positions, and the red arrows represent the detector positions in the movies.In the case of the separation with 1 cm, we can observe that the most part of detected photons were scattered only in the skull.Because the functional hemodynamic signals generate in the gray and white matters, the backscattered photons from scalp and skull are a kind of background signals.In the case of separation with 2 cm, the received intensity of light becomes lower from the skull layer than that with the case of 1 cm.Also, the received photons reached deeper in the brain.In the case of separation with 3 cm, the received intensity of light decrease more strongly.However, the ratio of received photons that are multiple-backscattered from the gray and white matters is increased; i.e., the signal-to-noise ratio of NIRS is improved.Thereafter, we located 10 detectors away from the light source with the arc length of 1 to 10 cm against the head in simulation.In Fig. 4, the curve of intensity distribution shows the received intensity was decreasing with distance of source-detector separation increasing.Because we recorded all the paths of the received photons in the simulations, the visited layers of each photon were marked.Figure 5 shows the ratios of the backscattered intensities from different layer versus the source-detector separation.Obviously, the ratio of skull layer decreases rapidly with the distance of source-detector separation increases.On the contrary, the ratio of white matter layer increases rapidly with the same condition.However, the total received intensity was decreasing with the separation increasing.Hence, it is a trade-off problem between the received intensity and the useful information in NIRS measurement.
The NIRS technique may provide information about an abnormality in a brain.To observe the abnormality of an absorbing target in a human brain, a target zone of 15 mm 3 in size on the surface of gray matter was defined in simulation (as shown in Fig. 6).The absorption coefficients of the target zone were set in amounts varying between 0.36 and 7.2 cm -1 .The values of the absorption coefficient may be exaggerated in this preliminary study.However, the simplified method can provide a variation trend for understanding.

Discussions
In this work, the simulation results presented demonstrate the photon migration with different source-detector separations in a head model.While the literature on the analysis about sourcedetector separation are voluminous, few previous efforts report dynamic propagating analysis combined with different source-detector separation.As mentioned before, we described the Monte Carlo method that is capable of performing source-detector separation analysis with dynamic photon migration movies.The transmitted photon path, so-called "banana shape," is showed in the movies of Fig. 3 The simulations we have performed in a realistic head model suggest that in the case of the adult human head, the penetration depth of the NIR light is several centimeters into the brain.The result of Fig. 1(a) is similar to that reported in [22], but we obtained a result of the vertical section in Fig. 1(b) to verify that our three-dimension head model is more realistic.
Figure 5 shows the ratio of multiple backscattered intensities of different layers for the brain model as a function of source-detector spacing.Our results show that the summed ratio of multiple backscattered light from the gray and white matter layer is greater than 50%, while the source-detector separation exceeds 2 cm and the summed ratio is over 90% for separations larger than 4 cm.However, in Fig. 4 we observed that the received intensity of the sourcedetector separation of 4 cm is roughly ten times smaller than that of 2 cm.In other words, the detected light carries more information of the cortex at lower intensity when the sourcedetector separation increases.That is the core of the trade-off problem in the NIRS measurement.The source-detector separation is usually chosen between 2 and 3 cm [3][4][5][6][7].The dynamic photon transportation with generally-used separations was also provided in our simulation results.
Figure 7 demonstrates the variation of surface intensity with an absorption target in the brain cortex.The absorption coefficient of the target region was increased gradually in order to observe the affected signals with the abnormality.We set the absorption coefficient of the target exaggeratedly large (7.2 cm -1 ) to make observation of the responsive region beyond the target region easier.It implies that the Monte Carlo simulation may offer a helpful method to the researchers for design of the configuration of the sources and detectors, and then for improvement of the signal-to-noise ratio in NIRS measurements.
Up to now, most Monte Carlo codes have been suited to a single wavelength.However, the NIRS system usually applies multi-wavelength sources, e.g., 730, 770, and 805 nm in Ref. [24], and 690 and 830 in Ref. [25].Therefore, we have reformed the Monte Carlo code for multi-wavelength sources to approach the practical NIRS measurement.

Fig. 1 .
Fig. 1.The nine MRI images on the left are corresponding to different depths of the brain.Here, we took 50 images to fill all the voxels.Each MRI image consists of 256x256 pixels.The schematic diagram on the right shows the anatomical structure of the human head.

Fig. 2 .
Fig. 2. The movies show the migration of a pulse of light with 800 nm wavelength through a brain structure.Movies are given for (a) the horizontal cross section, (b) the vertical cross section of human head.The sizes of MOV movie files are 2.08 and 1.1 Mega-Bytes, respectively.

Fig. 3 .
Fig. 3.The movies show the photon migration of the received photons with different distances of source-detector separation (a) 1 cm (b) 2 cm (c) 3 cm in the horizontal cross section of human head.The sizes of MOV movie files are 0.54, 0.51 and 0.49 Mega-Bytes, respectively.

Fig. 5 .
Fig. 5.The distributions of ratio of the received intensity from different layers of brain versus the distance of source-detector separation.

Fig. 6 . 1 Figure 7
Fig. 6.Target zone with the absorption coefficients at 0.36 (background), 0.72, 1.08, 1.44, 1.8, 3.6, 7.2 cm -1Figure7shows the received intensity distribution on the surface of the head.The red arrow indicates the position of the light source.The target zone is located 20 mm beneath the red circle.The received intensities between the source and the target do not change significantly.However, the intensity distribution beyond the target decreases when the absorption coefficient of the target increases.

Fig. 7 .
Fig. 7.The movie of received intensity distribution on the surface of head (The size of MOV movie files is 0.12 Mega-Bytes.) (a), 3(b), and 3(c).

Table 1 .
Optical properties of each tissue